Schätzen mehrstufiger logistischer Regressionsmodelle

Das folgende mehrstufige Logistikmodell mit einer erklärenden Variablen auf Ebene 1 (Einzelebene) und einer erklärenden Variablen auf Ebene 2 (Gruppenebene):

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

Dabei wird angenommen, dass die Residuen und eine multivariate Normalverteilung mit der Erwartung Null haben. Die Varianz der Restfehler wird als , und die Varianz der Restfehler u_ {1j} wird als \ sigma ^ 2_1 angegeben . u 1 j u 0 j σ 2 0 u 1 j σ 2 $u_{0j}$$u_{1j}$$u_{0j}$$\sigma^2_0$$u_{1j}$$\sigma^2_1$

Ich möchte den Parameter des Modells schätzen und verwende gerne den RBefehl glmmPQL.

Einsetzen von Gleichung (2) und (3) in Gleichung (1) ergibt,

Es gibt 30 Gruppen $(j=1,...,30)$ und 5 Einzelpersonen in jeder Gruppe.

R-Code:

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

Nun die Parameterschätzung.

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

AUSGABE :

      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 

• Warum dauert es nur Iteration, während ich erwähnte, dass das Argument Iterationen innerhalb der Funktion durchführt ?$1$$200$glmmPQLniter=200

• Auch der p-Wert der Variablen auf Gruppenebene und der Interaktion zwischen Ebenen zeigt, dass sie nicht signifikant sind. Warum behalten sie in diesem Artikel die Variable auf Gruppenebene und die Interaktion zwischen Ebenen zur weiteren Analyse bei?( X : Z ) ( Z ) ( X : Z $(Z)$$(X:Z)$$(Z)$$(X:Z)$

• Wie werden die Freiheitsgrade DFberechnet?

• Es stimmt nicht mit der relativen Verzerrung der verschiedenen Schätzungen der Tabelle überein . Ich habe versucht, die relative Verzerrung wie folgt zu berechnen:

   #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
  

• Warum stimmt die relative Verzerrung nicht mit der Tabelle überein?

Hier gibt es vielleicht zu viele Fragen. Einige Kommentare:

• Sie könnten in Betracht ziehen, glmeraus dem lme4Paket ( glmer(Y~X*Z+(1|cluster),family=binomial,data=sim_data)) zu verwenden; Es verwendet die Laplace-Näherung oder die Gauß-Hermite-Quadratur, die im Allgemeinen genauer als PQL sind (obwohl die Antworten in diesem Fall sehr ähnlich sind).
• Das niterArgument gibt die maximale Anzahl von Iterationen an. es war tatsächlich nur eine Iteration notwendig
• Ich bin mir nicht sicher, was Ihre Frage zum Interaktionsbegriff ist. Ob Sie nicht signifikante Interaktionsterme löschen sollten oder nicht, ist eine Dose Würmer und hängt sowohl von Ihrer statistischen Philosophie als auch von den Zielen Ihrer Analyse ab (siehe z. B. diese Frage ).
• Die Freiheitsgrade des Nenners werden nach einer einfachen "inner-äußeren" Heuristik berechnet, einer einfachen "inner-äußeren" Regel, die auf Seite 91 von Pinheiro und Bates (2000) beschrieben ist und in Google Books verfügbar ist Eine vernünftige Annäherung, aber die Berechnung von Freiheitsgraden ist komplex, insbesondere für GLMMs
• Wenn Sie versuchen, "Eine Simulationsstudie zur Stichprobengröße für mehrstufige logistische Regressionsmodelle" von Moineddin et al. (DOI: 10.1186 / 1471-2288-7-34) müssen Sie eine große Anzahl von Simulationen ausführen und Durchschnittswerte berechnen, nicht nur einen einzelnen Lauf vergleichen. Darüber hinaus sollten Sie wahrscheinlich versuchen, ihre Methoden näher zu kommen (zurück zu meinem ersten Punkt kommen, sie erklären , dass sie verwenden SAS PROC nlmixed mit adaptiver Gauß-Hermite Quadratur, so dass Sie mit zB besser dran wären glmer(...,nAGQ=10), wird es immer noch nicht genau übereinstimmen, aber es wird wahrscheinlich näher sein als glmmPQL.