kestrel<-read.table('d:/keith sockman/kestrel_data.csv',header=TRUE,sep=',',na.strings='.')
#restrict to first egg and nonmissing mass
kestrel1<-kestrel[kestrel$egg==1 & !is.na(kestrel$nstlmass),]
#get counts on
number of measurements per nestling
table(kestrel2$nest)
3 4
5 6 8 10 11 12 13 14 15 17 18 21 22 25 27 29 30
2
6 2 5
4 3 5 5 3 5
4 1 6 6 6 6 5 3 5
library(lattice)
trellis.par.set(col.whitebg())
#create
grouped data object for some nice default graphs
kestrel.grp<-groupedData(nstlmass~age|nest,data=kestrel1)
plot(kestrel.grp)
plot(kestrel.grp,outer=~sex*trtmnt)
plot(kestrel.grp,outer=~sex*trtmnt,aspect=.5)
Gompertz Parameterizations
|
Spangler 2005, version 1
|
|
|
Spangler
2005, version 2
|
|
|
Spangler
2005, version 3
|
|
|
Song
& Kuznetsova 2001
|
|
|
Chabot
and Stenson 2002
|
|
|
Lewis et al 2002
|
|
|
Generic
form
|
|
|
SSgompertz function in R 2004.
|
|
To understand how to proceed, fit
separate Gompertz models to each nest.
#get the nest
numbers that have 4 or more observations so can get confidence intervals
names(table(kestrel1$nest)[table(kestrel1$nest)>=4])->good.nests
#choose only these nests
kestrel1.short<-kestrel1[kestrel1$nest
%in% good.nests,]
#discard variables not needed
just.enough<-kestrel1.short[,c("nstlmass","age","nest")]
#load nlme
library
library(nlme)
#create a groupedData
object defining model structure
kestrel.grp<-groupedData(nstlmass~age|nest,data=just.enough)
#fit individual Gompertz
models to each nest
out.nls<-nlsList(nstlmass~SSgompertz(age,a0,b0,b1)|nest,data=kestrel.grp)
#plot parameter estimates
plot(intervals(out.nls))
#to change background to white first load
lattice library
library(lattice)
trellis.par.set(col.whitebg())
#produce plot again
plot(intervals(out.nls))
Interpreting the
parameters
#produce Fig. 2
#Gompertz
function
gomp<-function(x,a0,b0,b1) a0*exp(-b0*b1^x)
#x-values for plotting
xfit<-seq(0,25,.1)
#obtain parameter estimates
apply(coef(out.nls),2,function(x) c(mean(x),max(x),min(x)))->param.list
rownames(param.list)<-c('mean','max','min')
param.list
a0 b0
b1
mean 100.2154 3.506584 0.8184082
max 155.7684 9.960794 0.8909430
min 35.3333 1.242678 0.6472185
PARM graphs.txt online
#graph of a0
plot(xfit,sapply(xfit,function(x) gomp(x,param.list[1,1],
param.list[1,2],param.list[1,3])),type='n',
xlab='Age',ylab='Nestling
Mass',axes=FALSE,cex.lab=.9)
axis(1,cex.axis=.85)
axis(2,cex.axis=.85)
box()
lines(xfit,sapply(xfit,function(x) gomp(x,param.list[1,1],
param.list[1,2],param.list[1,3])),lty=1)
lines(xfit,sapply(xfit,function(x) gomp(x,param.list[1,1],
param.list[1,2],param.list[2,3])),lty=2,col=2)
lines(xfit,sapply(xfit,function(x) gomp(x,param.list[1,1],
param.list[1,2],param.list[3,3])),lty=2,col=4)
legend(12,30,expression(paste("average ", b[1]),
paste("maximum ", b[1]),
paste("minimum ",b[1])),lty=c(1,2,2),col=c(1,2,4),
cex=c(.8,.8,.8),bty='n')
mtext(expression(paste("Effect
of Varying ",b[1],
" on
Shape of Gompertz Curve")),side=3,line=.5,cex=.9)
#graph of b0
plot(xfit,sapply(xfit,function(x) gomp(x,param.list[1,1],
param.list[1,2],param.list[1,3])),type='n',
xlab='Age',ylab='Nestling
Mass',axes=FALSE,cex.lab=.9)
axis(1,cex.axis=.85)
axis(2,cex.axis=.85)
box()
lines(xfit,sapply(xfit,function(x) gomp(x,param.list[1,1],
param.list[1,2],param.list[1,3])),lty=1)
lines(xfit,sapply(xfit,function(x) gomp(x,param.list[1,1],
param.list[2,2],param.list[1,3])),lty=2,col=2)
lines(xfit,sapply(xfit,function(x) gomp(x,param.list[1,1],
param.list[3,2],param.list[1,3])),lty=2,col=4)
legend(15,30,expression(paste("average ", b[0]),
paste("maximum ", b[0]),
paste("minimum ",b[0])),lty=c(1,2,2),col=c(1,2,4),
cex=c(.8,.8,.8),bty='n')
mtext(expression(paste("Effect
of Varying ",b[0],
" on
Shape of Gompertz Curve")),side=3,line=.5,cex=.9)
#graph of b1
plot(xfit,sapply(xfit,function(x) gomp(x,param.list[2,1],
param.list[1,2],param.list[1,3])),type='n',
xlab='Age',ylab='Nestling
Mass',axes=FALSE,cex.lab=.9)
axis(1,cex.axis=.85)
axis(2,cex.axis=.85)
box()
lines(xfit,sapply(xfit,function(x) gomp(x,param.list[1,1],
param.list[1,2],param.list[1,3])),lty=1)
lines(xfit,sapply(xfit,function(x) gomp(x,param.list[2,1],
param.list[1,2],param.list[1,3])),lty=2,col=2)
lines(xfit,sapply(xfit,function(x) gomp(x,param.list[3,1],
param.list[1,2],param.list[1,3])),lty=2,col=4)
legend(0,155,expression(paste("average ",
a[0]),paste("maximum ",
a[0]),
paste("minimum ",a[0])),lty=c(1,2,2),col=c(1,2,4),
cex=c(.8,.8,.8),bty='n')
mtext(expression(paste("Effect
of Varying ",a[0],
" on
Shape of Gompertz Curve")),side=3,line=.5,cex=.9)
Model Fitting
kestrel2<-kestrel[kestrel$egg==1 & !is.na(kestrel$nstlmass) &
!is.na(kestrel$sex),]
> param.list
a0 b0
b1
mean 100.2154 3.506584 0.8184082
max 155.7684 9.960794 0.8909430
min 35.3333 1.242678
0.6472185
common.growth<-nlme(nstlmass~SSgompertz(age,a0,b0,b1),
fixed=a0+b0+b1~1,random=a0~1|nest,data=kestrel2,
start=c(100.2,3.5,.80),na.action=na.omit)
Nonlinear
mixed-effects model fit by maximum likelihood
Model: nstlmass ~ SSgompertz(age,
a0, b0, b1)
Data: kestrel2
AIC
BIC logLik
615.2492 627.2214 -302.6246
Random
effects:
Formula: a0 ~ 1 | nest
a0
Residual
StdDev: 21.94385
7.90911
Fixed
effects: a0 + b0 + b1 ~ 1
Value Std.Error
DF t-value
p-value
a0 100.47474 7.156029 61
14.04057 0
b0 2.73085 0.193302 61 14.12736 0
b1 0.86035 0.011884 61 72.39503 0
Correlation:
a0
b0
b0
-0.171
b1 0.520 -0.715
Standardized
Within-Group Residuals:
Min Q1 Med Q3 Max
-2.4137203
-0.4015826 0.3309273 0.6684498
2.0295954
Number
of Observations: 81
Number
of Groups: 18
Need initial values. Explain
how obtain them for full interaction model.
>
tapply(as.character(kestrel2$trtmnt),kestrel2$nest,function(x)
x[1])->trt.list
>
tapply(as.character(kestrel2$sex),kestrel2$nest,function(x)
x[1])-> gender.list
>
data1<-data.frame(rownames(trt.list),trt.list,gender.list)
>
colnames(data1)[1]<-"nest"
>
data2<-data.frame(rownames(coef(out.nls)),coef(out.nls))
>
colnames(data2)[1]<-"nest"
>
data3<-merge(data1,data2)
>
data3
nest trt.list gender.list a0
b0 b1
1 11 control m
95.06870 2.310242 0.8341338
2 12 control m
79.98548 9.960794 0.6472185
3 14
4 15
5 18 control f 103.98342 4.533675 0.7863282
6 21
7 22 control m 112.46085 2.598779 0.8627828
8 25 control m 113.91226 2.682545 0.8587332
9 27
10 30 control f 121.34002 3.610447 0.8274092
11 4 control f 155.76843 3.637684 0.8722734
12 6 control f 101.11714 4.770902 0.7852193
13 8 control f
58.81220 1.955768 0.8160963
>
tapply(data3$a0,list(data3$gender.list,data3$trt.list),mean)->a0.mean
>
a0.mean
f
125.57757 108.2042
m 54.59836 100.3568
>
tapply(data3$b0,list(data3$gender.list,data3$trt.list),mean)->b0.mean
>
b0.mean
f
3.074510 3.701695
m
1.687869 4.388090
>
tapply(data3$b1,list(data3$gender.list,data3$trt.list),mean)->b1.mean
>
b1.mean
f
0.8781828 0.8174653
m
0.7963734 0.8007171
and
A full-interaction means model
for parameter 
is
,
>
as.vector(a0.mean)
[1] 125.57757 54.59836
108.20424 100.35682
f
Û
> amat<-matrix(c(1,0,0,0,1,1,0,0,1,0,1,0,1,1,1,1),byrow=TRUE,ncol=4)
> amat
[,1] [,2] [,3] [,4]
[1,] 1 0 0 0
[2,] 1 1 0 0
[3,] 1 0 1 0
[4,] 1 1 1 1
> solve(amat,as.vector(a0.mean))
[1] 125.57757 -70.97921 -17.37333 63.13180
|
Level 1: |
|
|
Level 2: |
|
This is specified in R as follows.
a0.int<-nlme(nstlmass~SSgompertz(age,a0,b0,b1),
fixed=list(a0~sex*trtmnt,b0+b1~1),random=a0~1|nest,data=kestrel2,
start=c(125.57,-70.9,-17.3,63.1,2.7,.86),na.action=na.omit)
Observe that the fixed equation for uses the
shortcut notation sex*trtmnt
which
Sample size is small and number of parameters is large. We need AIC
N = 81
AICc<-function(AIC,K)
AIC+2*K*(K+1)/(dim(kestrel2)[1]-K-1)
>
AIC(a0.int)
606.0928
>
AICc(AIC(a0.int),8)
608.0928
Another model
a0.sex<-nlme(nstlmass~SSgompertz(age,a0,b0,b1),fixed=list(a0~sex,b0+b1~1),
random=a0~1|nest,data=kestrel2,start=c(113,-30,3.6,.8),na.action=na.omit)
Get initial
conditions
>
tapply(data3$a0,data3$gender.list,mean)
f m
113.1680 85.1040
>
a0.sex<-nlme(nstlmass~SSgompertz(age,a0,b0,b1),fixed=list(a0~sex,b0+b1~1),
+
random=a0~1|nest,data=kestrel2,start=c(113,-30,3.6,.8),na.action=na.omit)
>
fixef(a0.sex)
a0.(Intercept)
a0.sexm b0 b1
110.5616191 -20.0127043 2.7215707 0.8613758
>
AICc(AIC(a0.sex),6)
614.8257
Table 5 Model Comparison Results
|
|
Level-2
Model |
|
|
|
||
|
|
|
|
|
K |
AIC |
AICc |
|
|
|
|
|
5 |
615.25 |
616.05 |
|
|
|
|
|
6 |
613.69 |
614.83 |
|
|
|
|
|
6 |
610.54 |
611.67 |
|
|
|
|
|
7 |
608.01 |
609.55 |
|
|
|
|
|
8 |
606.09 |
608.09 |
|
A |
|
|
|
9 |
600.03 |
602.57 |
|
|
|
|
|
10 |
601.80 |
604.94 |
|
|
|
|
|
10 |
600.43 |
603.57 |
|
|
|
|
|
10 |
601.86 |
605.00 |
Fit best model
>
a0.int.b0sex<-nlme(nstlmass~SSgompertz(age,a0,b0,b1),fixed=list(a0~sex*trtmnt,b0~sex,b1~1),
random=a0~1|nest,
data=kestrel2,start=c(138.4,-84.9,-30.2,77.1,2.5,-1,.8),na.action=na.omit)
>
AICc(AIC(a0.int.b0sex),9)
602.5653
Examine Fit
final trellis 1.txt
final trellis2.txt
Diagnostics
library(car)
#fit final model
qq.plot(residuals(a0.int.b0sex,type='pearson'),
ylab='Pearson
residuals')
Heteroscedasticity
plot(predict(a0.int.b0sex),residuals(a0.int.b0sex,type='pearson'),
xlab='Predicted
values',ylab='Pearson residuals')
boxplot(residuals(a0.int.b0sex,type='pearson')~ factor(kestrel2$age),
xlab='Age',ylab='Pearson Residuals', col='gray85',axes=FALSE,
outline=FALSE)
points(jitter(as.numeric(factor(kestrel2$age)),amount=.1),
residuals(a0.int.b0sex,
type='pearson'),pch=19,cex=.7,col=4)
#add means to boxplots
points(1:6,tapply(residuals(a0.int.b0sex,type='pearson'),
kestrel2$age,mean),pch=16,col=2)
axis(1,labels=seq(0,25,5))
axis(2)
box()
#interactively
locate outliers
identify(factor(kestrel2$age),residuals(a0.int.b0sex,type='pearson'),
cex=.8,labels=kestrel2
$nest,offset=.75,col=4)
Correlation
ACF.Nfunction(residuals(a0.int.b0sex,type='pearson'),kestrel2$nest,5)->model.N
ACF.function(residuals(a0.int.b0sex,type='pearson'),kestrel2$nest,5)->model.acf
plot(model.acf[,1], model.acf[,2],type='h',xlab='Lag',
ylab='Autocorrelation',ylim=c(-1,1),col=4)
abline(h=0,lty=1)
lines(model.acf[,1],qnorm(.975)/sqrt(model.N),lty=2,col=2)
lines(model.acf[,1],-qnorm(.975)/sqrt(model.N),lty=2,col=2)
mtext('Autocorrelation of residuals',line=.5,side=3,font=2,cex=.9)
tapply(residuals(a0.int.b0sex),kestrel2$age,var)->vars
plot(1:6,vars,xlab='Age',ylab='Residual Variance',axes=FALSE)
axis(1,at=1:6,labels=seq(0,25,5))
axis(2)
box()
lines(lowess(vars~c(1,2,3,4,5,6)),lty=1,col=2)
abline(lm(vars~c(1,2,3,4,5,6)),lty=2,col=4)
legend(1,130,c('lowess','linear regression'),
col=c(2,4),lty=c(1,2),cex=c(.9,.9),bty='n')
