r-statistics-fanの日記

統計好きの現場の臨床医の覚書のようなもの

MCMCglmmを無理やり使ってみる

先日の続き

今度は混合効果モデルで推定する。

すなわち各大学にランダム効果を仮定する。

 

a <- matrix(scan("clipboard", quiet=TRUE), byrow=TRUE, ncol=3)

df <- data.frame(total = a[,1], yaba = a[,2], univ = a[,3])

df$id <- c(1:length(df$yaba))

 

fit3 <- glmmML(cbind(yaba,(total - yaba)) ~ factor(univ),cluster = id,  data = df,

               family = binomial, prior = c("logistic"))

summary(fit3)

Call:  glmmML(formula = cbind(yaba, (total - yaba)) ~ factor(univ),      family = binomial, data = df, cluster = id, prior = c("logistic")) 

 

                 coef se(coef)        z Pr(>|z|)

(Intercept)   -2.8002   0.1406 -19.9181   0.0000

factor(univ)1 -0.2332   0.4775  -0.4883   0.6250

factor(univ)2  0.6116   0.2572   2.3782   0.0174

 

Scale parameter in mixing distribution:  0.1361 logistic 

Std. Error:                              0.1061 

 

        LR p-value for H_0: sigma = 0:  0.5 

 

Residual deviance: 27.7 on 22 degrees of freedom AIC: 35.7 

 

やはり(univ)2=私立がヤバイ結果。

 

 

MCMCglmmもいじってみたが、収束が難かしい。

色々いじってそれなりの結果が出た組合せを貼る。

事前分布の設定は相当勉強しないと難しいぞ。

MCMCがふらついていい感じの動きにならんので

試行錯誤をするはめになり、意図的でないにせよ

恣意的になってしまう予感が。ぐぬぬ

 

#以下不勉強のため参考にすべきではないコード

library(MCMCglmm)

prior <- list(R=list(V=1,nu=0.02),G=list(G1=list(V=1, nu=0.02)))

fit4 <- MCMCglmm(cbind(yaba,(total - yaba)) ~ factor(univ), random =~id,

                 data = df,family =  "multinomial2", thin=10,burnin=10000,

                 nitt = 50000,prior = prior,verbose=FALSE)

summary(fit4)

plot(fit4$Sol)

posterior.mode(fit4$Sol)

predict.MCMCglmm(fit4)

Iterations = 10001:49991

 Thinning interval  = 10

 Sample size  = 4000 

 

 DIC: 761.3036 

 

 G-structure:  ~id

 

   post.mean l-95% CI u-95% CI eff.samp

id   0.05489 0.001906    0.172    721.6

 

 R-structure:  ~units

 

      post.mean l-95% CI u-95% CI eff.samp

units   0.05511 0.001891   0.1734     1105

 

 Location effects: cbind(yaba, (total - yaba)) ~ factor(univ) 

 

              post.mean l-95% CI u-95% CI eff.samp  pMCMC    

(Intercept)     -2.8198  -3.1134  -2.5264    790.3 <3e-04 ***

factor(univ)1   -0.3119  -1.3561   0.6969    760.1 0.5485    

factor(univ)2    0.6128   0.0786   1.1772   1035.4 0.0325 *  

 

いちおうこっちの解析でも私立がヤバイ結果。 

 

##下をクリップボードにコピー

62    6 0

59 5 0

68 5 0

68 5 0

51 4 0

51 4 0

55 4 0

56 4 0

60 4 0

61 4 0

67 4 0

88 4 0

57 3 0

60 2 0

61 2 0

65 2 0

70 2 0

79 2 0

41 4 1

92 2 1

53 9 2

45 8 2

41 4 2

54 4 2

54 3 2

56 3 2