Package install

Packages <- c("tidyverse", "car", "sme", "lattice", "lme4", "splines", "ggplot2", "visreg")
lapply(Packages, library, character.only = TRUE)

Data import

d1=read.csv("female rapa weight.csv")

Data structure

d1$id=as.factor(d1$subject)
d1$weight=as.numeric(d1$weight)
d1$group=as.factor(d1$group)
str(d1)
## 'data.frame':    280 obs. of  5 variables:
##  $ subject: int  1 2 3 4 5 6 7 8 9 10 ...
##  $ group  : Factor w/ 3 levels "0","3","5": 1 1 1 1 1 1 1 1 1 1 ...
##  $ age    : int  14 14 14 14 14 14 14 14 14 14 ...
##  $ weight : num  7.2 7.2 7.6 7 7.4 7.7 7.3 7.5 6.8 7.1 ...
##  $ id     : Factor w/ 35 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 9 10 ...
head(d1)

Explorative data analysis with graphics

Trellis plot using lattice package

print(xyplot(weight~age|id, data=d1,xlab="age", ylab="weight"))

Spaghetti plot

 print(xyplot(weight~age|group,type="l",d1, groups=as.factor(subject), strip=strip.custom(bg="white")))

Spline linear modeling

Spline degree of freedom is 1

par(mfrow=c(2,2))
mdA1<-lm(weight ~ bs(age,knots = c(18,22,28),degree=1),data = d1[d1$group=="0",] )
mdB1<-lm(weight ~ bs(age,knots = c(18,22,28),degree=1),data = d1[d1$group=="3",] )
mdC1<-lm(weight ~ bs(age,knots = c(18,22,28),degree=1),data = d1[d1$group=="5",] )

Spline degree of freedom is 2

par(mfrow=c(2,2))
mdA2<-lm(weight ~ bs(age,knots = c(18,22,28),degree=2),data = d1[d1$group=="0",] )
mdB2<-lm(weight ~ bs(age,knots = c(18,22,28),degree=2),data = d1[d1$group=="3",] )
mdC2<-lm(weight ~ bs(age,knots = c(18,22,28),degree=2),data = d1[d1$group=="5",] )

Spline degree of freedom is 3

par(mfrow=c(2,2))
mdA3<-lm(weight ~ bs(age,knots = c(18,22,28),degree=3),data = d1[d1$group=="0",] )
mdB3<-lm(weight ~ bs(age,knots = c(18,22,28),degree=3),data = d1[d1$group=="3",] )
mdC3<-lm(weight ~ bs(age,knots = c(18,22,28),degree=3),data = d1[d1$group=="5",] )

fitted line plotting

Plotting regression line with df=1

N=seq(14,35, by=0.1)
nd=data.frame(age=N)
plot(jitter(d1$age), d1$weight, col=c("blue", "red", "green")[d1$group], pch=20)
Ya1=predict(mdA1, nd)
Yb1=predict(mdB1, nd)
Yc1=predict(mdC1, nd)
lines(N, Ya1, col="blue")
lines(N, Yb1, col="red")
lines(N, Yc1, col="green")
legend(7,200,legend=c("0","3","5"), col=c("blue","red","green"), lty=c(1,1,1), ncol=1)

Plotting regression line with df=2

plot(jitter(d1$age), d1$weight, col=c("blue", "red", "green")[d1$group], pch=20)
Ya2=predict(mdA2, nd)
Yb2=predict(mdB2, nd)
Yc2=predict(mdC2, nd)
lines(N, Ya2, col="blue")
lines(N, Yb2, col="red")
lines(N, Yc2, col="green")
legend(7,200,legend=c("0","3","5"), col=c("blue","red","green"), lty=c(1,1,1), ncol=1)

Plotting regression line with df=3

plot(jitter(d1$age), d1$weight, col=c("blue", "red", "green")[d1$group], pch=20)
Ya3=predict(mdA3, nd)
Yb3=predict(mdB3, nd)
Yc3=predict(mdC3, nd)
lines(N, Ya3, col="blue")
lines(N, Yb3, col="red")
lines(N, Yc3, col="green")
legend(7,200,legend=c("0","3","5"), col=c("blue","red","green"), lty=c(1,1,1), ncol=1)

그래프 모양으로 보면 자유도 2인 그래프가 데이터에 잘 맞아 보입니다.

모델 요약

md1=lm(weight ~ group*bs(age,knots = c(18,22,28),degree=1),data = d1)
md2=lm(weight ~ group*bs(age,knots = c(18,22,28),degree=2),data = d1)
md3=lm(weight ~ group*bs(age,knots = c(18,22,28),degree=3),data = d1)
summary(md1);summary(md2);summary(md3)
## 
## Call:
## lm(formula = weight ~ group * bs(age, knots = c(18, 22, 28), 
##     degree = 1), data = d1)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.82066 -0.34570 -0.01538  0.36026  2.12500 
## 
## Coefficients:
##                                                    Estimate Std. Error
## (Intercept)                                          7.4081     0.1436
## group3                                              -0.4323     0.2030
## group5                                              -0.3304     0.2076
## bs(age, knots = c(18, 22, 28), degree = 1)1          0.7966     0.2932
## bs(age, knots = c(18, 22, 28), degree = 1)2          2.1734     0.2105
## bs(age, knots = c(18, 22, 28), degree = 1)3          6.9190     0.2202
## bs(age, knots = c(18, 22, 28), degree = 1)4         10.7169     0.2309
## group3:bs(age, knots = c(18, 22, 28), degree = 1)1  -0.6863     0.4147
## group5:bs(age, knots = c(18, 22, 28), degree = 1)1  -0.1729     0.4240
## group3:bs(age, knots = c(18, 22, 28), degree = 1)2  -1.1168     0.2977
## group5:bs(age, knots = c(18, 22, 28), degree = 1)2  -0.8089     0.3044
## group3:bs(age, knots = c(18, 22, 28), degree = 1)3  -1.4379     0.3114
## group5:bs(age, knots = c(18, 22, 28), degree = 1)3  -1.3760     0.3184
## group3:bs(age, knots = c(18, 22, 28), degree = 1)4  -1.6177     0.3265
## group5:bs(age, knots = c(18, 22, 28), degree = 1)4  -1.4492     0.3339
##                                                    t value Pr(>|t|)    
## (Intercept)                                         51.603  < 2e-16 ***
## group3                                              -2.129 0.034151 *  
## group5                                              -1.591 0.112705    
## bs(age, knots = c(18, 22, 28), degree = 1)1          2.717 0.007028 ** 
## bs(age, knots = c(18, 22, 28), degree = 1)2         10.324  < 2e-16 ***
## bs(age, knots = c(18, 22, 28), degree = 1)3         31.426  < 2e-16 ***
## bs(age, knots = c(18, 22, 28), degree = 1)4         46.417  < 2e-16 ***
## group3:bs(age, knots = c(18, 22, 28), degree = 1)1  -1.655 0.099134 .  
## group5:bs(age, knots = c(18, 22, 28), degree = 1)1  -0.408 0.683809    
## group3:bs(age, knots = c(18, 22, 28), degree = 1)2  -3.751 0.000216 ***
## group5:bs(age, knots = c(18, 22, 28), degree = 1)2  -2.657 0.008353 ** 
## group3:bs(age, knots = c(18, 22, 28), degree = 1)3  -4.618 6.05e-06 ***
## group5:bs(age, knots = c(18, 22, 28), degree = 1)3  -4.322 2.19e-05 ***
## group3:bs(age, knots = c(18, 22, 28), degree = 1)4  -4.954 1.29e-06 ***
## group5:bs(age, knots = c(18, 22, 28), degree = 1)4  -4.341 2.02e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.6264 on 265 degrees of freedom
## Multiple R-squared:  0.9677, Adjusted R-squared:  0.966 
## F-statistic: 567.6 on 14 and 265 DF,  p-value: < 2.2e-16
## 
## Call:
## lm(formula = weight ~ group * bs(age, knots = c(18, 22, 28), 
##     degree = 2), data = d1)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.94468 -0.33188 -0.01556  0.33454  2.12500 
## 
## Coefficients:
##                                                    Estimate Std. Error
## (Intercept)                                          7.3171     0.1711
## group3                                              -0.3776     0.2419
## group5                                              -0.2452     0.2473
## bs(age, knots = c(18, 22, 28), degree = 2)1          0.7575     0.3865
## bs(age, knots = c(18, 22, 28), degree = 2)2          1.0965     0.2829
## bs(age, knots = c(18, 22, 28), degree = 2)3          4.3474     0.3648
## bs(age, knots = c(18, 22, 28), degree = 2)4         10.1965     0.5987
## bs(age, knots = c(18, 22, 28), degree = 2)5         10.8079     0.2489
## group3:bs(age, knots = c(18, 22, 28), degree = 2)1  -0.5758     0.5466
## group5:bs(age, knots = c(18, 22, 28), degree = 2)1  -0.3928     0.5589
## group3:bs(age, knots = c(18, 22, 28), degree = 2)2  -0.8265     0.4000
## group5:bs(age, knots = c(18, 22, 28), degree = 2)2  -0.2743     0.4090
## group3:bs(age, knots = c(18, 22, 28), degree = 2)3  -1.6207     0.5159
## group5:bs(age, knots = c(18, 22, 28), degree = 2)3  -1.6222     0.5275
## group3:bs(age, knots = c(18, 22, 28), degree = 2)4  -1.2014     0.8467
## group5:bs(age, knots = c(18, 22, 28), degree = 2)4  -1.0849     0.8658
## group3:bs(age, knots = c(18, 22, 28), degree = 2)5  -1.6724     0.3520
## group5:bs(age, knots = c(18, 22, 28), degree = 2)5  -1.5343     0.3599
##                                                    t value Pr(>|t|)    
## (Intercept)                                         42.776  < 2e-16 ***
## group3                                              -1.561 0.119731    
## group5                                              -0.991 0.322407    
## bs(age, knots = c(18, 22, 28), degree = 2)1          1.960 0.051071 .  
## bs(age, knots = c(18, 22, 28), degree = 2)2          3.876 0.000134 ***
## bs(age, knots = c(18, 22, 28), degree = 2)3         11.916  < 2e-16 ***
## bs(age, knots = c(18, 22, 28), degree = 2)4         17.030  < 2e-16 ***
## bs(age, knots = c(18, 22, 28), degree = 2)5         43.423  < 2e-16 ***
## group3:bs(age, knots = c(18, 22, 28), degree = 2)1  -1.054 0.293078    
## group5:bs(age, knots = c(18, 22, 28), degree = 2)1  -0.703 0.482751    
## group3:bs(age, knots = c(18, 22, 28), degree = 2)2  -2.066 0.039814 *  
## group5:bs(age, knots = c(18, 22, 28), degree = 2)2  -0.671 0.503072    
## group3:bs(age, knots = c(18, 22, 28), degree = 2)3  -3.141 0.001875 ** 
## group5:bs(age, knots = c(18, 22, 28), degree = 2)3  -3.075 0.002328 ** 
## group3:bs(age, knots = c(18, 22, 28), degree = 2)4  -1.419 0.157139    
## group5:bs(age, knots = c(18, 22, 28), degree = 2)4  -1.253 0.211260    
## group3:bs(age, knots = c(18, 22, 28), degree = 2)5  -4.751 3.34e-06 ***
## group5:bs(age, knots = c(18, 22, 28), degree = 2)5  -4.263 2.82e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.6263 on 262 degrees of freedom
## Multiple R-squared:  0.9681, Adjusted R-squared:  0.966 
## F-statistic: 467.7 on 17 and 262 DF,  p-value: < 2.2e-16
## 
## Call:
## lm(formula = weight ~ group * bs(age, knots = c(18, 22, 28), 
##     degree = 3), data = d1)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.92742 -0.35204 -0.00206  0.32502  2.12500 
## 
## Coefficients:
##                                                    Estimate Std. Error
## (Intercept)                                         7.33043    0.17693
## group3                                             -0.38321    0.25022
## group5                                             -0.23949    0.25584
## bs(age, knots = c(18, 22, 28), degree = 3)1         0.50073    0.39752
## bs(age, knots = c(18, 22, 28), degree = 3)2         0.87142    0.46682
## bs(age, knots = c(18, 22, 28), degree = 3)3         1.99222    0.56436
## bs(age, knots = c(18, 22, 28), degree = 3)4         8.21427    1.12322
## bs(age, knots = c(18, 22, 28), degree = 3)5         9.08121    2.73339
## bs(age, knots = c(18, 22, 28), degree = 3)6        10.79457    0.25337
## group3:bs(age, knots = c(18, 22, 28), degree = 3)1 -0.38502    0.56218
## group5:bs(age, knots = c(18, 22, 28), degree = 3)1 -0.38070    0.57482
## group3:bs(age, knots = c(18, 22, 28), degree = 3)2 -0.73808    0.66018
## group5:bs(age, knots = c(18, 22, 28), degree = 3)2 -0.07895    0.67502
## group3:bs(age, knots = c(18, 22, 28), degree = 3)3 -1.09598    0.79813
## group5:bs(age, knots = c(18, 22, 28), degree = 3)3 -0.93863    0.81607
## group3:bs(age, knots = c(18, 22, 28), degree = 3)4 -2.16678    1.58847
## group5:bs(age, knots = c(18, 22, 28), degree = 3)4 -2.14478    1.62417
## group3:bs(age, knots = c(18, 22, 28), degree = 3)5  0.24008    3.86560
## group5:bs(age, knots = c(18, 22, 28), degree = 3)5  0.17628    3.95247
## group3:bs(age, knots = c(18, 22, 28), degree = 3)6 -1.66679    0.35833
## group5:bs(age, knots = c(18, 22, 28), degree = 3)6 -1.54006    0.36638
##                                                    t value Pr(>|t|)    
## (Intercept)                                         41.431  < 2e-16 ***
## group3                                              -1.531 0.126872    
## group5                                              -0.936 0.350112    
## bs(age, knots = c(18, 22, 28), degree = 3)1          1.260 0.208937    
## bs(age, knots = c(18, 22, 28), degree = 3)2          1.867 0.063071 .  
## bs(age, knots = c(18, 22, 28), degree = 3)3          3.530 0.000492 ***
## bs(age, knots = c(18, 22, 28), degree = 3)4          7.313 3.26e-12 ***
## bs(age, knots = c(18, 22, 28), degree = 3)5          3.322 0.001021 ** 
## bs(age, knots = c(18, 22, 28), degree = 3)6         42.603  < 2e-16 ***
## group3:bs(age, knots = c(18, 22, 28), degree = 3)1  -0.685 0.494042    
## group5:bs(age, knots = c(18, 22, 28), degree = 3)1  -0.662 0.508369    
## group3:bs(age, knots = c(18, 22, 28), degree = 3)2  -1.118 0.264606    
## group5:bs(age, knots = c(18, 22, 28), degree = 3)2  -0.117 0.906988    
## group3:bs(age, knots = c(18, 22, 28), degree = 3)3  -1.373 0.170883    
## group5:bs(age, knots = c(18, 22, 28), degree = 3)3  -1.150 0.251131    
## group3:bs(age, knots = c(18, 22, 28), degree = 3)4  -1.364 0.173730    
## group5:bs(age, knots = c(18, 22, 28), degree = 3)4  -1.321 0.187822    
## group3:bs(age, knots = c(18, 22, 28), degree = 3)5   0.062 0.950525    
## group5:bs(age, knots = c(18, 22, 28), degree = 3)5   0.045 0.964460    
## group3:bs(age, knots = c(18, 22, 28), degree = 3)6  -4.652 5.26e-06 ***
## group5:bs(age, knots = c(18, 22, 28), degree = 3)6  -4.203 3.62e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.6283 on 259 degrees of freedom
## Multiple R-squared:  0.9683, Adjusted R-squared:  0.9658 
## F-statistic: 395.2 on 20 and 259 DF,  p-value: < 2.2e-16

모델 세가지를 비교하는 방법은 아직 찾지 못했고 결정계수만 놓고 본다면 md1과 md2가 동일하나 그래프를 보면 md2이 그나마 가장 나은듯

Mixed effects modeling with lme4 packages

Spline degree 1

lmm1=lmer(weight~group*bs(age,knots = c(18,22,28) ,degree=1)+ (1|group/subject), data = d1, na.action=na.exclude, REML=F)
## singular fit
lmm11=lmer(weight~group*bs(age,knots = c(18,22,28) ,degree=1)+ (1+age|group/subject), data = d1, na.action=na.exclude, REML=F)
## singular fit
lmm12=lmer(weight~group*bs(age,knots = c(18,22,28) ,degree=1)+ (1|group/subject) + (0+age|group/subject), data = d1, na.action=na.exclude, REML=F)
## singular fit
anova(lmm1,lmm11,lmm12)

Spline degree 2

lmm2=lmer(weight~group*bs(age,knots = c(18,22,28) ,degree=2)+ (1|group/subject), data = d1, na.action=na.exclude, REML=F)
## singular fit
lmm21=lmer(weight~group*bs(age,knots = c(18,22,28) ,degree=2)+ (1+age|group/subject), data = d1, na.action=na.exclude, REML=F)
## singular fit
lmm22=lmer(weight~group*bs(age,knots = c(18,22,28) ,degree=2)+ (1|group/subject) + (0+age|group/subject), data = d1, na.action=na.exclude, REML=F)
## singular fit
anova(lmm2,lmm21,lmm22)

Random effects plotting

기대한 그래프가 안나옵니다. 모델에 문제가 있는데 아직 찾아내지 못했습니다. ### degree 1 #### Group 0

d1$pred.lmm12=fitted(lmm12)
ggplot(data=subset(d1, group=="0")) + geom_point(aes(x=age,y=weight), color="red", size=1) +   geom_line(aes(x=age,y=pred.lmm12), col="blue") + facet_wrap(~subject, ncol=4)

Group 3

d1$pred.lmm12=fitted(lmm12)
ggplot(data=subset(d1, group=="3")) + geom_point(aes(x=age,y=weight), color="red", size=1) +   geom_line(aes(x=age,y=pred.lmm12), col="blue") + facet_wrap(~subject, ncol=4)

Group 5

d1$pred.lmm12=fitted(lmm12)
ggplot(data=subset(d1, group=="5")) + geom_point(aes(x=age,y=weight), color="red", size=1) +   geom_line(aes(x=age,y=pred.lmm12), col="blue") + facet_wrap(~subject, ncol=4)

degree 2

Group 0

d1$pred.lmm22=fitted(lmm22)
ggplot(data=subset(d1, group=="0")) + geom_point(aes(x=age,y=weight), color="red", size=1) +   geom_line(aes(x=age,y=pred.lmm22), col="blue") + facet_wrap(~subject, ncol=4)

Group 3

d1$pred.lmm22=fitted(lmm22)
ggplot(data=subset(d1, group=="3")) + geom_point(aes(x=age,y=weight), color="red", size=1) +   geom_line(aes(x=age,y=pred.lmm22), col="blue") + facet_wrap(~subject, ncol=4)

Group 5

d1$pred.lmm22=fitted(lmm22)
ggplot(data=subset(d1, group=="5")) + geom_point(aes(x=age,y=weight), color="red", size=1) +   geom_line(aes(x=age,y=pred.lmm22), col="blue") + facet_wrap(~subject, ncol=4)