d=read.delim("http://www.public.iastate.edu/~dnett/S511/IonLeakage.txt") d d\$geno=as.factor(d\$geno) plot(d\$temp,d\$y,col=4,pch=19, xlab="Temperature (degrees Celsius)", ylab="Measure of Ion Leakage") #Fit the 4-parameter Gompertz model #to the data. Initially, we will ignore #the genotype factor. o=nls(y~b1+b2*exp(-exp(-((temp-b3)/b4))), data=d, start=c(b1=7,b2=55,b3=-5,b4=-2)) summary(o) #By default, a numerical algorithm is used #to compute the derivatives needed for #the Gauss-Newton algorithm. In some cases, #there is an advantage to providing the #derivatives along with the function. #This can be done using the deriv function. f=deriv(y~b1+b2*exp(-exp(-((x-b3)/b4))), c("b1","b2","b3","b4"), function(b1,b2,b3,b4,x){}) f o=nls(y~f(b1,b2,b3,b4,temp), data=d, start=c(b1=6,b2=55,b3=-5,b4=-2)) summary(o) #Store the least squares estimator of beta. b=coef(o) b #Add the estimated mean curve to the plot. x=seq(-16,6,by=.1) lines(x,f(b[1],b[2],b[3],b[4],x)) #Examine a plot of standardized #residuals versus fitted values. plot(o) #Get an estimate of the variance of #for the squares estimator bhat: #MSE(Dhat'Dhat)^{-1} v=vcov(o) round(v,3) #Obtain the error degrees of freedom. summary(o)\$df df=summary(o)\$df[2] df #Form a 95% confidence interval for beta_3. b[3]-qt(.975,df)*sqrt(v[3,3]) b[3]+qt(.975,df)*sqrt(v[3,3]) #The following function can be used to obtain #approximate 100(1-a)% confidence intervals #for each element of C%*%beta. ci=function(nlsout,C,a=0.05) { b=coef(nlsout) V=vcov(nlsout) df=summary(nlsout)\$df[2] Cb=C%*%b se=sqrt(diag(C%*%V%*%t(C))) tval=qt(1-a/2,df) low=Cb-tval*se up=Cb+tval*se m=cbind(C,Cb,se,low,up) dimnames(m)[[2]]=c(paste("c",1:ncol(C),sep=""), "estimate","se", paste(100*(1-a),"% Conf.",sep=""), "limits") m } ci(o,matrix(c(1,-1,0,0, 0,0,1,0),byrow=T,nrow=2)) #We can also use the profile function #to invert the reduced versus full model #F test to obtain an approximate 100(1-alpha)% #confidence interval for any one of the #components of beta. op=profile(o,which=3,alpha=0.05) b3=as.matrix(op\$b3)[,c(1,4)] b3 cbind(pt(b3[,1],df),b3) #An even simpler way to do this is to use #the function confint, which gives less #information but produces the profile-based #confidence intervals. confint(o) #Now suppose we want to estimate the #temperature at which the measure of #ion leakage is 50. x50=deriv(y~-b4*log(-log((50-b1)/b2))+b3, c("b1","b2","b3","b4"), function(b1,b2,b3,b4){}) x50 x50(b[1],b[2],b[3],b[4]) est=as.numeric(x50(b[1],b[2],b[3],b[4])) est #Use the delta method to obtain an #approximate standard error associated #with the estimate. der=attr(x50(b[1],b[2],b[3],b[4]),"gradient") der se=drop(sqrt(der%*%v%*%t(der))) se #Obtain an approximate 95% #confidence interval for the temperature #at which the measure of ion leakage is 50. est-qt(.975,df)*se est+qt(.975,df)*se #Now let's consider the genotype variable. plot(d\$temp,d\$y,col=as.numeric(d\$geno),pch=as.numeric(d\$geno), xlab="Temperature (degrees Celsius)", ylab="Measure of Ion Leakage") legend("topright",c("Genotype 1","Genotype 2"),col=1:2,pch=1:2) #Let's fit a model that allows a separate #4-parameter Gompertz model for each genotype. ogeno=nls(y~b1[geno]+b2[geno]*exp(-exp(-((temp-b3[geno])/b4[geno]))), data=d, start=list(b1=c(7,7),b2=c(55,55),b3=c(-5,-5),b4=c(-2,-2))) summary(ogeno) b=coef(ogeno) b #Add the estimated mean curves to the plot. lines(x,f(b[1],b[3],b[5],b[7],x),col=1) lines(x,f(b[2],b[4],b[6],b[8],x),col=2) #Compare reduced versus full model. anova(o,ogeno) #Test for a lack of fit by comparing to #a cell means model with different means #for each combination of genotype and temperature ocellmeans=lm(y~geno*as.factor(temp),data=d) anova(ofull,ocellmeans)