Inclass Exercise3
tjlee
2020-10-19
1 load and see
## ID Age Age_m Time PBF
## 1 1 9.32 13.19 -3.87 7.94
## 2 1 10.33 13.19 -2.86 15.65
## 3 1 11.24 13.19 -1.95 13.51
## 4 1 12.19 13.19 -1.00 23.23
## 5 1 13.24 13.19 0.05 10.52
## 6 1 14.24 13.19 1.05 20.45# construct indicator variable for after menarche
dta <- dta %>%
mutate(ID = factor(ID),
Time_a = ifelse(Time > 0, Time, 0),
Menarche = as.factor(ifelse(Time > 0, "T", "F")))## 'data.frame': 1049 obs. of 7 variables:
## $ ID : Factor w/ 162 levels "1","2","3","4",..: 1 1 1 1 1 1 2 2 2 2 ...
## $ Age : num 9.32 10.33 11.24 12.19 13.24 ...
## $ Age_m : num 13.2 13.2 13.2 13.2 13.2 ...
## $ Time : num -3.87 -2.86 -1.95 -1 0.05 ...
## $ PBF : num 7.94 15.65 13.51 23.23 10.52 ...
## $ Time_a : num 0 0 0 0 0.05 ...
## $ Menarche: Factor w/ 2 levels "F","T": 1 1 1 1 2 2 1 1 1 1 ...
## - attr(*, "datalabel")= chr ""
## - attr(*, "time.stamp")= chr "23 Mar 2011 15:46"
## - attr(*, "formats")= chr "%9.0g" "%9.0g" "%9.0g" "%9.0g" ...
## - attr(*, "types")= int 254 254 254 254 254
## - attr(*, "val.labels")= chr "" "" "" "" ...
## - attr(*, "var.labels")= chr "" "" "" "" ...
## - attr(*, "version")= int 122 MC plot
#before
ggplot(subset(dta, Menarche=="F"),
aes(Time, PBF, group=ID)) +
#geom_point()+
geom_line(col="#FC766AFF")+
stat_smooth(aes(group=1),method='lm',formula=y~x,col="#FC766AFF")+
labs(x="Age (in year)",
y="% Body fat")
# baseline model if you want one
mb_0 <- lmer(PBF ~ Time + (1 | ID), data = subset(dta, Menarche=='F'))
tab_model(mb_0)| PBF | |||
|---|---|---|---|
| Predictors | Estimates | CI | p |
| (Intercept) | 20.20 | 19.06 – 21.34 | <0.001 |
| Time | -0.09 | -0.33 – 0.16 | 0.494 |
| Random Effects | |||
| σ2 | 10.02 | ||
| τ00 ID | 43.81 | ||
| ICC | 0.81 | ||
| N ID | 162 | ||
| Observations | 497 | ||
| Marginal R2 / Conditional R2 | 0.000 / 0.814 | ||
Before menarche, the percent body fat for girls averaged about 20.2 give or take 7.34 (= √43.81+10.02) with no apparent change over age in years.
Before menarche, any two body fat percent measurements taken from the same girl were quite similar and had a correlation of 0.81.
#After
ggplot(subset(dta, Menarche=="T"),
aes(Time, PBF, group=ID)) +
# geom_point()+
geom_line(col="#00DDDD")+
stat_smooth(aes(group=1),method='lm',formula=y~x,col="#00DDDD")+
labs(x="Age (in year)",
y="% Body fat")
| PBF | |||
|---|---|---|---|
| Predictors | Estimates | CI | p |
| (Intercept) | 22.59 | 21.42 – 23.75 | <0.001 |
| Time | 1.97 | 1.68 – 2.25 | <0.001 |
| Random Effects | |||
| σ2 | 7.82 | ||
| τ00 ID | 46.21 | ||
| τ11 ID.Time | 1.42 | ||
| ρ01 ID | -0.55 | ||
| ICC | 0.82 | ||
| N ID | 160 | ||
| Observations | 552 | ||
| Marginal R2 / Conditional R2 | 0.123 / 0.844 | ||
After menarche, the percent body fat for girls averaged about 22.59 give or take 7.35 (= √46.21+7.82 ) with an average annual change of 1.97 give or take 1.19 (= √1.42 ).
After menarche, any two measurements of body fat percent from the same girl had a correlation of 0.82.
betas_a <- na.omit(coef(lmList(PBF ~ Time | ID, data=subset(dta, Menarche=='T'))))
car::dataEllipse(betas_a[,1],betas_a[,2])## Registered S3 methods overwritten by 'car':
## method from
## influence.merMod lme4
## cooks.distance.influence.merMod lme4
## dfbeta.influence.merMod lme4
## dfbetas.influence.merMod lme4
3 Integration
m4 <- lmer(PBF ~ Time + Time_a + (Time + Time_a | ID), data=dta)
dta_full <- dta %>%
mutate(yhat = fitted(m4))
# individual connected lines different color for before and after
# add group regression lines
ggplot(dta_full, aes(Time, yhat, color = Menarche)) +
geom_line(aes(group = ID)) +
#geom_point(alpha = 0.5) +
#stat_smooth(method="lm", aes(group = Menarche)) +
guides(color = F) +
labs(x = "Age relative to menarche (in years)",
y = "Percent body fat")
# model with 3 random effects associated with an individual
summary(m1 <- lmer(PBF ~ Time + Time_a + (Time + Time_a | ID), data = dta))## Linear mixed model fit by REML ['lmerMod']
## Formula: PBF ~ Time + Time_a + (Time + Time_a | ID)
## Data: dta
##
## REML criterion at convergence: 6062.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.7742 -0.5900 -0.0359 0.5946 3.3798
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## ID (Intercept) 45.940 6.778
## Time 1.631 1.277 0.29
## Time_a 2.750 1.658 -0.54 -0.83
## Residual 9.473 3.078
## Number of obs: 1049, groups: ID, 162
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 21.3614 0.5645 37.838
## Time 0.4171 0.1572 2.654
## Time_a 2.0471 0.2280 8.980
##
## Correlation of Fixed Effects:
## (Intr) Time
## Time 0.351
## Time_a -0.515 -0.872## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0049345 (tol = 0.002, component 1)## Linear mixed model fit by REML ['lmerMod']
## Formula: PBF ~ Time + (Time | ID)
## Data: dta
##
## REML criterion at convergence: 6196.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.75772 -0.59624 0.02202 0.61693 3.09807
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## ID (Intercept) 38.2248 6.1826
## Time 0.6935 0.8328 -0.24
## Residual 11.6218 3.4091
## Number of obs: 1049, groups: ID, 162
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 23.04984 0.49989 46.11
## Time 1.55480 0.08483 18.33
##
## Correlation of Fixed Effects:
## (Intr)
## Time -0.200
## convergence code: 0
## Model failed to converge with max|grad| = 0.0049345 (tol = 0.002, component 1)
# residual plot
plot(m1, resid(., type = "pearson") ~ fitted(.),
abline = 0, pch = 20, cex = .8, id = 0.05,
xlab = "Ftted values", ylab = "Pearson Residuals")
The end