require(plyr) simfun <- function(n,beta=c(1,1),tau=0.001, Cor=matrix(c(1,0.25,0.25,1),nrow=2), debug=FALSE) { D <- expand.grid(id=factor(1:n),t=1:3) D$x <- rnorm(3*n) # x_{it} ~ N(0,1) X <- model.matrix(~x,data=D) Z <- model.matrix(~id-1+id:x,data=D) b <- MASS::mvrnorm(n,mu=c(0,0),Sigma=tau^2*Cor) D$y <- rbinom(nrow(D), prob=plogis(X %*% beta + Z %*% c(b)), size=1) D <- D[order(D$id),] ## ADMB expects data in this format (!!) if (debug) save("D",file="current_data.RData") D } cmsg <- c("singular convergence (7)", "false convergence (8)") nMCMC <- 1000 tauvec <- c(0.001,0.5,2) nvec <- c(50,100,500) sfun2 <- function(fun,nMCMC=1000,progress=!debug, debug=FALSE, debugrpt=10, savename="checkpoint",...) { chkfile <- paste(savename,"RData",sep=".") dumpdat <- function(fn) { save("i","j","k","fun","nvec","tauvec", ".Random.seed","res","warn", file=fn) } res <- array(dim=c(2,4,nMCMC,length(tauvec),length(nvec))) dimnames(res) <- list(par=c("b1","b2"), val=c("est","stderr","Z","P"), rep=seq(nMCMC), tau=tauvec, n=nvec) NA_ans <- matrix(NA,nrow=2,ncol=4) warn <- array("",dim=c(nMCMC,length(tauvec),length(nvec))) dimnames(warn) <- dimnames(res)[3:5] for (i in seq_along(tauvec)) { for (j in seq_along(nvec)) { cat("*",i,nvec[i],j,tauvec[j],"\n") ## could more compactly call raply(), but it becomes awkward to save ## the warning info if (progress) pb <- txtProgressBar(style=3) for (k in seq(nMCMC)) { if (debug) { if (k %% debugrpt==0) cat(k,"\n") ## for restoration to current state -- ## assumes only a single run is being debugged at a time ## (could use tmpfile() ...) dumpdat(paste(savename,"current.RData",sep="_")) } if (progress) setTxtProgressBar(pb,k/nMCMC) res[,,k,i,j] <- withCallingHandlers(tryCatch(fun(n=nvec[j],tau=tauvec[i],...), error = function(e) { if (debug) { cat("** caught error\n") ## fixme: try to dump current vals?? save("i","j","k","fun","nvec","tauvec", ".Random.seed","res","warn", file=paste(savename,"current.RData",sep="_")) ## save.image() or save()? dumpdat(paste(savename,i,j,k,"RData",sep=".")) } warn[k,i,j] <<- paste("ERROR:",e$message) NA_ans}), warning = function(w) { if (debug) { cat("** caught warning\n") ## fixme: try to dump current vals?? ## save.image() or save()? dumpdat(paste(savename,i,j,k,"RData",sep=".")) } warn[k,i,j] <<- w$message invokeRestart("muffleWarning") }) } ## for (k in seq(nMCMC)) if (progress) close(pb) ## assumes a single run at a time (use tempfile()?) save("res","warn",file=paste(savename,"RData",sep=".")) } ## for j } ## for i unlink(paste(savename,".RData",sep="")) list(res=res,warnings=warn) } ## can make this a little more simpler/general if we just save results ## in a list format -- then we have to hack them into appropriate shape later ... ## structure as list; get info from first run sfun3 <- function(fun,nMCMC=1000,progress=!debug, debug=FALSE, debugrpt=10, savename="checkpoint",...) { chkfile <- paste(savename,"RData",sep=".") res <- vector("list",nMCMC) warn <- vector("list",nMCMC) dframe <- expand.grid(rep=seq(nMCMC),tau=tauvec,n=nvec) ## try to run once to get appropriate structure ... ? ## ans0 <- fun(n=nvec[1],tau=tauvec[1]) ctr <- 0 for (i in seq_along(tauvec)) { for (j in seq_along(nvec)) { cat("*",i,nvec[i],j,tauvec[j],"\n") if (progress) pb <- txtProgressBar(style=3) for (k in seq(nMCMC)) { if (debug) { if (k %% debugrpt==0) cat(k,"\n") save.image(file=paste(savename,"current.RData",sep="_")) } if (progress) setTxtProgressBar(pb,k/nMCMC) ctr <- ctr+1 res[[ctr]] <- withCallingHandlers(tryCatch(fun(n=nvec[j],tau=tauvec[i],...), error = function(e) { if (debug) { cat("** caught error\n") save.image(file=paste(savename,i,j,k,"RData",sep=".")) } warn[[ctr]] <<- paste("ERROR:",e$message) NA}), warning = function(w) { if (debug) { cat("** caught warning\n") save.image(file=paste(savename,i,j,k,"RData",sep=".")) } warn[[ctr]] <<- w$message invokeRestart("muffleWarning") }) } ## for (k in seq(nMCMC)) if (progress) close(pb) ## assumes a single run at a time (use tempfile()?) save("res","warn","dframe",file=paste(savename,"RData",sep=".")) } ## for j } ## for i unlink(paste(savename,".RData",sep="")) list(res=res,warnings=warn,dframe=dframe) } ## quicker/dirtier -- fails often for smaller n admb.L_fun <- function(...,n,clean=TRUE,verbose=FALSE,is=0, tpl="analyzer_df") { d <- simfun(n=n,...) write_dat(tpl,list(n=n,x=d$x,Y=d$y)) unlink(paste(tpl,".par",sep="")) xargs <- "-mno 10000 -ams 10000000 -noinit" if (is>0) { xargs <- paste(xargs,"-is",is) } run_admb(tpl,extra.args=xargs) g0 <- read_admb(tpl) if (clean) clean_admb(tpl) coef(summary(g0))[4:5,] } ## The idea is that first you fit the model and then check for very small variances ## tau1 and tau2. If one is small you make it inactive and then calculate the ## Hessian for the rest. ## Easiest is to run the model twice. Before running it the first time ## get rid of analyzer.par admb.LX_fun <- function(...,n,clean=TRUE,verbose=FALSE,is=0, tpl="analyzer") { xx <- "-mno 100000 -gbs 200000000 -cbs 100000000" datfile <- paste(tpl,".dat",sep="") parfile <- paste(tpl,".par",sep="") xargs1 <- paste("-ind",datfile, xx,"-nohess") xargs2 <- paste("-ind",datfile, xx,"-phase 10 -ainp",parfile) d <- simfun(n=n,...) write_dat(tpl,list(n=n,x=d$x,Y=d$y)) unlink(parfile) if (is>0) { xargs1 <- paste(xargs1,"-is",is) xargs2 <- paste(xargs1,"-is",is) } run_admb(tpl,extra.args=xargs1,verbose=verbose) run_admb(tpl,extra.args=xargs2,verbose=verbose) g0 <- read_admb(tpl) if (clean) clean_admb(tpl) coef(summary(g0))[4:5,] }