1 data input

dta <- read.table("alcohol_use.txt", header = T)
head(dta)

##   sid coa sex age14  alcuse    peer    cpeer  ccoa
## 1   1   1   0     0 1.73205 1.26491  0.24691 0.549
## 2   1   1   0     1 2.00000 1.26491  0.24691 0.549
## 3   1   1   0     2 2.00000 1.26491  0.24691 0.549
## 4   2   1   1     0 0.00000 0.89443 -0.12357 0.549
## 5   2   1   1     1 0.00000 0.89443 -0.12357 0.549
## 6   2   1   1     2 1.00000 0.89443 -0.12357 0.549

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 ...

2 data management

dta$sid <- factor(dta$sid)
dta$coa <- factor(dta$coa)
dta$sex <- factor(dta$sex, 0:1, labels = c("F", "M"))
str(dta)

## 'data.frame':    246 obs. of  8 variables:
##  $ sid   : Factor w/ 82 levels "1","2","3","4",..: 1 1 1 2 2 2 3 3 3 4 ...
##  $ coa   : Factor w/ 2 levels "0","1": 2 2 2 2 2 2 2 2 2 2 ...
##  $ sex   : Factor w/ 2 levels "F","M": 1 1 1 2 2 2 2 2 2 2 ...
##  $ 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 ...

head(dta)

##   sid coa sex age14  alcuse    peer    cpeer  ccoa
## 1   1   1   F     0 1.73205 1.26491  0.24691 0.549
## 2   1   1   F     1 2.00000 1.26491  0.24691 0.549
## 3   1   1   F     2 2.00000 1.26491  0.24691 0.549
## 4   2   1   M     0 0.00000 0.89443 -0.12357 0.549
## 5   2   1   M     1 0.00000 0.89443 -0.12357 0.549
## 6   2   1   M     2 1.00000 0.89443 -0.12357 0.549

3 m0, method=MLE

library(lme4)

## Loading required package: Matrix

m0 <- lmer(alcuse ~ coa + peer + age14 + peer*age14 + ( 1 |sid), dta, REML=FALSE)
summary(m0)

## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: alcuse ~ coa + peer + age14 + peer * age14 + (1 | sid)
##    Data: dta
## 
##      AIC      BIC   logLik deviance df.resid 
##    620.5    645.0   -303.2    606.5      239 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.3210 -0.6694 -0.0541  0.5612  2.6356 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  sid      (Intercept) 0.3214   0.5669  
##  Residual             0.4765   0.6903  
## Number of obs: 246, groups:  sid, 82
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept) -0.31177    0.16970  -1.837
## coa1         0.56513    0.15585   3.626
## peer         0.69586    0.13028   5.341
## age14        0.42469    0.09289   4.572
## peer:age14  -0.15138    0.07435  -2.036
## 
## Correlation of Fixed Effects:
##            (Intr) coa1   peer   age14 
## coa1       -0.310                     
## peer       -0.726 -0.133              
## age14      -0.547  0.000  0.465       
## peer:age14  0.446  0.000 -0.571 -0.814

sjPlot::tab_model(m0, show.obs=F, show.ngroups=F, show.r2=T, show.re.var=T,show.icc=T, show.se=T)

	alcuse
Predictors	Estimates	std. Error	CI	p
(Intercept)	-0.31	0.17	-0.64 – 0.02	0.066
coa [1]	0.57	0.16	0.26 – 0.87	<0.001
peer	0.70	0.13	0.44 – 0.95	<0.001
age14	0.42	0.09	0.24 – 0.61	<0.001
peer * age14	-0.15	0.07	-0.30 – -0.01	0.042
Random Effects
σ²	0.48
τ₀₀ _sid	0.32
ICC	0.40
Marginal R² / Conditional R²	0.292 / 0.577

4 m1, method=MLE

m1 <- lmer(alcuse ~ coa + peer + age14 + peer*age14 + ( age14 |sid), dta, REML=FALSE)
summary(m1)

## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: alcuse ~ coa + peer + age14 + peer * age14 + (age14 | sid)
##    Data: dta
## 
##      AIC      BIC   logLik deviance df.resid 
##    606.7    638.3   -294.4    588.7      237 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.59554 -0.40414 -0.08352  0.45550  2.29976 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr 
##  sid      (Intercept) 0.2409   0.4908        
##           age14       0.1392   0.3730   -0.03
##  Residual             0.3373   0.5808        
## Number of obs: 246, groups:  sid, 82
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept) -0.31382    0.14612  -2.148
## coa1         0.57120    0.14623   3.906
## peer         0.69518    0.11126   6.249
## age14        0.42469    0.10559   4.022
## peer:age14  -0.15138    0.08451  -1.791
## 
## Correlation of Fixed Effects:
##            (Intr) coa1   peer   age14 
## coa1       -0.338                     
## peer       -0.709 -0.146              
## age14      -0.410  0.000  0.351       
## peer:age14  0.334  0.000 -0.431 -0.814

sjPlot::tab_model(m0, show.obs=F, show.ngroups=F, show.r2=T, show.re.var=T,show.icc=T, show.se=T)

	alcuse
Predictors	Estimates	std. Error	CI	p
(Intercept)	-0.31	0.17	-0.64 – 0.02	0.066
coa [1]	0.57	0.16	0.26 – 0.87	<0.001
peer	0.70	0.13	0.44 – 0.95	<0.001
age14	0.42	0.09	0.24 – 0.61	<0.001
peer * age14	-0.15	0.07	-0.30 – -0.01	0.042
Random Effects
σ²	0.48
τ₀₀ _sid	0.32
ICC	0.40
Marginal R² / Conditional R²	0.292 / 0.577

5 Compare model

sjPlot::tab_model(m0, m1, show.obs=F, show.ngroups=F, show.r2=T, show.re.var=T,show.icc=T, show.se=T)

	alcuse				alcuse
Predictors	Estimates	std. Error	CI	p	Estimates	std. Error	CI	p
(Intercept)	-0.31	0.17	-0.64 – 0.02	0.066	-0.31	0.15	-0.60 – -0.03	0.032
coa [1]	0.57	0.16	0.26 – 0.87	<0.001	0.57	0.15	0.28 – 0.86	<0.001
peer	0.70	0.13	0.44 – 0.95	<0.001	0.70	0.11	0.48 – 0.91	<0.001
age14	0.42	0.09	0.24 – 0.61	<0.001	0.42	0.11	0.22 – 0.63	<0.001
peer * age14	-0.15	0.07	-0.30 – -0.01	0.042	-0.15	0.08	-0.32 – 0.01	0.073
Random Effects
σ²	0.48				0.34
τ₀₀	0.32 _sid				0.24 _sid
τ₁₁					0.14 _sid.age14
ρ₀₁					-0.03 _sid
ICC	0.40				0.58
Marginal R² / Conditional R²	0.292 / 0.577				0.293 / 0.701

anova(m0,m1)

## Data: dta
## Models:
## m0: alcuse ~ coa + peer + age14 + peer * age14 + (1 | sid)
## m1: alcuse ~ coa + peer + age14 + peer * age14 + (age14 | sid)
##    npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)    
## m0    7 620.47 645.01 -303.23   606.47                         
## m1    9 606.70 638.25 -294.35   588.70 17.765  2  0.0001388 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

1012hw3

jenny hsu

1 data input

2 data management

3 m0, method=MLE

4 m1, method=MLE

5 Compare model

6 The end