Estimation des modèles de régression logistique à plusieurs niveaux

Le modèle logistique à plusieurs niveaux suivant avec une variable explicative au niveau 1 (niveau individuel) et une variable explicative au niveau 2 (niveau groupe):

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

où, les résidus au niveau du groupe et u 1 j sont supposés avoir une distribution normale multivariée avec une espérance nulle. La variance des erreurs résiduelles est spécifiée comme , et la variance des erreurs résiduelles est spécifiée comme σ 2 1 .$u_{0j}$$u_{1j}$ σ 2 0 u 1 $u_{0j}$$\sigma^2_0$$u_{1j}$$\sigma^2_1$

Je veux estimer le paramètre du modèle et j'aime utiliser la Rcommande glmmPQL.

La substitution de l'équation (2) et (3) dans l'équation (1) donne,

$(j=1,...,30)$

Code 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)

Maintenant, l'estimation des paramètres.

  #### 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))

PRODUCTION :

      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 

• $1$$200$glmmPQLniter=200

• $(Z)$$(X:Z)$$(Z)$$(X:Z)$

• Aussi Comment les degrés de liberté DFsont-ils calculés?

• Il ne correspond pas au biais relatif des différentes estimations du tableau . J'ai essayé de calculer le biais relatif comme suit:

   #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

• Pourquoi le biais relatif ne correspond-il pas au tableau?

Il y a peut-être trop de questions ici. Certains commentaires:

• vous pourriez envisager d'utiliser à glmerpartir du lme4package ( glmer(Y~X*Z+(1|cluster),family=binomial,data=sim_data)); il utilise l'approximation de Laplace ou la quadrature de Gauss-Hermite, qui sont généralement plus précises que PQL (bien que les réponses soient très similaires dans ce cas).
• L' niterargument spécifie le nombre maximal d'itérations; une seule itération était réellement nécessaire
• Je ne sais pas quelle est votre question sur le terme d'interaction. Que vous deviez ou non supprimer des termes d'interaction non significatifs est un peu une boîte de vers et dépend à la fois de votre philosophie statistique et des objectifs de votre analyse (par exemple, voir cette question )
• les degrés de liberté du dénominateur sont calculés selon une simple heuristique `` intérieur-extérieur '', une simple règle `` intérieur-extérieur '' décrite à la page 91 de Pinheiro et Bates (2000), qui est disponible sur Google Books ... c'est généralement une approximation raisonnable mais le calcul des degrés de liberté est complexe, notamment pour les GLMM
• si vous essayez de reproduire "Une étude de simulation de la taille de l'échantillon pour les modèles de régression logistique à plusieurs niveaux" par Moineddin et al. (DOI: 10.1186 / 1471-2288-7-34), vous devez exécuter un grand nombre de simulations et calculer des moyennes, pas seulement comparer une seule exécution. De plus, vous devriez probablement essayer de vous rapprocher de leurs méthodes (pour en revenir à mon premier point, ils déclarent qu'ils utilisent SAS PROC NLMIXED avec une quadrature adaptative de Gauss-Hermite, donc vous feriez mieux, par exemple glmer(...,nAGQ=10); il ne le sera toujours pas correspondent exactement, mais ce sera probablement plus proche que glmmPQL.