Имитация повторных измерений двоичных данных с определенным соотношением шансов - PullRequest
Новая
       102

Имитация повторных измерений двоичных данных с определенным соотношением шансов

0 голосов
/ 15 сентября 2018

Я пытаюсь смоделировать бинарный результат, в котором у меня есть N субъектов (с вероятностями для каждого субъекта), измеренных в два разных периода (скажем, до и после). Я хочу увеличить предметные вероятности на определенное значение отношения шансов (ИЛИ) между двумя периодами.

После моделирования я использовал glm и lme4::glmer, чтобы проверить, правильно ли была оценена моя заранее заданная пропорция шансов. Я ожидал, что только ИЛИ, оцененный в glm, будет предвзятым. Тем не менее, OR, оцененное в lme4::glmer, также было смещено, так как мои предопределенные значения OR увеличиваются. Как я могу исправить это смещение?

Большое спасибо,

Ниже мое моделирование 

rm(list=ls(all=TRUE))
library(lme4)
library(ggplot2)

N = 2000                                                #Number of subjects
X = 1:20                                                #Odds ratio values tested
set.seed(20)    
P = runif(N,-4,4)                                       #Subject-specific probability (in logit scale)

#Vectors that will be used to create a data frame
ind = rep(paste0("Sub",1:N),2)                          #Vector of individuals
x1 = c(rep(0,N),rep(1,N))                               #Categorical Predictor Variable x1
OR.glm = NULL;OR.glmer = NULL

#Loop over X
for (OR in X){
  value = rbinom(N,1,plogis(P))                         #Simulating values for x1=0
  value.simu = rbinom(N,1,plogis(P+log(OR)))            #Simulating values for x1=1
  df = data.frame(ind=ind,y=c(value,value.simu),x1=x1)  #Creating data frame

  #Using glm
  GLM = glm(y~factor(x1),data=df,family="binomial")
  OR.glm = c(OR.glm,exp(GLM$coef[2]))

  #Using glmer for each subject
  GLMER = glmer(y~factor(x1)+(1|ind),data=df,family="binomial")
  OR.glmer = c(OR.glmer,exp(summary(GLMER)$coef[2,1]))
}

DF = data.frame(method = rep(c("glm","glmer"),each=length(X)),
        data = c(OR.glm,OR.glmer),x = rep(X,2))
ggplot(DF,aes(x = x,y = data,group=method, colour=method))+ theme_bw()+
    geom_point() + stat_smooth(method = 'loess') +
    geom_abline(slope=1, intercept=0) + ylim(0, max(X)) + xlim(0, max(X)) +
    xlab("Expected OR") + ylab("Observed OR")

1 Ответ

0 голосов
/ 15 сентября 2018

Насколько я вижу, вы не моделируете нормальные случайные эффекты, что является предположением для модели логистической регрессии смешанных эффектов, установленной glmer().

Приведенный ниже код имитирует данные с обычным случайным эффектом и соответствует модели с glmer() из lme4 и mixed_model() из GLMMadaptive , который по умолчанию использует адаптивную квадратуру Гаусса в оценке (специально код работает с матрицами проектирования для фиксированных и случайных эффектов, чтобы облегчить его расширение, если вы хотите): 

set.seed(1234)
n <- 100 # number of subjects
K <- 8 # number of measurements per subject

# we constuct a data frame with the design: 
DF <- data.frame(id = rep(seq_len(n), each = K),
                 sex = rep(gl(2, n/2, labels = c("male", "female")), each = K))

# design matrices for the fixed and random effects
X <- model.matrix(~ sex, data = DF)
Z <- model.matrix(~ 1, data = DF)

betas <- c(-2.13, 1) # fixed effects coefficients
D11 <- 0.48 # variance of random intercepts

# we simulate random effects
b <- rnorm(n, sd = sqrt(D11))
# linear predictor
eta_y <- drop(X %*% betas + rowSums(Z * b[DF$id]))
# we simulate binary longitudinal data
DF$y <- rbinom(n * K, 1, plogis(eta_y))

###############################################################################

library("lme4")

fm <- glmer(y ~ sex + (1 | id), data = DF, family = binomial())
summary(fm)

library("GLMMadaptive")

gm <- mixed_model(y ~ sex, random = ~ 1 | id, data = DF, family = binomial())
summary(gm)
...