MLM Chicks
Olesya Volchenko
May 1, 2019
library(sjPlot)
library(lme4)
library(ICC)
library(ggplot2)
library(ggthemes)ChickWeight <- ChickWeight
table(ChickWeight$Chick)##
## 18 16 15 13 9 20 10 8 17 19 4 6 11 3 1 12 2 5 14 7 24 30 22 23 27
## 2 7 8 12 12 12 12 11 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12
## 28 26 25 29 21 33 37 36 31 39 38 32 40 34 35 44 45 43 41 47 49 46 50 42 48
## 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 10 12 12 12 12 12 12 12 12 12table(ChickWeight$Time)##
## 0 2 4 6 8 10 12 14 16 18 20 21
## 50 50 49 49 49 49 49 48 47 47 46 45hist(ChickWeight$Time)
hist(ChickWeight$weight)
barplot(table(ChickWeight$Diet))
ggplot(data = ChickWeight, aes(x = Time, y = weight, group = Chick)) + geom_line() 
fit1 <- lmer(weight ~ (1 | Chick), data = ChickWeight)
summary(fit1)## Linear mixed model fit by REML ['lmerMod']
## Formula: weight ~ (1 | Chick)
## Data: ChickWeight
##
## REML criterion at convergence: 6544
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.8191 -0.7669 -0.1873 0.6167 3.1114
##
## Random effects:
## Groups Name Variance Std.Dev.
## Chick (Intercept) 541.2 23.26
## Residual 4534.2 67.34
## Number of obs: 578, groups: Chick, 50
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 120.988 4.339 27.89ICCbare(Chick, weight, ChickWeight)## [1] 0.1077609ranef(fit1)## $Chick
## (Intercept)
## 18 -16.1849326
## 16 -32.4423320
## 15 -29.7281441
## 13 -31.3003428
## 9 -23.4489431
## 20 -25.0682943
## 10 -22.3203044
## 8 -16.4546580
## 17 -16.7752533
## 19 -20.1611694
## 4 -12.7514109
## 6 -4.2620849
## 11 5.2577373
## 3 -3.0353037
## 1 -5.4888661
## 12 -4.0657999
## 2 -0.6308125
## 5 3.3439586
## 14 17.8690481
## 7 17.0839082
## 24 -32.2326966
## 30 -10.2978485
## 22 -9.8562073
## 23 -5.6360799
## 27 -6.2249349
## 28 5.2577373
## 26 5.8956635
## 25 13.0109945
## 29 12.2749258
## 21 37.3994050
## 33 -6.6175049
## 37 -10.8867035
## 36 8.2020122
## 31 4.4725973
## 39 7.8094422
## 38 12.5693533
## 32 21.5493918
## 40 21.5493918
## 34 28.1740103
## 35 42.5028149
## 44 -10.2771645
## 45 -0.8270975
## 43 12.9619233
## 41 4.3744548
## 47 4.0800273
## 49 9.8704346
## 46 7.7112997
## 50 15.6117707
## 42 16.5441244
## 48 21.5984630plot_model(fit1, type = "re")
confint(fit1)## 2.5 % 97.5 %
## .sig01 15.07037 31.89203
## .sigma 63.46621 71.62446
## (Intercept) 112.34043 129.53104FSRI model
adding Diet as a fixed predictor
fit2 <- lmer(weight ~ Diet + (1 | Chick), data = ChickWeight)
summary(fit2)## Linear mixed model fit by REML ['lmerMod']
## Formula: weight ~ Diet + (1 | Chick)
## Data: ChickWeight
##
## REML criterion at convergence: 6506.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.8495 -0.7499 -0.1835 0.6363 3.0852
##
## Random effects:
## Groups Name Variance Std.Dev.
## Chick (Intercept) 305.8 17.49
## Residual 4526.6 67.28
## Number of obs: 578, groups: Chick, 50
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 101.63 6.04 16.826
## Diet2 20.98 10.24 2.050
## Diet3 41.32 10.24 4.036
## Diet4 33.40 10.27 3.253
##
## Correlation of Fixed Effects:
## (Intr) Diet2 Diet3
## Diet2 -0.590
## Diet3 -0.590 0.348
## Diet4 -0.588 0.347 0.347anova(fit1, fit2)## Data: ChickWeight
## Models:
## fit1: weight ~ (1 | Chick)
## fit2: weight ~ Diet + (1 | Chick)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## fit1 3 6554.8 6567.9 -3274.4 6548.8
## fit2 6 6542.3 6568.4 -3265.1 6530.3 18.523 3 0.0003431 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1adding Time as a fixed predictor
fit3 <- lmer(weight ~ Time + (1 | Chick), data = ChickWeight)
summary(fit3)## Linear mixed model fit by REML ['lmerMod']
## Formula: weight ~ Time + (1 | Chick)
## Data: ChickWeight
##
## REML criterion at convergence: 5619.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.0086 -0.5528 -0.0888 0.4975 3.4588
##
## Random effects:
## Groups Name Variance Std.Dev.
## Chick (Intercept) 717.9 26.79
## Residual 799.4 28.27
## Number of obs: 578, groups: Chick, 50
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 27.8451 4.3877 6.346
## Time 8.7261 0.1755 49.716
##
## Correlation of Fixed Effects:
## (Intr)
## Time -0.422anova(fit1, fit3)## Data: ChickWeight
## Models:
## fit1: weight ~ (1 | Chick)
## fit3: weight ~ Time + (1 | Chick)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## fit1 3 6554.8 6567.9 -3274.4 6548.8
## fit3 4 5630.3 5647.8 -2811.2 5622.3 926.47 1 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Note: Since fit2 and fit3 are non-nested, you cannot compare those with anova
let’s try to randomize Time and visualise the results (and exclude Diet, for now)
fit31 <- lmer(weight ~ Time + (1 + Time | Chick), data = ChickWeight)
summary(fit31)## Linear mixed model fit by REML ['lmerMod']
## Formula: weight ~ Time + (1 + Time | Chick)
## Data: ChickWeight
##
## REML criterion at convergence: 4827.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.6731 -0.5563 -0.0277 0.5010 3.4939
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## Chick (Intercept) 140.53 11.855
## Time 14.14 3.761 -0.95
## Residual 163.51 12.787
## Number of obs: 578, groups: Chick, 50
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 29.1780 1.9573 14.91
## Time 8.4531 0.5408 15.63
##
## Correlation of Fixed Effects:
## (Intr)
## Time -0.871anova(fit3, fit31)## Data: ChickWeight
## Models:
## fit3: weight ~ Time + (1 | Chick)
## fit31: weight ~ Time + (1 + Time | Chick)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## fit3 4 5630.3 5647.8 -2811.2 5622.3
## fit31 6 4841.8 4868.0 -2414.9 4829.8 792.5 2 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1plot_model(fit31, type = "re")
Let’s make a plot combining both random slopes and random intercepts
ranef(fit31)## $Chick
## (Intercept) Time
## 18 2.4034442 -0.91840755
## 16 19.7772171 -7.52516438
## 15 18.1847484 -6.40855717
## 13 17.4510920 -6.36128924
## 9 17.6676807 -5.43632956
## 20 12.9229658 -4.94955040
## 10 12.2466174 -4.51738492
## 8 10.9758864 -3.36619512
## 17 11.7239346 -3.76673155
## 19 8.6302088 -3.73309794
## 4 6.3910342 -2.50790394
## 6 7.3100869 -1.60652069
## 11 5.0563017 -0.12224630
## 3 -0.9565009 -0.26541669
## 1 0.4672194 -0.76946547
## 12 -1.0208768 -0.38163536
## 2 -1.3077270 0.07804658
## 5 -5.7524847 1.20286456
## 14 -9.9974553 3.58375216
## 7 -15.1588362 4.23219176
## 24 21.4138185 -7.04609739
## 30 7.7694579 -2.40780856
## 22 7.9613838 -2.38169106
## 23 5.7199436 -1.54461971
## 27 2.7577612 -1.18930911
## 28 -3.8515666 1.16190027
## 26 -5.3170772 1.45085621
## 25 -8.4511049 2.76919421
## 29 -13.2912821 3.37730565
## 21 -19.3675882 7.31304297
## 33 8.8526561 -2.11575028
## 37 4.4606946 -2.00253530
## 36 -4.2902173 1.58370413
## 31 -5.3554121 1.28307439
## 39 -7.4202754 1.98712006
## 38 -11.6221523 3.17254295
## 32 -13.9581164 4.60292655
## 40 -15.0199815 4.75600363
## 34 -19.7142557 6.23950329
## 35 -25.3546520 8.79764905
## 44 7.0100037 -1.47731752
## 45 2.8421429 -0.54409777
## 43 3.6751100 1.01511887
## 41 2.2021387 0.18163551
## 47 1.2556671 0.28222075
## 49 -2.7712230 1.56791678
## 46 -3.3875267 1.39381203
## 50 -7.9157157 3.00874876
## 42 -9.7325898 3.38421365
## 48 -16.1145979 4.91977824fit31coef <- coef(fit31)
fit31coef <- fit31coef$`Chick`
names(fit31coef) <- c("Intercept", "Time")
range(ChickWeight$Time)## [1] 0 21range(ChickWeight$weight)## [1] 35 373plot(range(ChickWeight$Time), range(ChickWeight$weight), type = "n", xlab = "Time", ylab = "Weight")
for (i in 1:50){
abline(a = fit31coef$Intercept[i], b = fit31coef$Time[i])
}
fit4 <- lmer(weight ~ Time + Diet + (1 | Chick), data = ChickWeight)
summary(fit4)## Linear mixed model fit by REML ['lmerMod']
## Formula: weight ~ Time + Diet + (1 | Chick)
## Data: ChickWeight
##
## REML criterion at convergence: 5584
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.0591 -0.5779 -0.1182 0.4962 3.4515
##
## Random effects:
## Groups Name Variance Std.Dev.
## Chick (Intercept) 525.4 22.92
## Residual 799.4 28.27
## Number of obs: 578, groups: Chick, 50
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 11.2438 5.7887 1.942
## Time 8.7172 0.1755 49.684
## Diet2 16.2100 9.4643 1.713
## Diet3 36.5433 9.4643 3.861
## Diet4 30.0129 9.4708 3.169
##
## Correlation of Fixed Effects:
## (Intr) Time Diet2 Diet3
## Time -0.307
## Diet2 -0.550 -0.015
## Diet3 -0.550 -0.015 0.339
## Diet4 -0.550 -0.011 0.339 0.339anova(fit3, fit4)## Data: ChickWeight
## Models:
## fit3: weight ~ Time + (1 | Chick)
## fit4: weight ~ Time + Diet + (1 | Chick)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## fit3 4 5630.3 5647.8 -2811.2 5622.3
## fit4 7 5619.2 5649.7 -2802.6 5605.2 17.143 3 0.0006603 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(fit2, fit4)## Data: ChickWeight
## Models:
## fit2: weight ~ Diet + (1 | Chick)
## fit4: weight ~ Time + Diet + (1 | Chick)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## fit2 6 6542.3 6568.4 -3265.1 6530.3
## fit4 7 5619.2 5649.7 -2802.6 5605.2 925.09 1 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1ranef(fit4)## $Chick
## (Intercept)
## 18 9.6771423
## 16 -11.3629452
## 15 -10.1994528
## 13 -34.2329027
## 9 -22.3998948
## 20 -24.8404527
## 10 -20.6988999
## 8 -5.6362980
## 17 -12.3418381
## 19 -17.4448227
## 4 -6.2774215
## 6 6.5170183
## 11 20.8645404
## 3 8.3659258
## 1 4.6681108
## 12 6.8128435
## 2 11.9897845
## 5 17.9802447
## 14 39.8713094
## 7 38.6880086
## 24 -50.0240410
## 30 -16.9655751
## 22 -16.2999684
## 23 -9.9397267
## 27 -10.8272022
## 28 6.4785718
## 26 7.4400037
## 25 18.1636672
## 29 17.0543227
## 21 54.9199480
## 33 -29.4641897
## 37 -35.8983878
## 36 -7.1293873
## 31 -12.7500660
## 39 -7.7210377
## 38 -0.5472766
## 32 12.9867262
## 40 12.9867262
## 34 22.9708266
## 35 44.5660661
## 44 -15.2858319
## 45 -14.9417416
## 43 5.8399785
## 41 -7.1023739
## 47 -7.5461117
## 49 1.1807317
## 46 -2.0733455
## 50 9.8336187
## 42 11.2387884
## 48 18.8562872fit4coef <- coef(fit4)
fit4coef <- fit4coef$`Chick`
names(fit4coef) <- c("Intercept", "Time", "Diet2", "Diet3", "Diet4")
range(ChickWeight$Time)## [1] 0 21range(ChickWeight$weight)## [1] 35 373plot(range(ChickWeight$Time), range(ChickWeight$weight), type = "n", xlab = "Time", ylab = "Weight")
for (i in 1:50){
abline(a = fit4coef$Intercept[i], b = fit4coef$Time[i])
}
fit5 <- lmer(weight ~ Time + Diet + (Time + 1 | Chick), data=ChickWeight, REML=F)
summary(fit5)## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: weight ~ Time + Diet + (Time + 1 | Chick)
## Data: ChickWeight
##
## AIC BIC logLik deviance df.resid
## 4834.1 4873.3 -2408.0 4816.1 569
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.7927 -0.5834 -0.0388 0.4759 3.5200
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## Chick (Intercept) 147.70 12.153
## Time 13.85 3.721 -0.99
## Residual 163.44 12.784
## Number of obs: 578, groups: Chick, 50
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 26.3563 2.2434 11.748
## Time 8.4439 0.5349 15.787
## Diet2 2.8382 2.2655 1.253
## Diet3 2.0075 2.2655 0.886
## Diet4 9.2547 2.2686 4.079
##
## Correlation of Fixed Effects:
## (Intr) Time Diet2 Diet3
## Time -0.806
## Diet2 -0.344 -0.004
## Diet3 -0.344 -0.004 0.344
## Diet4 -0.343 -0.005 0.344 0.344anova(fit4, fit5)## Data: ChickWeight
## Models:
## fit4: weight ~ Time + Diet + (1 | Chick)
## fit5: weight ~ Time + Diet + (Time + 1 | Chick)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## fit4 7 5619.2 5649.7 -2802.6 5605.2
## fit5 9 4834.1 4873.3 -2408.0 4816.1 789.12 2 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1fit51 <- lmer(weight ~ Time + Diet + (Time + 1 || Chick), data=ChickWeight, REML=F)
summary(fit51)## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: weight ~ Time + Diet + ((1 | Chick) + (0 + Time | Chick))
## Data: ChickWeight
##
## AIC BIC logLik deviance df.resid
## 4900.6 4935.5 -2442.3 4884.6 570
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.5589 -0.5759 0.0104 0.5013 3.5397
##
## Random effects:
## Groups Name Variance Std.Dev.
## Chick (Intercept) 84.13 9.172
## Chick.1 Time 12.11 3.480
## Residual 165.94 12.882
## Number of obs: 578, groups: Chick, 50
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 32.7541 2.6182 12.510
## Time 8.4615 0.5033 16.813
## Diet2 -4.0813 4.4865 -0.910
## Diet3 -13.7193 4.4865 -3.058
## Diet4 -0.5887 4.4937 -0.131
##
## Correlation of Fixed Effects:
## (Intr) Time Diet2 Diet3
## Time -0.061
## Diet2 -0.582 0.006
## Diet3 -0.582 0.006 0.339
## Diet4 -0.581 0.005 0.339 0.339anova(fit5, fit51)## Data: ChickWeight
## Models:
## fit51: weight ~ Time + Diet + ((1 | Chick) + (0 + Time | Chick))
## fit5: weight ~ Time + Diet + (Time + 1 | Chick)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## fit51 8 4900.6 4935.5 -2442.3 4884.6
## fit5 9 4834.1 4873.3 -2408.0 4816.1 68.528 1 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Now, let’s show how different diets work together with time
ranef(fit5)## $Chick
## (Intercept) Time
## 18 4.1668785 -1.29591773
## 16 23.7784392 -7.57079433
## 15 20.2959381 -6.33338499
## 13 20.1593290 -6.34159960
## 9 17.5189125 -5.25438992
## 20 15.5150916 -4.92271512
## 10 14.1831208 -4.45322840
## 8 10.6546572 -3.16300204
## 17 11.9519838 -3.60547220
## 19 11.4638613 -3.71945927
## 4 7.6470366 -2.40438333
## 6 5.1634268 -1.30984247
## 11 0.5975071 0.30603394
## 3 0.2502925 -0.15820763
## 1 1.8840758 -0.67436099
## 12 0.5773555 -0.29665198
## 2 -0.7847380 0.22413806
## 5 -4.5908894 1.31322401
## 14 -11.9945795 3.87404131
## 7 -14.5173423 4.37323644
## 24 22.9989183 -7.13056193
## 30 7.8960691 -2.40777267
## 22 7.8450038 -2.36787675
## 23 5.1473284 -1.50461895
## 27 3.7480449 -1.23771996
## 28 -3.8689533 1.17152819
## 26 -4.8876421 1.43528310
## 25 -9.0862450 2.81446946
## 29 -11.4520848 3.28262842
## 21 -23.5718432 7.56236937
## 33 7.0842582 -1.95919323
## 37 6.2024790 -2.04482065
## 36 -5.2473758 1.69578131
## 31 -4.5149541 1.29317745
## 39 -6.8066294 2.01035700
## 38 -10.7568740 3.18199823
## 32 -15.1777495 4.73108770
## 40 -15.7535025 4.85668299
## 34 -20.6257195 6.35086050
## 35 -28.7170758 9.04890528
## 44 6.1596481 -1.85045878
## 45 2.7194734 -0.90789798
## 43 -1.6723635 0.94798339
## 41 0.5453797 -0.09495143
## 47 0.1356136 -0.02473568
## 49 -4.0817198 1.27226706
## 46 -3.6574757 1.03913470
## 50 -8.8877909 2.69450242
## 42 -10.2004118 3.04154692
## 48 -15.4361630 4.51278076fit5coef <- coef(fit5)
fit5coef <- fit5coef$`Chick`
names(fit5coef) <- c("Intercept", "Time", "Diet2", "Diet3", "Diet4")
fit5coef$Chick <- rownames(fit5coef)
range(ChickWeight$Time)## [1] 0 21range(ChickWeight$weight)## [1] 35 373diets <- unique(data.frame(Chick = ChickWeight$Chick, Diet = ChickWeight$Diet))
fit5coef <- merge(fit5coef, diets, all.x = T)
dummies <- model.matrix(~fit5coef$Diet)
names(dummies) <- c("Intercept", "Diet2d", "Diet3d", "Diet4d")
fit5coef <- cbind(fit5coef, dummies)
names(fit5coef) <- c("Chick", "Intercept", "Time", "Diet2",
"Diet3", "Diet4", "Diet", "(Intercept)",
"Diet2d", "Diet3d", "Diet4d")
plot(range(ChickWeight$Time), range(ChickWeight$weight), type = "n", xlab = "Time", ylab = "Weight")
for (i in 1:50){
abline(a = fit5coef$Intercept[i] + fit5coef$Diet4d[i]* fit5coef$Diet4[i] + fit5coef$Diet3d[i]* fit5coef$Diet3[i] + fit5coef$Diet2d[i]* fit5coef$Diet2[i],
b = fit5coef$Time[i], col = fit5coef$Diet[i])
}
Cross-level interaction
fit6 <- lmer(weight ~ Time * Diet + (1 + Time | Chick), data=ChickWeight, REML=F)
summary(fit6)## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: weight ~ Time * Diet + (1 + Time | Chick)
## Data: ChickWeight
##
## AIC BIC logLik deviance df.resid
## 4824.2 4876.5 -2400.1 4800.2 566
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.7508 -0.5693 -0.0401 0.4694 3.5415
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## Chick (Intercept) 103.61 10.179
## Time 10.01 3.165 -0.99
## Residual 163.36 12.781
## Number of obs: 578, groups: Chick, 50
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 33.6541 2.8023 12.009
## Time 6.2799 0.7303 8.598
## Diet2 -5.0205 4.8072 -1.044
## Diet3 -15.4038 4.8072 -3.204
## Diet4 -1.7475 4.8145 -0.363
## Time:Diet2 2.3293 1.2508 1.862
## Time:Diet3 5.1430 1.2508 4.112
## Time:Diet4 3.2528 1.2515 2.599
##
## Correlation of Fixed Effects:
## (Intr) Time Diet2 Diet3 Diet4 Tm:Dt2 Tm:Dt3
## Time -0.881
## Diet2 -0.583 0.513
## Diet3 -0.583 0.513 0.340
## Diet4 -0.582 0.513 0.339 0.339
## Time:Diet2 0.514 -0.584 -0.882 -0.300 -0.299
## Time:Diet3 0.514 -0.584 -0.300 -0.882 -0.299 0.341
## Time:Diet4 0.514 -0.584 -0.300 -0.300 -0.882 0.341 0.341anova(fit5, fit6)## Data: ChickWeight
## Models:
## fit5: weight ~ Time + Diet + (Time + 1 | Chick)
## fit6: weight ~ Time * Diet + (1 + Time | Chick)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## fit5 9 4834.1 4873.3 -2408.0 4816.1
## fit6 12 4824.2 4876.5 -2400.1 4800.2 15.85 3 0.001217 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1fit6coef <- coef(fit6)
fit6coef <- fit6coef$`Chick`
names(fit6coef) <- c("Intercept", "Time", "Diet2", "Diet3", "Diet4",
"Time:Diet2", "Time:Diet3", "Time:Diet4")
fit6coef$Chick <- rownames(fit6coef)
fit6coef <- merge(fit6coef, diets, all.x = T)
dummies <- model.matrix(~fit6coef$Diet)
names(dummies) <- c("Intercept", "Diet2d", "Diet3d", "Diet4d")
fit6coef <- cbind(fit6coef, dummies[,2:4])
names(fit6coef) <- c("Chick", "Intercept", "Time", "Diet2",
"Diet3", "Diet4", "TimeDiet2",
"TimeDiet3", "TimeDiet4", "Diet",
"Diet2d", "Diet3d", "Diet4d")plot(range(ChickWeight$Time), range(ChickWeight$weight), type = "n", xlab = "Time", ylab = "Weight")
for (i in 1:50){
abline(a = fit6coef$Intercept[i] + fit6coef$Diet4d[i]* fit6coef$Diet4[i] + fit6coef$Diet3d[i]* fit6coef$Diet3[i] + fit6coef$Diet2d[i]* fit6coef$Diet2[i],
b = fit6coef$Time[i] + fit6coef$Diet4d[i]* fit6coef$TimeDiet4[i] + fit6coef$Diet3d[i]* fit6coef$TimeDiet3[i] + fit6coef$Diet2d[i]* fit6coef$TimeDiet2[i],
col = fit6coef$Diet[i])
}
legend(x = 0, y = 250, legend = c("Diet1", "Diet2", "Diet3", "Diet4"), col = 1:4, lty = 1, bty = "n")
ggplot(fortify(fit6), aes(Time, weight, color=Diet)) +
stat_summary(fun.data=mean_se, geom="pointrange") +
stat_summary(aes(y=.fitted), fun.y=mean, geom="line")
plot_model(fit6, type = "int")