pdf("F:\\Rmisc\\graphs\\salmonella.pdf") salmonella <- read.table("http://www.stat.ufl.edu/~winner/data/salmonella.dat",header=F, col.names=c("expt","dose","resp")) attach(salmonella) for (i in 1:length(dose)) {if(dose[i] == 0) dose[i] <- 0.000001} expt01 <- expt-1 (vardoseexp <- aggregate(resp,by=list(dose,expt),FUN=var)) (ndoseexp <- aggregate(resp,by=list(dose,expt),FUN=length)) (varwt <- rep(1/vardoseexp[[3]],each=ndoseexp[[3]])) salmonella <- data.frame(salmonella,varwt) expt1 <- subset(salmonella,expt==1) salm.mod1 <- nls(resp ~ (b0+exp(b1+b2*log(dose)))*exp(-b3*dose), start=c(b0=20,b1=10,b2=1,b3=1),weight=varwt,data=expt1) summary(salm.mod1) deviance(salm.mod1) expt2 <- subset(salmonella,expt==2) salm.mod2 <- nls(resp ~ (b0+exp(b1+b2*log(dose)))*exp(-b3*dose), start=c(b0=20,b1=10,b2=1,b3=1),weight=varwt,data=expt2) summary(salm.mod2) deviance(salm.mod2) salm.mod3 <- nls(resp ~ (b0+exp(b1+b2*log(dose)))*exp(-b3*dose), start=c(b0=20,b1=10,b2=1,b3=1),weight=varwt,data=salmonella) summary(salm.mod3) deviance(salm.mod3) salm.mod4 <- nls(resp ~ expt01*(b02+exp(b12+b22*log(dose)))*exp(-b32*dose) + (1-expt01)*(b01+exp(b11+b21*log(dose)))*exp(-b31*dose), start=c(b01=20,b11=10,b21=1,b31=1,b02=20,b12=10,b22=1,b32=1), weight=varwt,data=salmonella) summary(salm.mod4) deviance(salm.mod4) salm.mod5 <- nls(resp ~ expt01*(b02+exp(b1+b22*log(dose)))*exp(-b3*dose) + (1-expt01)*(b01+exp(b1+b21*log(dose)))*exp(-b3*dose), start=c(b01=20,b1=10,b21=1,b3=1,b02=20,b22=1), weight=varwt,data=salmonella) summary(salm.mod5) deviance(salm.mod5) salm.mod6 <- nls(resp ~ expt01*(b0+exp(b1+b22*log(dose)))*exp(-b3*dose) + (1-expt01)*(b0+exp(b1+b21*log(dose)))*exp(-b3*dose), start=c(b0=20,b1=10,b21=1,b3=1,b22=1), weight=varwt,data=salmonella) summary(salm.mod6) deviance(salm.mod6) salm.mod7 <- nls(resp ~ expt01*(b02+exp(b1+b2*log(dose)))*exp(-b3*dose) + (1-expt01)*(b01+exp(b1+b2*log(dose)))*exp(-b3*dose), start=c(b01=20,b1=10,b2=1,b3=1,b02=20), weight=varwt,data=salmonella) summary(salm.mod7) deviance(salm.mod7) salm.mod8 <- nls(resp ~ expt01*(b0+exp(b12+b2*log(dose)))*exp(-b3*dose) + (1-expt01)*(b0+exp(b11+b2*log(dose)))*exp(-b3*dose), start=c(b0=20,b11=10,b2=1,b3=1,b12=10), weight=varwt,data=salmonella) summary(salm.mod8) deviance(salm.mod8) salm.mod9 <- nls(resp ~ expt01*(b02+exp(b1+b2*log(dose)))*exp(-b32*dose) + (1-expt01)*(b01+exp(b1+b2*log(dose)))*exp(-b31*dose), start=c(b01=20,b1=10,b2=1,b31=1,b02=20,b32=1), weight=varwt,data=salmonella) summary(salm.mod9) deviance(salm.mod9) anova(salm.mod3,salm.mod4) anova(salm.mod5,salm.mod4) anova(salm.mod6,salm.mod4) anova(salm.mod7,salm.mod4) anova(salm.mod8,salm.mod4) anova(salm.mod9,salm.mod4) plot(dose,resp,pch=expt) dosev <- seq(.001,3.200,.001) expt1v <- rep(0,length(seq(.001,3.200,.001))) expt2v <- rep(1,length(seq(.001,3.200,.001))) yhatexp1 <- predict(salm.mod6,list(dose=dosev, expt01=expt1v)) yhatexp2 <- predict(salm.mod6,list(dose=dosev, expt01=expt2v)) lines(dosev,yhatexp1,lty=1) lines(dosev,yhatexp2,lty=2) legend("topleft",c("Experiment1","Experiment2"),pch=c(1,2),lty=c(1,2)) dev.off()