fitmodel <- function(bs,d) { frm <- "DV ~ " first <- 1 for (i in 1:m) { if (bs[i]==1) { if (first) { frm <- paste(frm,"E",i,sep="") first <- 0 } else { frm <- paste(frm,"+E",i,sep="") } } } frm <- paste(frm,"-1",sep="") frm <- as.formula(frm) lm(frm,data=d,weights=1/WT^2) } AICs <- function(fm,d) { pd <- predict.lm(fm) mssq <- mean(((d$DV-pd)/d$WT)^2) m2ll <- log(mssq) + log(2*pi) + 1 k <- length(fm$coefficients)+1 aic <- m2ll + 2*k/mn if (mn-k-1 >= 1) { aicct <- k*(1+(k+1)/(mn-k-1)) aicc <- m2ll + 2*aicct/mn } else { aicc <- NA } mpe2 <- mean(((d$DVV-pd)/d$WT)^2) aicv <- log(mssq) + log(2*pi) + mpe2/mssq c(m2ll, aic, aicc, aicv, mpe2, mssq) } run <- function(i,d,bs) { b <- bs[[i]] fm <- fitmodel(b,d) afm <- AICs(fm,d) c(0,afm[1],afm[2],afm[3],afm[4],afm[5],afm[6],0) } # the subdirectories should have been created with ``mwd.sh'' runnm <- function(irun,bs,nd,ndv) { b <- bs[[irun]] setwd(paste(wd,"/",irun,sep="")) cat("; CONSTRUCT YN\n",file="aicyn.inc") yn <- "YNE=" np <- 0 for (j in 1:m) { if (b[j]==1) { np <- np+1 cat(paste("LAEP",j,"=THETA(",np,")*E",j,"\n",sep=""),file="aicyn.inc",append=T) if (np==1) { yn <- paste(yn,"LAEP",j,sep="") } else { yn <- paste(yn,"+LAEP",j,sep="") } } } cat(paste(yn,"\n",sep=""),file="aicyn.inc",append=T) if (np==1) { cat(paste("(0.01)\n",sep=""),file="aictheta.inc") } else { cat(paste("(0.01)x",np,"\n",sep=""),file="aictheta.inc") } if (what) { system("rm -f aicwh.xml") system("n73dx aicwh.nmt",ignore.stdout=T) x <- try(xmlTreeParse("aicwh.xml", getDTD=F), silent=T) } else if (wobs) { system("rm -f aicwo.xml") system("n73dx aicwo.nmt",ignore.stdout=T) x <- try(xmlTreeParse("aicwo.xml", getDTD=F), silent=T) } else if (naive) { system("rm -f aicn.xml") system("n73dx aicn.nmt",ignore.stdout=T) x <- try(xmlTreeParse("aicn.xml", getDTD=F), silent=T) } else { system("rm -f aic.xml") system("n73dx aic.nmt",ignore.stdout=T) x <- try(xmlTreeParse("aic.xml", getDTD=F), silent=T) } system("mv report aic.rep") iere <- -1 if (class(x)[[1]]!="try-error") { r <- xmlRoot(x) iere <- as.double(xmlValue(r[['nonmem']][['problem']][['estimation']][['termination_status']])) mvof <- as.double(xmlValue(r[['nonmem']][['problem']][['estimation']][['final_objective_function']])) if (is.na(mvof)) iere <- -2 } if (iere>=0) { ethetas <- as.double(sapply(1:np, function(i) xmlValue(r[['nonmem']][['problem']][['estimation']][['theta']][[i]]))) eomega2 <- as.double(xmlValue(r[['nonmem']][['problem']][['estimation']][['omega']][[1]])) esigma2 <- as.double(xmlValue(r[['nonmem']][['problem']][['estimation']][['sigma']][[1]])) if (w2>0 && eomega2<0.0001 && !naive) iere <- -3 # the predictions calculated next may be used as weights ypr <<- rep(0,mn) jr <- 1 for (j in 1:m) { if (b[j]==1) { ypr <<- ypr + ndv[[paste("E",j,sep="")]]*ethetas[[jr]] jr <- jr+1 } } for (j in 1:np) { if (j==1) { cat(paste(ethetas[1]," FIX\n",sep=""),file="aicthetav.inc") } else { cat(paste(ethetas[j]," FIX\n",sep=""),file="aicthetav.inc",append=T) } } cat(paste(eomega2," FIX\n",sep=""),file="aicomega2v.inc") cat(paste(esigma2," FIX\n",sep=""),file="aicsigma2v.inc") if (what) { system("n73dx aicwhv.nmt",ignore.stdout=T) x <- xmlTreeParse("aicwhv.xml", getDTD=F) } else if (wobs) { system("n73dx aicwov.nmt",ignore.stdout=T) x <- xmlTreeParse("aicwov.xml", getDTD=F) } else if (naive) { system("n73dx aicnv.nmt",ignore.stdout=T) x <- xmlTreeParse("aicnv.xml", getDTD=F) } else { system("n73dx aicv.nmt",ignore.stdout=T) x <- xmlTreeParse("aicv.xml", getDTD=F) } r <- xmlRoot(x) mvofv <- as.double(xmlValue(r[['nonmem']][['problem']][['estimation']][['final_objective_function']])) # calculate mean prediction error if (what) { idx <- which(ypr!=0) if (length(idx)= 1) { aicc <- mvof + 2*k*(1+(k+1)/(mn-k-1))/mn } else { aicc <- NA } aicv <- mvofv # when w^2 this is equal to; but let NONMEM calculate the marginal likelihood # aicv <- mlw + log(esigma2) + log(2*pi) + mpe2/esigma2 } else { mvof <- 0 aic <- 0 aicc <- 0 aicv <- 0 mpe2 <- 0 eomega2 <- 0 esigma2 <- 0 } setwd(wd) c(iere,mvof,aic,aicc,aicv,mpe2,esigma2,eomega2) } runm <- function(i,d,bs) { print(i) d$DV <- ynp*(rep(rlnorm(n,0,w),each=m)+rnorm(n*m,0,s)) d$DVV <- ynp*(rep(rlnorm(n,0,w),each=m)+rnorm(n*m,0,s)) lapply(1:maxms, run, d, bs) } runmnm <- function(i,d,bs) { print(i) if (what) { if (tnrm) { ynpi <- ynp*rep(rlnorm(n,0,w),each=m) yp <- ynpi*(1+rtruncnorm(n*m,-1,1,0,s)) ynpiv <- ynp*rep(rlnorm(n,0,w),each=m) ypv <- ynpiv*(1+rtruncnorm(n*m,-1,1,0,s)) } else { ynpi <- ynp*rep(rlnorm(n,0,w),each=m) yp <- ynpi*(1+rnorm(n*m,0,s)) ynpiv <- ynp*rep(rlnorm(n,0,w),each=m) ypv <- ynpiv*(1+rnorm(n*m,0,s)) } nd$DV <- yp ndv$DV <- ypv } else { if (tnrm) { nd$DV <- ynp*(rep(rlnorm(n,0,w),each=m)+rtruncnorm(n*m,-1,1,0,s)) ndv$DV <- ynp*(rep(rlnorm(n,0,w),each=m)+rtruncnorm(n*m,-1,1,0,s)) } else { nd$DV <- ynp*(rep(rlnorm(n,0,w),each=m)+rnorm(n*m,0,s)) ndv$DV <- ynp*(rep(rlnorm(n,0,w),each=m)+rnorm(n*m,0,s)) } } write.table(nd,"aic.d",quote=F,row.names=F) write.table(ndv,"aicv.d",quote=F,row.names=F) if (whats) { what <- T rrm <- list(runnm(maxms, bs, nd, ndv)) nd$WT <- ypr ndv$WT <- ypr write.table(nd,"aic.d",quote=F,row.names=F) write.table(ndv,"aicv.d",quote=F,row.names=F) what <- F rr <- mclapply(1:maxms1, runnm, bs, nd, ndv) rr[maxms] <- rrm rr } else { mclapply(1:maxms, runnm, bs, nd, ndv) } }