Я пытаюсь сделать учебный пример для демонстрации предположений о причинно-следственной связи, используя различные вероятности, особенности чувствительности, неправильную классификацию и т. Д. c. Тем не менее, я не настолько знаком с r и делаю такие процедуры. Любая помощь приветствуется:
Имеет
#constants
n_combos <- 216 # all possible combinations of the parameters we will set below
num_sims <- 1 # number of simulations (for now)
num_measures <- 7 # measures of interest
n <- 10^6 # samples size
trials <- 1 # data generated treated as 1 trial
matrix_sim_output216 <- matrix(NA, ncol=num_measures, nrow=n_combos) # This creates the master matrix
# starts at 0 rows because all sims will concatenate onto it - matrix of (216*10^4)x10; n_combos*num_sims x num_measures
set.seed(456) # set the PRNG so this can be reproduced
# vectors with chosen parameter values for variable relationships
v_X_Y=c(2) # column 1, relationship of X and Y, always stays the same
v_C_Y=c(0.5, 1, 2) # column 2, C and Y
v_C_X=c(0.5, 1, 2) # column 3, C and X
v_Prev_C=c(0.1, 0.2, 0.5) # column 4, Prevalence of C
v_se=c(0.8, 0.6) # column 5, sensitivity
v_sp=c(0.9, 0.6) # column 6, specificity
v_intercept=c(-1, -4) # column 7
# create dataset containing all possible combinations of parameters chosen
df_combos <- expand.grid(
v_X_Y,
v_C_Y,
v_C_X,
v_Prev_C,
v_se,
v_sp,
v_intercept,
KEEP.OUT.ATTRS = FALSE
)
# check math on this, 3 vars with 3 values & 3 vars with 2 values
# nCr; 3C1*3C1*3C1*2C1*2C1*2C1 = [(3C1)^3 * (2C1)^3] = 216
dim(df_combos) # 216x7
# order parameters by largest in each column so once results are obtained, we can rbind() to show what the parameters were that created the results for each individual row i.e. columns of the variables, where every row has a different value
+-------+-------+-------+----------+------+------+-------------+
| v_X_Y | v_C_Y | v_C_X | v_Prev_C | v_se | v_sp | v_intercept |
+-------+-------+-------+----------+------+------+-------------+
| 2 | 2 | 2 | 0.5 | 0.8 | 0.9 | -4 |
| 2 | 2 | 2 | 0.5 | 0.8 | 0.9 | -1 |
| . | . | . | . | . | . | . |
+-------+-------+-------+----------+------+------+-------------+
df_combos <- df_combos[order(
df_combos[,1], # v_X_Y
df_combos[,2], # v_C_Y
df_combos[,3], # v_C_X
df_combos[,4], # v_Prev_C
df_combos[,5], # v_se
df_combos[,6], # v_sp
df_combos[,7] # v_intercept
), ]
# sorting by column index
# sorted by largest from left to right so we know what the parameters are
# name columns
var.labels <- c("v_X_Y",
"v_C_Y",
"v_C_X",
"v_Prev_C",
"v_se",
"v_sp",
"v_intercept"
)
df_combos <- setNames(df_combos, var.labels)
Теперь я хочу создать переменные, используя эти значения, но я хочу, чтобы они были
- Создано
- Анализировано (показано ниже)
- Сохранено в основной кадр данных
Все по каждой строке и соответствующим значениям столбца. Другими словами, я хочу выполнить один и тот же анализ в 216 разных строках, где мои переменные создаются с использованием значения столбца каждой указанной строки c, где результаты анализа хранятся в кадре данных, который непосредственно соответствует этой строке и ее значения столбца. Например, если (B = бета), используя df_combos сверху,
Создание переменной
for (row in 1:nrow(df_combos)) { # should grab the values from df_combos to create variables by column of row 1 then loop to go through row 2 through row n=216
X_Y <- df_combos[row, "v_X_Y"]
C_Y <- df_combos[row, "v_C_Y"]
C_X <- df_combos[row, "v_C_X"]
Prev_C <- df_combos[row, "v_Prev_C"]
se <- df_combos[row, "v_se"]
sp <- df_combos[row, "v_sp"]
intercept <- df_combos[row, "v_intercept"]
# confounder
C1 <- rbinom(n, trials, Prev_C)
# exposure
X <- rbinom(n, trials, inv.logit(intercept+C_X*C1))
# outcome
Y <- rbinom(n, trials, inv.logit(intercept+C_Y*C1+X_Y*X))
# Notice that X is a vector such that X = [0,1,1,0,....] etc
# we want to multiply these values by values of se and 1-sp
# now we create a vector such that p = [1-sp, se, se, 1-sp,....] etc
# start with sensitivity
p <- matrix(X*df_combos[row, “se”], ncol = 1, nrow = n)
# now 0 is replaced by 1-sp
p[X==0] <- (1-sp)
# misclassification variable, X*=S
S <- rbinom(n=n, trials, p)
# create df of variables to use in GLMs and by-hand lnORs
dset <- data.frame(X, C1, Y, S)
Анализ, сравнение оценок эффекта (бета-версии)
# Estimate 1) X -> Y where C1 exists
B1 <- summary(glm(Y~X+C1, data = dset, family=binomial(link=logit)))$coefficients[2, 1]
# Estimate 2) should be equal to Observed Effect, X* -> Y i.e. log((a/b))/((c/d))
B2 <- summary(glm(Y~S, data = dset, family = binomial(link = logit)))$coefficients[2,1]
# Estimate 3) create 2x2 observed cells (exposure misclassification is uncontrolled)
a = sum(dset$S==1 & dset$Y==1)
b = sum(dset$S==1 & dset$Y==0)
c = sum(dset$S==0 & dset$Y==1)
d = sum(dset$S==0 & dset$Y==0)
N1 = a+c
N0 = b+d
#B2.1 =log(((a/b))/((c/d))) should match glm B2
# adjusted
A = (a-((1-sp)*N1))/(se+sp-1)
B = (b-((1-sp)*N0))/(se+sp-1)
C = N1-A
D = N0-B
B3 <- log(((A/B))/((C/D)))
# Estimate 4) Observed
B4 <- summary(glm(Y~S+C1, data = dset, family = binomial(link = logit)))$coefficients[2,1]
# Estimate 5) Difference
B_diff = abs(X_Y – B4)
# compare prevalence of confounder (C1) among exposure (X) and misclassified exposure (S)
C1_prev_in_X = sum(dset$X==1 & dset$C1==1)/n
C1_prev_in_S = sum(dset$S==1 & dset$C1==1)/n
# create vector for easy concatenation later
matrix_sim_output216[row,1] <- B1
matrix_sim_output216[row,2] <- B2
matrix_sim_output216[row,3] <- B3
matrix_sim_output216[row,4] <- B4
matrix_sim_output216[row,5] <- B_diff
matrix_sim_output216[row,6] <- C1_prev_in_X
matrix_sim_output216[row,7] <- C1_prev_in_S
}
df_output216 <- as.data.frame(matrix_sim_output216)
df_params_output216 <- rbind(df_combos, df_output216)
, где df_params_output216 будет выглядеть следующим образом
Want
+-----+-----+-----+--------+-----+-----+-----------+----+----+----+----+----------+-----------+-----------+
| X_Y | C_Y | C_X | Prev_C | se | sp | intercept | B1 | B2 | B3 | B4 | |X_Y-B4| | C1_Prev_X | C1_Prev_S |
+-----+-----+-----+--------+-----+-----+-----------+----+----+----+----+----------+-----------+-----------+
| 2 | 0.5 | 2 | 0.5 | 0.8 | 0.9 | -1 | . | . | . | . | . | . | . |
| 2 | 0.5 | 2 | 0.5 | 0.8 | 0.9 | -4 | . | . | . | . | . | . | . |
| . | . | . | . | . | . | . | . | . | . | . | . | . | . |
+-----+-----+-----+--------+-----+-----+-----------+----+----+----+----+----------+-----------+-----------+
Но они не совпадают соответствующим образом. Под этим я подразумеваю, что если я запускаю этот код без l oop и просто вставляю эти настройки параметров слева (ниже), проверяя каждую строку отдельно, мои бета-версии не равны тем, что когда я запускаю l oop. Вот так
Запуск l oop:
+-----+-----+-----+--------+-----+-----+-----------+----+----+----+----+----------+-----------+-----------+
| X_Y | C_Y | C_X | Prev_C | se | sp | intercept | B1 | B2 | B3 | B4 | |X_Y-B4| | C1_Prev_X | C1_Prev_S |
+-----+-----+-----+--------+-----+-----+-----------+----+----+----+----+----------+-----------+-----------+
| 2 | 0.5 | 2 | 0.5 | 0.8 | 0.9 | -1 | q | w | e | r | t | a | s |
+-----+-----+-----+--------+-----+-----+-----------+----+----+----+----+----------+-----------+-----------+
По каждой строке нет l oop:
+-----+-----+-----+--------+-----+-----+-----------+----+----+----+----+----------+-----------+-----------+
| X_Y | C_Y | C_X | Prev_C | se | sp | intercept | B1 | B2 | B3 | B4 | |X_Y-B4| | C1_Prev_X | C1_Prev_S |
+-----+-----+-----+--------+-----+-----+-----------+----+----+----+----+----------+-----------+-----------+
| 2 | 0.5 | 2 | 0.5 | 0.8 | 0.9 | -1 | p | l | m | o | k | n | i |
+-----+-----+-----+--------+-----+-----+-----------+----+----+----+----+----------+-----------+-----------+
Я не думаю, что это случайно, потому что я использую одно и то же семя для обоих, и увеличение моего n не меняет этого. Кажется, что значения из df_combos не соответствуют тому, что происходит в моих строках и столбцах для результатов. Я также пытался использовать аналогичные подходы с некоторыми из семейства apply (), но использование margin = 1, 2 или c (1,2) не меняет эти проблемы, и я получаю ошибки, такие как «error in object, xi [j] замыкание »Иногда он запускается, но дает неправильные результаты. Запускается, но говорит, что« длина объекта не кратна замене », и результаты не имеют смысла. У меня также есть вложенный l oop, такой как
for (r in 1:nrow(mat))
for (c in 1:ncol(mat))
. Я попытался создать отдельную матрицу в l oop, чтобы сохранить результаты каждой строки, а затем связать в конце мастер, но это тоже не сработало. Опять же, любая помощь приветствуется.