Memperkirakan Model Regresi Logistik Multilevel

Model logistik bertingkat berikut dengan satu variabel penjelas pada level 1 (level individu) dan satu variabel penjelas pada level 2 (level grup):

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

di mana, residu level grup dan diasumsikan memiliki distribusi normal multivariat dengan ekspektasi nol. kesalahan sisa ditentukan sebagai , dan varian dari kesalahan sisa ditentukan sebagai . 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$

Saya ingin memperkirakan parameter model dan saya suka menggunakan Rperintah glmmPQL.

Mengganti persamaan (2) dan (3) dalam persamaan (1) menghasilkan,

Ada 30 kelompok dan 5 individu di setiap kelompok.$(j=1,...,30)$

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

Sekarang estimasi parameter.

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

OUTPUT:

      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 

• Mengapa hanya butuh iterasi sementara saya menyebutkan untuk mengambil iterasi di dalam fungsi oleh argumen ?$1$$200$glmmPQLniter=200

• Juga nilai-p dari variabel tingkat-kelompok dan interaksi lintas-tingkat menunjukkan mereka tidak signifikan. Masih mengapa dalam artikel ini , mereka menjaga variabel level grup dan interaksi lintas level untuk analisis lebih lanjut?( X : Z ) ( Z ) ( X : Z $(Z)$$(X:Z)$$(Z)$$(X:Z)$

• Juga Bagaimana derajat kebebasan DFdihitung?

• Itu tidak cocok dengan bias relatif dari berbagai estimasi tabel . Saya mencoba menghitung bias relatif sebagai:

   #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
  

• Mengapa bias relatif tidak cocok dengan tabel?

Mungkin ada terlalu banyak pertanyaan di sini. Beberapa komentar:

• Anda dapat mempertimbangkan menggunakan glmerdari lme4paket ( glmer(Y~X*Z+(1|cluster),family=binomial,data=sim_data)); menggunakan pendekatan Laplace atau quadrature Gauss-Hermite, yang umumnya lebih akurat daripada PQL (meskipun jawabannya sangat mirip dalam kasus ini).
• The niterArgumen menentukan maksimum jumlah iterasi; hanya satu iterasi yang benar-benar diperlukan
• Saya tidak yakin apa pertanyaan Anda tentang istilah interaksi. Apakah Anda harus meninggalkan istilah interaksi yang tidak signifikan atau tidak adalah sedikit dari cacing, dan tergantung pada filosofi statistik Anda dan pada tujuan analisis Anda (mis. Lihat pertanyaan ini )
• derajat kebebasan penyebut dihitung berdasarkan heuristik sederhana 'luar-dalam', aturan 'luar-dalam' sederhana yang dijelaskan pada halaman 91 dari Pinheiro dan Bates (2000), yang tersedia di Google Books ... umumnya perkiraan yang masuk akal tetapi perhitungan derajat kebebasan itu kompleks, terutama untuk GLMM
• jika Anda mencoba mereplikasi "Sebuah studi simulasi ukuran sampel untuk model regresi logistik bertingkat" oleh Moineddin et al. (DOI: 10.1186 / 1471-2288-7-34), Anda perlu menjalankan sejumlah besar simulasi dan menghitung rata-rata, tidak hanya membandingkan satu jalankan. Selain itu, Anda mungkin harus mencoba untuk lebih dekat dengan metode mereka (kembali ke poin pertama saya, mereka menyatakan bahwa mereka menggunakan SAS PROC NLMIXED dengan quadrature Gauss-Hermite adaptif, sehingga Anda akan lebih baik dengan misalnya glmer(...,nAGQ=10); itu masih tidak akan lebih baik dengan misalnya ; cocok persis, tetapi mungkin akan lebih dekat daripada glmmPQL.