load("../DATA/Owls.rda")
library(R2jags)
load.module("glm")
source("owls_BUGS_funs.R")
## ooo <- owls_BUGS_fit(Owls,ni=1000,nb=200,nt=10,nc=3)
## unique(names(list.samplers(ooo$model)))
st1 <- system.time(ofit <- owls_BUGS_fit(Owls))
st1 <- unname(st1["elapsed"])

dropdev <- function(x) {
  ## FIXME: doesn't change last.values ?
  if (inherits(x,"bugs")) {
    w <- which(x$root.short=="deviance")
    for (i in grep("\\.short$",names(x))) {
      x[[i]] <- x[[i]][-w]
    }
    x <- lapply(x,dropdev)
    class(x) <- "bugs"
    return(x)
  }
  if (!is.null(dim(x))) {
    dn <- dimnames(x)
    for (i in seq_along(dn)) {
      if ("deviance" %in% dn[[i]]) {
        w <- which(dn[[i]]=="deviance")
        x2 <- x[slice.index(x,i)!=w]
        dx <- dim(x)
        dx[i] <- dx[i]-1
        dn[[i]] <- dn[[i]][-w]
        return(array(x2,dim=dx,dimnames=dn))
      }
    }
  } else if ("deviance" %in% names(x)) {
    return(x[names(x)!="deviance"])
  }
  x
}

save("ofit",file="owls_BUGS_fits.RData")

library(coda)
library(emdbook)
ofit_mcmc <- as.mcmc(ofit$BUGSoutput)
c1 <- dropdev(unlist(ofit$BUGSoutput$mean))
c1["psi"] <- -qlogis(c1["psi"])
c2 <- dropdev(HPDinterval(lump.mcmc.list(ofit_mcmc)))
c1["sigma"] <- c1["sigma"]^2
c2["psi",] <- -qlogis(c2["psi",])
c2["sigma",] <- c2["sigma",]^2
results <- list(obj=NULL,
                coef=c1,
                confint.mcmc=c2,
                time=c(fit=unname(st1),mcmc=unname(st1)),
                terminfo=NULL,
                convinfo=list(effsize=dropdev(effectiveSize(ofit_mcmc)),
                  gelman=dropdev(gelman.diag(ofit_mcmc)$psrf[,"Point est."])))

save("results",file="fit.RData")





### center (and scale?) covariates
## orthogonalize design matrix?
## load glm module?
## design matrix vs