rm(list=ls(all=TRUE)) # clear off .GlobalEnv graphics.off() # clear off old graphics windows options(stringsAsFactors = FALSE) wt=structure(list(day = c(2, 4, 6, 9, 13, 15), K2p5 = c(54, 30.75, 26.875, 23.375, 45.625, 90.5), K5 = c(57.875, 34.75, 29.125, 21.625, 40.875, 71.25), K10 = c(48.25, 29.125, 26.375, 15.875, 31.875, 69)), .Names = c("day", "K2p5", "K5", "K10"), row.names = 2:7, class = "data.frame") matplot(wt[,1],wt[,2:4],type="b", xlab="Days",ylab="10e4 cells") legend(7,80,bty="n",pch=as.character(1:3),c("2.5 uM","5 uM","10 uM")) wt library(odesolve) g=NULL g$d1=wt[,c(1,2)] g$d2=wt[,c(1,3)] g$d3=wt[,c(1,4)] names(g$d1)[2]="N" names(g$d2)[2]="N" names(g$d3)[2]="N" g$nData=dim(g$d1)[1]+dim(g$d2)[1]+dim(g$d3)[1] # -3 # minus 3 if IC fixed to data g delta0=1 parICs=c(beta=.8, delta1=delta0,delta2=delta0,delta3=delta0,mu=5e-5,S1=wt[1,2],S2=wt[1,3],S3=wt[1,4]) # #parICs=c(beta=0.9,mu=5e-6, delta1=0.8,delta2=1,delta3=1.2,S1=wt[1,2],S2=wt[1,3],S3=wt[1,4]) g$params=data.frame(initial=parICs,final=parICs,opt=TRUE,CI95prct="not fitted") # this chunk shows that the DLL option is ~100x faster DLL=FALSE DLL=TRUE source("C:/Users/radivot/src/htdocs/EPBI473/files/wk13tumorTherapy/baCommon.R") if (DLL) dyn.load("C:/Users/radivot/src/htdocs/EPBI473/files/wk13tumorTherapy/baMod.dll") t0 = Sys.time(); g=simData(g) difftime(Sys.time(),t0,units="secs")[[1]] plotN(g) g$d1 g$params[c("S1","S2","S3"),"opt"]=FALSE # comment to fit all g$params p0=t(g$params[g$params[,"opt"],"initial",drop=F]) p0=unlist(as.data.frame(p0)) p0=sapply(p0,log) g$optNames<-names(p0) g g$params t0 = Sys.time(); opt<-optim(p0,fopt,hessian=TRUE,model=g,control=list(trace=F,maxit=5000)) difftime(Sys.time(),t0,units="mins")[[1]] opar<-opt$par Hess=opt$hessian value=opt$value sg<-sqrt(value/(g$nData-length(p0))) if (abs(det(Hess))>1e-8) { sig=sg*sqrt(diag(solve(Hess/2))) upper=signif(opar+1.96*sig,3) lower=signif(opar-1.96*sig,3) g$CI=exp(cbind(lower,upper)) } g$params[g$optNames,"final"]=signif(exp(opar),3) g$params[g$optNames,"CI95prct"]=paste("(",signif(g$CI[,"lower"],3),", ",signif(g$CI[,"upper"],3),")",sep="") g$params # now look at final parameter estimate situation g=simData(g) plotN(g) if (DLL) dyn.unload("C:/Users/radivot/src/htdocs/EPBI473/files/wk13tumorTherapy/baMod.dll") #install.packages("hwriter") library(hwriter) hwrite(g$params, 'C:/Users/radivot/src/htdocs/EPBI473/files/wk13tumorTherapy/baMod.html')