Q4
dta <-read.table("alcohol_use.txt",h=T)
str(dta)
## 'data.frame': 246 obs. of 8 variables:
## $ sid : int 1 1 1 2 2 2 3 3 3 4 ...
## $ coa : int 1 1 1 1 1 1 1 1 1 1 ...
## $ sex : int 0 0 0 1 1 1 1 1 1 1 ...
## $ age14 : int 0 1 2 0 1 2 0 1 2 0 ...
## $ alcuse: num 1.73 2 2 0 0 ...
## $ peer : num 1.265 1.265 1.265 0.894 0.894 ...
## $ cpeer : num 0.247 0.247 0.247 -0.124 -0.124 ...
## $ ccoa : num 0.549 0.549 0.549 0.549 0.549 0.549 0.549 0.549 0.549 0.549 ...
require(lme4)
## Loading required package: lme4
## Loading required package: Matrix
summary(m1<- lmer(data=dta, alcuse ~ coa + peer*age14 +(1|sid)))
## Linear mixed model fit by REML ['lmerMod']
## Formula: alcuse ~ coa + peer * age14 + (1 | sid)
## Data: dta
##
## REML criterion at convergence: 621.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.30588 -0.66394 -0.05233 0.55906 2.60999
##
## Random effects:
## Groups Name Variance Std.Dev.
## sid (Intercept) 0.3377 0.5811
## Residual 0.4823 0.6945
## Number of obs: 246, groups: sid, 82
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) -0.31177 0.17225 -1.810
## coa 0.56513 0.15878 3.559
## peer 0.69586 0.13219 5.264
## age14 0.42469 0.09346 4.544
## peer:age14 -0.15138 0.07481 -2.024
##
## Correlation of Fixed Effects:
## (Intr) coa peer age14
## coa -0.312
## peer -0.725 -0.134
## age14 -0.543 0.000 0.461
## peer:age14 0.442 0.000 -0.566 -0.814
summary(m2 <- lmer(data=dta,alcuse ~ coa + peer*age14 + ( 1 + age14| sid )))
## Linear mixed model fit by REML ['lmerMod']
## Formula: alcuse ~ coa + peer * age14 + (1 + age14 | sid)
## Data: dta
##
## REML criterion at convergence: 603.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.59995 -0.40432 -0.07739 0.44372 2.27436
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## sid (Intercept) 0.2595 0.5094
## age14 0.1469 0.3832 -0.05
## Residual 0.3373 0.5808
## Number of obs: 246, groups: sid, 82
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) -0.31382 0.14871 -2.110
## coa 0.57120 0.14898 3.834
## peer 0.69518 0.11322 6.140
## age14 0.42469 0.10690 3.973
## peer:age14 -0.15138 0.08556 -1.769
##
## Correlation of Fixed Effects:
## (Intr) coa peer age14
## coa -0.339
## peer -0.708 -0.146
## age14 -0.408 0.000 0.350
## peer:age14 0.332 0.000 -0.429 -0.814
AIC(m1, m2)
## df AIC
## m1 7 635.6048
## m2 9 621.6977
Q5
dta<-as.data.frame(sleepstudy)
str(dta)
## 'data.frame': 180 obs. of 3 variables:
## $ Reaction: num 250 259 251 321 357 ...
## $ Days : num 0 1 2 3 4 5 6 7 8 9 ...
## $ Subject : Factor w/ 18 levels "308","309","310",..: 1 1 1 1 1 1 1 1 1 1 ...
#plot look over the data of sub
require(ggplot2)
## Loading required package: ggplot2
ggplot(dta, aes(Days, Reaction, color = Subject))+
geom_point() +
stat_smooth(aes(group = Subject), method = "lm", se = F) +
labs(x = "Days", y = "Reaction time") +
theme_bw()
ggplot(dta, aes(Days, Reaction)) +
geom_point(alpha = 0.4,size = 1) +
stat_smooth(aes(group = 1), method = "lm", se = F, col = "skyblue") +
facet_wrap( ~ Subject)+theme_bw()
summary(m0 <- lmer(Reaction ~ Days+ (Days | Subject), data = dta))
## Linear mixed model fit by REML ['lmerMod']
## Formula: Reaction ~ Days + (Days | Subject)
## Data: dta
##
## REML criterion at convergence: 1743.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.9536 -0.4634 0.0231 0.4634 5.1793
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## Subject (Intercept) 612.09 24.740
## Days 35.07 5.922 0.07
## Residual 654.94 25.592
## Number of obs: 180, groups: Subject, 18
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 251.405 6.825 36.84
## Days 10.467 1.546 6.77
##
## Correlation of Fixed Effects:
## (Intr)
## Days -0.138