New subject: Predicted values for lm.mixed + patches

6 Jan 2008

Hello!

I noticed that the predicted vaues (PR) in ls.mixed are not caulculated as in 
other {poisson,gamma,probit,logit}.mixed models. In ls.mixed the PR are equal
to 

Y = Xbeta + Zb

while PR in other *.mixed models are equal to

rdist(linkinv(eta))

where 

eta = Xbeta + Zb

Current implementation for PR in ls.mixed is not predicting the values from
the normal distribution, but just adding variability from beta and b posteriors
to fitted values (Xbeta + Zb).

If I got the whole thing right, the code in qi.mixed.R (from line 95)
should be

else if (class(object) == "lmer"){
        ev <- betas %*% t(as.matrix(fTerms))
        mu <- ev
        ## For predicted values, add in random effects draws
        for (i in 1:length(rTerms)){
                mu <- mu + gammas[[names(rTerms[i])]] %*% 
                           t(as.matrix(rTerms[[i]]))
        }
        ## Sample from normal distribution with given mean and variance
        n <- length(mu[,1])
        for (i in 1:ncol(mu)){
                pr[,i] <- rnorm(n, mean=mu[,i], sd=attr(VarCorr(object),
"sc"))
        }
}

Of course ls.mixed vignette should be updated. My proposal is: 

The predicted values (qi$pr) are draws from the gaussian (normal) distribution
defined by mean $\mu_{ij}$ and variance $\sigma^2$,

                              $$\mu_{ij} = X_{ij}\beta + Z_{ij} b_{i}$$

given $X_{ij}$ and $Z_{ij}$ and simulations of $\beta$ and $b_{i}$ from 
their posterior distributions. The estimated variance covariance matrices 
are taken as correct and are themselves not simulated.

Additionally the code in qi.mixed.R uses length(mut[,i]) or length(mut[i,])
within each loop of prediction. I propose the following change starting on
line 65

n <- length(mut[, 1])
if (family == "binomial"){
        ev <- theta
        for (i in 1:ncol(mut)){
                pr[,i] <- as.character(rbinom(n, 1, mut[,i]))
        }
        if (!is.null(y)) {
                if (NCOL(y) > 1) {
                        y <- y[,1]
                }
        }
}
else if (family == "Gamma"){
        ev <- theta * 1/alpha
        n <- length(mut[i,])
        for (i in 1:nrow(mut)){
                pr[i,] <- rgamma(n, shape = mut[i,], scale= 1/alpha)
        }
}
else if (family == "poisson"){
        ev <- theta
        for (i in 1:ncol(mut)){
                pr[,i] <- rpois(n, lambda = mut[,i])
        }
}

Regards, Gregor

-
Zelig Mailing List, served by Harvard-MIT Data Center
Send messages: zelig(a)lists.gking.harvard.edu
[un]subscribe Options: http://lists.gking.harvard.edu/?info=zelig
Zelig program information: http://gking.harvard.edu/zelig/