Estimando modelos de regressão logística multinível

9

O seguinte modelo logístico multinível com uma variável explicativa no nível 1 (nível individual) e uma variável explicativa no nível 2 (nível do grupo):

π 0 j = γ 00 + γ 01 z j + u 0 j( 2 ) π 1 j = γ 10 + γ 11 z j + u 1 j( 3 )

logit(pij)=π0j+π1jxij(1)
π0j=γ00+γ01zj+u0j(2)
π1 1j=γ10+γ11zj+você1 1j(3)

onde, os resíduos no nível do grupo e são assumidos como tendo uma distribuição normal multivariada com expectativa zero. A variação dos erros residuais é especificada como e a variação dos erros residuais é especificada como .você0 0jvocê1 1jvocê0 0jσ0 02você1 1jσ1 12

Quero estimar o parâmetro do modelo e gosto de usar o Rcomando glmmPQL.

Substituindo as equações (2) e (3) na equação (1), obtém-se:

logit(pEuj)=γ00+γ10xEuj+γ01zj+γ11xEujzj+você0 0j+você1 1jxEuj(4)

Existem 30 grupos e 5 indivíduos em cada grupo.(j=1 1,...,30)

Código R:

   #Simulating data from multilevel logistic distribution 
   library(mvtnorm)
   set.seed(1234)

   J <- 30             ## number of groups
   n_j <- rep(5,J)     ## number of individuals in jth group
   N <- sum(n_j)

   g_00 <- -1
   g_01 <- 0.3
   g_10 <- 0.3
   g_11 <- 0.3

   s2_0 <- 0.13  ##variance corresponding to specific ICC
   s2_1 <- 1     ##variance standardized to 1
   s01  <- 0     ##covariance assumed zero

   z <- rnorm(J)
   x <- rnorm(N)

   #Generate (u_0j,u_1j) from a bivariate normal .
   mu <- c(0,0)
  sig <- matrix(c(s2_0,s01,s01,s2_1),ncol=2)
  u <- rmvnorm(J,mean=mu,sigma=sig,method="chol")

  pi_0 <- g_00 +g_01*z + as.vector(u[,1])
  pi_1 <- g_10 + g_11*z + as.vector(u[,2])
  eta <- rep(pi_0,n_j)+rep(pi_1,n_j)*x
  p <- exp(eta)/(1+exp(eta))

  y <- rbinom(N,1,p)

Agora, a estimativa de parâmetro.

  #### estimating parameters 
  library(MASS)
  library(nlme)

  sim_data_mat <- matrix(c(y,x,rep(z,n_j),rep(1:30,n_j)),ncol=4)
  sim_data <- data.frame(sim_data_mat)
  colnames(sim_data) <- c("Y","X","Z","cluster")
  summary(glmmPQL(Y~X*Z,random=~1|cluster,family=binomial,data=sim_data,,niter=200))

RESULTADO :

      iteration 1
      Linear mixed-effects model fit by maximum likelihood
      Data: sim_data 

      Random effects:
      Formula: ~1 | cluster
              (Intercept)  Residual
      StdDev: 0.0001541031 0.9982503

      Variance function:
      Structure: fixed weights
      Formula: ~invwt 
      Fixed effects: Y ~ X * Z 
                      Value Std.Error  DF   t-value p-value
      (Intercept) -0.8968692 0.2018882 118 -4.442404  0.0000
      X            0.5803201 0.2216070 118  2.618691  0.0100
      Z            0.2535626 0.2258860  28  1.122525  0.2712
      X:Z          0.3375088 0.2691334 118  1.254057  0.2123
      Correlation: 
           (Intr) X      Z     
      X   -0.072              
      Z    0.315  0.157       
      X:Z  0.095  0.489  0.269

      Number of Observations: 150
      Number of Groups: 30 
  • Por que leva apenas iteração enquanto mencionei para fazer iterações dentro da função pelo argumento ?2001 1200glmmPQLniter=200

  • Também o valor p da variável em nível de grupo e a interação em nível cruzado mostram que eles não são significativos. Ainda por que, neste artigo , eles mantêm a variável em nível de grupo e a interação em nível cruzado para análises adicionais?( X : Z ) ( Z ) ( X : Z )(Z)(X:Z)(Z)(X:Z)

  • Como também estão DFsendo calculados os graus de liberdade ?

  • Não corresponde ao viés relativo das várias estimativas da tabela . Tentei calcular o viés relativo como:

     #Estimated Fixed Effect parameters :
    
     hat_g_00 <- -0.8968692 #overall intercept
     hat_g_10 <- 0.5803201  # X
     hat_g_01 <-0.2535626   # Z
     hat_g_11 <-0.3375088   #X*Z
    
    fixed <-c(g_00,g_10,g_01,g_11)
    hat_fixed <-c(hat_g_00,hat_g_10,hat_g_01,hat_g_11)
    
    
    #Estimated Random Effect parameters :
    
    hat_s_0 <-0.0001541031  ##Estimated Standard deviation of random intercept 
    hat_s_1 <-  0.9982503 
    
    std  <- c(sqrt(0.13),1) 
    hat_std  <- c(0.0001541031,0.9982503) 
    
    ##Relative bias of Fixed Effect :
    rel_bias_fixed <- ((hat_fixed-fixed)/fixed)*100
    [1] -10.31308  93.44003 -15.47913  12.50293
    
    ##Relative bias of Random Effect :
    rel_bias_Random <- ((hat_std-std)/std)*100
    [1] -99.95726  -0.17497
  • Por que o viés relativo não corresponde à tabela?
abc
fonte

Respostas:

7

Talvez haja muitas perguntas aqui. Alguns comentários:

  • você pode considerar usar a glmerpartir do lme4pacote ( glmer(Y~X*Z+(1|cluster),family=binomial,data=sim_data)); ele usa a aproximação de Laplace ou a quadratura de Gauss-Hermite, que geralmente são mais precisas que o PQL (embora as respostas sejam muito semelhantes nesse caso).
  • O niterargumento especifica o número máximo de iterações; apenas uma iteração era realmente necessária
  • Não sei ao certo qual é sua pergunta sobre o termo de interação. Se você deve ou não descartar termos de interação não significativos, é uma espécie de lata de vermes e depende tanto da sua filosofia estatística quanto dos objetivos da sua análise (por exemplo, veja esta pergunta )
  • os graus de liberdade do denominador estão sendo calculados de acordo com uma simples heurística 'interior-exterior', uma regra simples 'interior-exterior' descrita na página 91 de Pinheiro e Bates (2000), disponível no Google Livros ... geralmente é uma aproximação razoável, mas o cálculo dos graus de liberdade é complexo, especialmente para GLMMs
  • se você estiver tentando replicar "Um estudo de simulação do tamanho da amostra para modelos de regressão logística multinível" por Moineddin et al. (DOI: 10.1186 / 1471-2288-7-34), é necessário executar um grande número de simulações e calcular médias, e não apenas comparar uma única execução. Além disso, você provavelmente deve tentar se aproximar de seus métodos (voltando ao meu primeiro ponto, eles afirmam que usam SAS PROC NLMIXED com quadratura adaptável de Gauss-Hermite, para que você fique melhor com, por exemplo glmer(...,nAGQ=10); ainda não corresponde exatamente, mas provavelmente estará mais perto do que glmmPQL.
Ben Bolker
fonte
I need to run a large number of simulations and compute averages300E[θ^]=θ
glmer()σ0 02σ1 12summary(glmer(Y~X*Z+(1|cluster),family=binomial,data=sim_data,nAGQ=10))
2
você está assumindo que as aproximações que usamos para a estimativa do GLMM são imparciais. Provavelmente isso não é verdade; a maioria das melhores aproximações (não PQL) é assintoticamente imparcial, mas ainda é tendenciosa para amostras de tamanho finito.
precisa
11
@ABC: Sim, os dois links contêm exemplos de como replicar um pedaço de código várias vezes. Deve ser fácil agrupar seu código em uma função e executar o comando replicate, por exemplo.
Ryan Simmons
11
@ABC: Quanto à outra parte da sua pergunta, estou um pouco confuso com o que está incomodando. Você está gerando números aleatórios; sem arredondamento ou um número infinitamente grande de repetições, você nunca obterá exatamente 0 com o viés (ou, de fato, uma estimativa exatamente precisa de QUALQUER parâmetro). No entanto, com um número suficientemente grande de replicações (por exemplo, 1000), é provável que você tenha um viés muito pequeno (próximo a 0). O artigo que você cita que está tentando replicar demonstra isso.
Ryan Simmons