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
