beerfoam <- read.table("http://www.stat.ufl.edu/~winner/data/beer_foam2.dat", header=F,col.names=c("foam.time","brand1","brand2","brand3")) attach(beerfoam) foam.time for (i in 1:length(foam.time)) { if (foam.time[i] == 0) foam.time[i] <- 0.0001 } foam.time Y.foam <- c(brand1,brand2,brand3) X1 <- c(rep(1,15),rep(0,15),rep(0,15)) X2 <- c(rep(0,15),rep(1,15),rep(0,15)) X3 <- c(rep(0,15),rep(0,15),rep(1,15)) t.foam <- c(foam.time,foam.time,foam.time) brand <- rep(1:3,each=15) foam.mod1 <- nls(Y.foam ~ b01*X1*exp(-b11*X1*t.foam) + b02*X2*exp(-b12*X2*t.foam) + b03*X3*exp(-b13*X3*t.foam), start=c(b01=10,b02=10,b03=10,b11=0.01,b12=0.01,b13=0.01)) summary(foam.mod1) time.x <- 0:360 yhat.b1 <- coef(foam.mod1)[1] * exp(-coef(foam.mod1)[4]*time.x) yhat.b2 <- coef(foam.mod1)[2] * exp(-coef(foam.mod1)[5]*time.x) yhat.b3 <- coef(foam.mod1)[3] * exp(-coef(foam.mod1)[6]*time.x) plot(t.foam,Y.foam,pch=brand) lines(time.x,yhat.b1,lty=1) lines(time.x,yhat.b2,lty=2) lines(time.x,yhat.b3,lty=5) legend(240,16,c("Erd","Aug","Bud"),pch=c(1,2,3),lty=c(1,2,5)) (beta.foam1 <- coef(foam.mod1)) vcov.foam1 <- vcov(foam.mod1) (mse <- deviance(foam.mod1)/(45-6)) GpGinv <- vcov.foam1/mse Kp <- matrix(c(0,1,-1,0,0,0,0,0,0,0,1,-1),byrow=T,ncol=6) (F.foam1a <- t(Kp %*% beta.foam1) %*% solve(Kp %*% GpGinv %*% t(Kp)) %*% (Kp %*% beta.foam1)) (F.foam1b <- nrow(Kp)*mse) (F.foam1 <- F.foam1a / F.foam1b) (crit.F <- qf(.95,2,39)) (P.F <- 1-pf(F.foam1,2,39)) ### Fit Reduced Model X23 <- X2 + X3 foam.mod2 <- nls(Y.foam ~ b01*X1*exp(-b11*X1*t.foam) + b023*X23*exp(-b123*X23*t.foam), start=c(b01=10,b023=10,b11=0.01,b123=0.01)) summary(foam.mod2) ### Compare Reduced and Complete Models anova(foam.mod2,foam.mod1)