
## sh*sc = mean
## sh*sc^2 = var
## sc = var/mean = phi
## sh = m^2/var = l^2/(l*phi) = l/phi
## horrible hack to generate "quasipoisson" data
rqpois <- function(n,lambda,phi,roundvals=FALSE) {
  sc <- phi
  sh <- lambda/phi
  r <- rgamma(n,shape=sh,scale=sc)
  if (roundvals) round(r) else r
}

## or via negative binomial with mean-dependent k:
## mean = mu
## var = mu*(1+mu/k)
## var = phi*mu
## (1+mu/k) = phi*mu
## k= mu/(phi*mu-1)
## provided (phi*mu>1) we're OK ...
rqpois2 <- function(n,lambda,phi,min.k=1e-3,roundvals=FALSE) {
  mu <- lambda
  k <- mu/(phi*mu-1)
  if (any(k<min.k)) warning("some overdispersion parameters rounded up to",min.k)
  k <- pmax(min.k,k)
  r <- rnbinom(n,mu=mu,size=k)
}

my.mer.sim <- function (object, nsim = 1, seed = NULL, ...) 
{
    if (!is.null(seed)) 
        set.seed(seed)
    if (!exists(".Random.seed", envir = .GlobalEnv)) 
        runif(1)
    ## FIXME: implement offset, should be fairly easy (?)
    if (!is.null(object@call$offset))
      warning("offset ignored -- not yet implemented")
    RNGstate <- .Random.seed  ## why saved?
    dims <- object@dims
    sigma <- lme4:::sigma(object) ## FIXME: can remove prefix once inside namespace
    etasim.fix <- as.vector(object@X %*% fixef(object))  ## fixed-effect contribution
    etasim.reff <- (as(t(object@A) %*% matrix(rnorm(nsim * dims["q"]), nc = nsim), 
                       "matrix"))  ## UNSCALED random-effects contribution
    if (length(object@V) == 0 && length(object@muEta) == 0)  {
      etasim.resid <- matrix(rnorm(nsim * dims["n"]), nc = nsim) ## UNSCALED residual
      etasim <- etasim.fix + sigma*(etasim.reff+etasim.resid)
      return(etasim)
    }
    if (length(object@muEta)>0) {
      ## GLMM
      ## n.b. DON'T scale random-effects
      ## if sigma!=1, it applies to the "quasi"- part of the model
      etasim <- etasim.fix+etasim.reff 
      family <- object@call$family
      if(is.character(family))
        family <- get(family, mode = "function", envir = parent.frame(2))
      if(is.function(family)) family <- family()
      if(is.null(family$family)) stop("'family' not recognized")
      musim <- family$linkinv(etasim)
      n <- length(musim) ## FIXME: or could be dims["n"]?
      vsim <- matrix(switch(family$family,
                     poisson=rpois(n*nsim,lambda=musim),
                     quasipoisson=rqpois(n*nsim,musim,sigma), ## quasipois hack,
                     binomial={
                       ## FIXME: not general/robust?
                       ## what is the best way to extract
                       size <- rowSums(object@frame[,1:2])
                       rbinom(n*nsim,prob=musim,size=size)
                     },
                     stop("simulation not implemented for family",family$family)),
                     nc=nsim)
      return(vsim)
    }
    stop("simulate method for NLMMs not yet implemented")
}

if (FALSE) {
  ## simulate "data"
  b_eff = rnorm(10,0,sd=0.1)
  a = 1
b = 1

zz = expand.grid(block=1:10,x=1:10,rep=1:10)
zz <- within(zz,
       { eta0 <- a+b*x+b_eff[block];
         eta1 <- rnorm(length(eta0),eta0,0.1)
         y <- rpois(length(eta0),exp(eta0))
         })

zz2 <- within(zz,
       { eta0 <- a+b*x+b_eff[block];
         eta1 <- rnorm(length(eta0),eta0,0.1)
         y <- rbinom(length(eta0),prob=plogis(eta0),size=20)
         })

library(lme4)

m3 <- lmer(cbind(y,20-y)~x+(1|block),data=zz2,family="binomial")
my.mer.sim(m3)
m1 <- lmer(y~x+(1|block),data=zz,family="poisson")
m0 <- lmer(eta1~x+(1|block),data=zz)

s0 <- simulate(m0,seed=1001)
s1 <- my.mer.sim(m0,seed=1001)
all(abs(s0-s1)<1e-10)

s2 <- my.mer.sim(m1)
m2 <- refit(m1,s2)

zzz = replicate(100,refit(m1,my.mer.sim(m1)))
library(reshape)
coefs <- t(sapply(zzz,fixef))

densityplot(~value|X2,data=melt(coefs),
            scales=list(relation="free"))

plot(jitter(s2),jitter(zz$y))
lines(lowess(zz$y~s2))
dotplot(ranef(m1))
lme4ns <- environment(getMethod("simulate","mer"))
setMethod("simulate", "mer", my.mer.sim,
               where = lme4ns)




zz = mlist4[[1]]
slotNames(zz)

      

roundgammavar <- function(shape,scale) {
  maxv <- round(qgamma(0.999,shape=shape,scale=scale))
  vec <- 0:maxv
  v <- c(0,vec+0.5)
  pdfdvals <- pgamma(v[-1],shape=shape,scale=scale)-
    pgamma(v[-length(v)],shape=shape,scale=scale)
  c(sh=shape,sc=scale,
    mx=sum(pdfdvals*vec),vx=sum(pdfdvals*(vec-shape*scale)^2),
    m0=shape*scale,v0=shape*scale^2)
}

## Think about quasipois hack: why does rounding *increase* variance?
## Ask the R list???

dotplot(ranef(xxx))

sf <- function(x) { c(mean(x),var(x))}
rr = replicate(50,sf(round(rgamma(10000,1,2))))
rr2 = replicate(250,sf(rgamma(10000,1,2)))
plot(density(rr2[2,]))
rug(rr2[2,])
plot(density(rr[2,]))
}



