GLMM
Ching Wen, Su
2016/12/4
require(pacman)## Loading required package: pacmanpacman::p_load(tidyverse, lme4, ltm)Q1
#data
dta<- read.table("http://titan.ccunix.ccu.edu.tw/~psycfs/lmm/Data/starter.txt", h=T)
head(dta, 10) sid item resp
1 1 1 1
2 1 2 1
3 1 3 1
4 1 4 1
5 1 5 1
6 1 6 1
7 1 7 1
8 1 8 1
9 1 9 0
10 2 1 1dta$sid<-factor(dta$sid)
dta$item<- factor(dta$item)
str(dta) 'data.frame': 1350 obs. of 3 variables:
$ sid : Factor w/ 150 levels "1","2","3","4",..: 1 1 1 1 1 1 1 1 1 2 ...
$ item: Factor w/ 9 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 9 1 ...
$ resp: int 1 1 1 1 1 1 1 1 0 1 ...#plot
ggplot(dta, aes(item, resp))+
stat_summary(fun.data = mean_se)+
theme_bw()+
labs(x = "Item", y = "Accuracy")
#model
summary(m0 <- glmer(resp ~ -1 + item + (1 | sid), data = dta, family = binomial)) Generalized linear mixed model fit by maximum likelihood (Laplace
Approximation) [glmerMod]
Family: binomial ( logit )
Formula: resp ~ -1 + item + (1 | sid)
Data: dta
AIC BIC logLik deviance df.resid
1430.0 1482.1 -705.0 1410.0 1340
Scaled residuals:
Min 1Q Median 3Q Max
-6.0216 -0.5795 0.3395 0.5451 3.6577
Random effects:
Groups Name Variance Std.Dev.
sid (Intercept) 1.165 1.079
Number of obs: 1350, groups: sid, 150
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
item1 2.9031 0.3405 8.526 < 2e-16 ***
item2 1.7032 0.2383 7.147 8.85e-13 ***
item3 1.3373 0.2214 6.039 1.55e-09 ***
item4 1.2102 0.2167 5.584 2.36e-08 ***
item5 1.0017 0.2102 4.765 1.89e-06 ***
item6 0.9241 0.2082 4.439 9.03e-06 ***
item7 1.8649 0.2475 7.534 4.91e-14 ***
item8 0.2030 0.1973 1.029 0.303
item9 -0.8778 0.2077 -4.225 2.39e-05 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
item1 item2 item3 item4 item5 item6 item7 item8
item2 0.108
item3 0.113 0.158
item4 0.115 0.160 0.171
item5 0.116 0.164 0.175 0.178
item6 0.117 0.165 0.176 0.179 0.184
item7 0.105 0.144 0.153 0.155 0.158 0.159
item8 0.116 0.167 0.180 0.183 0.189 0.191 0.160
item9 0.101 0.148 0.161 0.166 0.172 0.174 0.142 0.189
convergence code: 0
unable to evaluate scaled gradient
Model failed to converge: degenerate Hessian with 2 negative eigenvaluesQ2
#data
dta<-read.table("http://titan.ccunix.ccu.edu.tw/~psycfs/lmm/Data/adolescentSmoking.txt" ,h=T)
head(dta) newid sex parsmk wave smkreg
1 1 1 0 1 0
2 1 1 0 2 0
3 1 1 0 4 0
4 1 1 0 5 0
5 1 1 0 6 0
6 2 0 0 1 0dta$sex<-factor(dta$sex,levels = c(0,1), labels = c("boys", "girls"))
dta$parsmk<-factor(dta$parsmk)
dta$wave<-factor(dta$wave)
str(dta) 'data.frame': 8730 obs. of 5 variables:
$ newid : int 1 1 1 1 1 2 2 2 2 2 ...
$ sex : Factor w/ 2 levels "boys","girls": 2 2 2 2 2 1 1 1 1 1 ...
$ parsmk: Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 1 1 ...
$ wave : Factor w/ 6 levels "1","2","3","4",..: 1 2 4 5 6 1 2 3 4 5 ...
$ smkreg: int 0 0 0 0 0 0 0 0 0 0 ...dtam<- aggregate(smkreg ~ wave + sex, FUN = mean, data = dta)
#plot
ggplot(dtam, aes(wave, smkreg, group = sex, linetype = sex))+
geom_point()+
geom_line()+
theme_bw()+
ylim(0,.2)+
labs(x = "Wave of study", y = "Propotion of smokers")
#model
summary(m0 <- glmer(smkreg~sex+parsmk+wave+(1|newid), family = "binomial", data = dta)) Generalized linear mixed model fit by maximum likelihood (Laplace
Approximation) [glmerMod]
Family: binomial ( logit )
Formula: smkreg ~ sex + parsmk + wave + (1 | newid)
Data: dta
AIC BIC logLik deviance df.resid
3803.4 3867.1 -1892.7 3785.4 8721
Scaled residuals:
Min 1Q Median 3Q Max
-3.4829 -0.0328 -0.0209 -0.0144 4.2244
Random effects:
Groups Name Variance Std.Dev.
newid (Intercept) 40.98 6.402
Number of obs: 8730, groups: newid, 1760
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -8.0360 0.3551 -22.630 < 2e-16 ***
sexgirls 0.2196 0.2694 0.815 0.414944
parsmk1 1.1726 0.2762 4.246 2.18e-05 ***
wave2 -0.8710 0.2454 -3.549 0.000386 ***
wave3 -0.3562 0.2394 -1.488 0.136775
wave4 0.1828 0.2352 0.777 0.437049
wave5 0.6205 0.2339 2.653 0.007983 **
wave6 1.1000 0.2330 4.720 2.36e-06 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) sxgrls prsmk1 wave2 wave3 wave4 wave5
sexgirls -0.415
parsmk1 -0.261 0.033
wave2 -0.381 -0.009 0.012
wave3 -0.422 -0.007 0.009 0.642
wave4 -0.465 -0.007 0.008 0.649 0.667
wave5 -0.504 -0.007 0.007 0.647 0.670 0.693
wave6 -0.548 -0.012 0.005 0.639 0.666 0.693 0.711
convergence code: 0
Model failed to converge with max|grad| = 4.75774 (tol = 0.001, component 1)