FW599_Day3_MixedModel
R Markdown
1. Streams example
library(nlme)
library(lattice)Read in the data and convert to a groupedData object
TheData <- scan("mixed1a.txt",what=list(Stream=0,Tst=0,Density=0),n=3*18)identical read in of data here to teach scan()
TheData <- scan("mixed1a.txt",what=list(Stream=0,Tst=0,Density=0),nlines=18)
note that "=0" tells R that we are dealing with something numeric
Streams <- data.frame(TheData)
print(Streams)## Stream Tst Density
## 1 1 97.958 95.762
## 2 1 97.958 88.589
## 3 1 97.958 101.885
## 4 2 92.684 85.250
## 5 2 92.684 98.803
## 6 2 92.684 90.473
## 7 3 90.199 94.500
## 8 3 90.199 84.409
## 9 3 90.199 91.775
## 10 4 66.194 57.824
## 11 4 66.194 70.890
## 12 4 66.194 73.666
## 13 5 73.195 72.888
## 14 5 73.195 68.265
## 15 5 73.195 75.438
## 16 6 54.037 50.678
## 17 6 54.037 59.399
## 18 6 54.037 54.823print("Linear Fixed Effects Model")## [1] "Linear Fixed Effects Model"lm1 <- lm(Density~1,data=Streams)
par(mfrow=c(2,2), mar=c(5,5,1,1))
plot(sample(lm1$residuals,18));abline(0,0,lwd=3)# doesn't look so bad, but wait
boxplot(split(lm1$residuals,Streams$Stream),ylab="Residual",xlab="Stream")
abline(0,0,lwd=3)
print(summary(lm1))##
## Call:
## lm(formula = Density ~ 1, data = Streams)
##
## Residuals:
## Min 1Q Median 3Q Max
## -27.951 -9.707 1.295 12.821 23.256
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 78.629 3.778 20.81 1.56e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 16.03 on 17 degrees of freedomprint(paste("AIC ",AIC(lm1)))## [1] "AIC 153.934277919409"print("Linear Fixed Effects Model")## [1] "Linear Fixed Effects Model"lm2 <- lm(Density~factor(Stream)-1,data=Streams)
boxplot(split(lm2$residuals,Streams$Stream),ylab="Residual",xlab="Stream")
abline(0,0,lwd=3)
print(summary(lm2))##
## Call:
## lm(formula = Density ~ factor(Stream) - 1, data = Streams)
##
## Residuals:
## Min 1Q Median 3Q Max
## -9.6360 -4.1995 0.5205 4.0615 7.2943
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## factor(Stream)1 95.412 3.513 27.16 3.82e-12 ***
## factor(Stream)2 91.509 3.513 26.05 6.25e-12 ***
## factor(Stream)3 90.228 3.513 25.69 7.38e-12 ***
## factor(Stream)4 67.460 3.513 19.21 2.24e-10 ***
## factor(Stream)5 72.197 3.513 20.55 1.01e-10 ***
## factor(Stream)6 54.967 3.513 15.65 2.39e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 6.084 on 12 degrees of freedom
## Multiple R-squared: 0.9962, Adjusted R-squared: 0.9942
## F-statistic: 518.8 on 6 and 12 DF, p-value: 8.925e-14print(paste("AIC ",AIC(lm2)))## [1] "AIC 122.786553623177"print("Linear Mixed Effects Model")## [1] "Linear Mixed Effects Model"lm3 <- lme(fixed = Density ~ 1,data=Streams,random = ~ 1 | Stream)
boxplot(split(lm3$residuals,Streams$Stream),ylab="Residual",xlab="Stream")
abline(0,0,lwd=3)
print(summary(lm3))## Linear mixed-effects model fit by REML
## Data: Streams
## AIC BIC logLik
## 133.7973 136.2969 -63.89865
##
## Random effects:
## Formula: ~1 | Stream
## (Intercept) Residual
## StdDev: 15.78877 6.083881
##
## Fixed effects: Density ~ 1
## Value Std.Error DF t-value p-value
## (Intercept) 78.62872 6.60332 12 11.90745 0
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -1.6704314 -0.8239512 0.1256738 0.5295573 1.2987990
##
## Number of Observations: 18
## Number of Groups: 6hint, use coef for "best" estimate of random effects and "intervals" for CI on variance
coef(lm3)nlme::intervals(lm3, 0.95)## Approximate 95% confidence intervals
##
## Fixed effects:
## lower est. upper
## (Intercept) 64.24132 78.62872 93.01612
##
## Random Effects:
## Level: Stream
## lower est. upper
## sd((Intercept)) 8.23611 15.78877 30.26734
##
## Within-group standard error:
## lower est. upper
## 4.077843 6.083881 9.0767602. Relationship between weight and length example
Read in the data and convert to a groupedData object
TheData <- scan("mixed1b.txt",what=list(Subject=0,TST1=0,TST2=0,Length=0,Weight=0,TST3=0),n=6*100)
LenW <- data.frame(TheData)
LenW$Subject=as.factor(LenW$Subject)
par(mfrow=c(1,2))
plot(Weight~Length, data=LenW)
plot(Weight~Length, data=subset(LenW, Subject==1),pch=2)
points(Weight~Length, data=subset(LenW, Subject==2),pch=3)
points(Weight~Length, data=subset(LenW, Subject==3),pch=4)
points(Weight~Length, data=subset(LenW, Subject==4),pch=5)
points(Weight~Length, data=subset(LenW, Subject==5),pch=6)
points(Weight~Length, data=subset(LenW, Subject==5),pch=7)
RANDOM INTERCEPT AND REML
lm1 <- lme(Weight~Length,data=LenW,random = ~1 |Subject,method="REML")
plot(lm1,Weight~fitted(.)|Subject, abline = c(0,1))
print("Model 1")## [1] "Model 1"summary(lm1)## Linear mixed-effects model fit by REML
## Data: LenW
## AIC BIC logLik
## 6.72662 17.06649 0.63669
##
## Random effects:
## Formula: ~1 | Subject
## (Intercept) Residual
## StdDev: 0.2113517 0.205754
##
## Fixed effects: Weight ~ Length
## Value Std.Error DF t-value p-value
## (Intercept) -0.0266742 0.09560888 89 -0.27899 0.7809
## Length 3.0171512 0.02958709 89 101.97527 0.0000
## Correlation:
## (Intr)
## Length -0.682
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -2.41911627 -0.57997864 -0.05204277 0.57190963 2.32084355
##
## Number of Observations: 100
## Number of Groups: 10intervals(lm1)## Approximate 95% confidence intervals
##
## Fixed effects:
## lower est. upper
## (Intercept) -0.216647 -0.02667423 0.1632986
## Length 2.958362 3.01715121 3.0759401
##
## Random Effects:
## Level: Subject
## lower est. upper
## sd((Intercept)) 0.1274321 0.2113517 0.3505361
##
## Within-group standard error:
## lower est. upper
## 0.1776430 0.2057540 0.2383134coef(lm1)##################################### # ON YOUR OWN TIME IF YOU WANT TO FIT NONLINEAR MIXED MODELS #####################################
#3. Length at age data. Nonlinear fit of Von Bertalanffy growth with and without random effect on L-infinity
##################################### #####################################
library(nlme)
library(graphics)
FileName <- "week6b.txt"
TheData <- scan(FileName,what=list(Subject=0,Age=0,Length=0,NULL),skip=1,n=4*100)
xx <- as.data.frame(cbind(Subject=TheData$Subject,Age=TheData$Age,Length=TheData$Length))
AgeLen <- groupedData(Length~Age|Subject,data=as.data.frame(xx),labels=list(x="Age",y="Length"),units=list(x="(yr)",y="(cm)"))
print(AgeLen)## Grouped Data: Length ~ Age | Subject
## Subject Age Length
## 1 1 4.903 54.131
## 2 1 9.996 75.290
## 3 1 3.922 47.212
## 4 1 2.915 38.468
## 5 1 6.351 62.542
## 6 1 2.382 33.122
## 7 1 6.348 62.532
## 8 1 6.598 63.705
## 9 1 7.851 68.932
## 10 1 8.113 69.728
## 11 2 2.832 40.777
## 12 2 7.794 74.026
## 13 2 7.703 73.966
## 14 2 1.687 27.155
## 15 2 6.619 68.902
## 16 2 2.467 36.542
## 17 2 3.987 51.649
## 18 2 2.015 31.243
## 19 2 4.827 58.135
## 20 2 9.319 79.450
## 21 3 1.901 35.221
## 22 3 2.494 43.375
## 23 3 4.209 62.764
## 24 3 8.969 92.024
## 25 3 4.268 63.403
## 26 3 9.919 95.324
## 27 3 4.924 69.250
## 28 3 8.826 91.455
## 29 3 1.189 23.392
## 30 3 6.083 77.767
## 31 4 9.311 85.416
## 32 4 3.258 48.285
## 33 4 7.878 80.316
## 34 4 8.129 81.154
## 35 4 3.450 50.291
## 36 4 6.619 74.198
## 37 4 3.217 48.093
## 38 4 6.093 71.099
## 39 4 2.429 38.961
## 40 4 3.065 46.288
## 41 5 9.371 75.962
## 42 5 1.845 27.736
## 43 5 2.648 36.984
## 44 5 4.436 52.803
## 45 5 8.476 73.429
## 46 5 8.304 72.820
## 47 5 7.670 70.541
## 48 5 7.518 69.841
## 49 5 5.963 62.537
## 50 5 3.515 45.247
## 51 6 9.689 84.446
## 52 6 6.551 71.911
## 53 6 7.473 76.593
## 54 6 3.532 50.075
## 55 6 6.481 71.795
## 56 6 2.612 40.092
## 57 6 9.217 82.942
## 58 6 9.593 84.287
## 59 6 7.928 78.614
## 60 6 4.884 61.778
## 61 7 5.708 69.809
## 62 7 5.614 69.172
## 63 7 3.785 54.617
## 64 7 6.505 74.862
## 65 7 8.994 85.618
## 66 7 5.746 70.142
## 67 7 2.173 35.952
## 68 7 8.228 82.839
## 69 7 4.430 60.091
## 70 7 7.054 77.657
## 71 8 4.068 60.608
## 72 8 5.491 72.706
## 73 8 6.142 76.985
## 74 8 9.967 94.134
## 75 8 5.032 69.261
## 76 8 9.910 94.080
## 77 8 5.629 73.679
## 78 8 3.514 55.122
## 79 8 8.033 87.116
## 80 8 8.162 87.682
## 81 9 3.958 49.271
## 82 9 3.758 47.395
## 83 9 7.276 69.011
## 84 9 1.092 17.536
## 85 9 3.747 47.414
## 86 9 1.006 16.450
## 87 9 1.508 23.506
## 88 9 6.342 64.464
## 89 9 8.019 71.616
## 90 9 5.550 60.159
## 91 10 5.383 57.338
## 92 10 1.060 16.513
## 93 10 2.460 33.829
## 94 10 2.042 29.203
## 95 10 3.291 42.074
## 96 10 1.893 27.429
## 97 10 9.812 75.047
## 98 10 9.697 74.578
## 99 10 1.824 26.741
## 100 10 4.462 51.614 par(pch=16)
plot(AgeLen)
plot(AgeLen,outer=~1)
Non-linear model
par(mfrow=c(1,1))
lm1 <- nls(formula=Length~Linf*(1-exp(-1*Kappa*(Age-Tzero))),data=AgeLen,start=c(Linf=100,Kappa=0.2,Tzero=0))
print(summary(lm1))##
## Formula: Length ~ Linf * (1 - exp(-1 * Kappa * (Age - Tzero)))
##
## Parameters:
## Estimate Std. Error t value Pr(>|t|)
## Linf 97.79915 4.58527 21.329 < 2e-16 ***
## Kappa 0.20258 0.02558 7.919 3.99e-12 ***
## Tzero 0.08890 0.18218 0.488 0.627
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.417 on 97 degrees of freedom
##
## Number of iterations to convergence: 2
## Achieved convergence tolerance: 4.382e-06 print(AIC(lm1))## [1] 626.6415 plot(lm1,pch=16,csi=0.3)
boxplot(split(residuals(lm1),AgeLen$Subject),ylab="Residual",xlab="Subject",csi=0.2)
abline(0,0,lty=1,lwd=4)
Non-linear mixed model
print("non linear model")## [1] "non linear model" lm2 <- nlme(model=Length~Linf*(1-exp(-1*Kappa*(Age-Tzero))),data=AgeLen,random=Linf~1,fixed=Linf+Kappa+Tzero~1,start=c(Linf=100,Kappa=0.2,Tzero=0))
#print(summary(lm2))
LinfEst <- lm2$coefficients$fixed[1]
#print(LinfEst+lm2$coefficients$random$Subject)
plot(lm2,pch=16,csi=0.3)
plot(lm2,Subject~resid(.),abline=0)
plot(LinfEst+lm2$coefficients$random$Subject,pch=16,cex=1.5)
plot(lm2, form = resid(., type = "p") ~ fitted(.) | Subject, abline = 0,pch=16)
plot(augPred(lm2),csi=0.3,pch=16,lwd=2)