How do I calculate the likelihood of out-of-sample values for a mixed effects model? I'm trying to use this method for calculating the Information Coefficient using bootstrapping. The advantage of using bootstrapping is that I can compare models that are not nested. But to do this, I need to be able to calculating the likelihood of out-of-sample data (because I'm bootstrapping). 
I have tried several different methods, which give me wildly different results. This is easiest to illustrate when calculating the log-likelihood for the in-sample data. The easiest option is to use logLik:
    data(Orthodont,package="MEMSS")
    mod<-lmer(distance~age+(1+age|Subject), data=Orthodont)
    logLik(mod)

    > -221.3183. 

But I get a different result using the residuals: 
    resid<-residuals(mod)
    sum(dnorm(resid,sd=sd(resid),log=TRUE))

    > -162.1903

I also tried using the residual variance given by lmer:
    sum(dnorm(resid,sd=sigma(mod),log=TRUE))

    > -165.5434

I know that log-likelihood is sometimes calculated by integrating over values for the parameters, whereas by using residuals, I am conditioning on the point-estimates for the parameters. However, according to the help for logLik.merMod, logLik returns "log-likelihood at the fitted value of the parameters." I think that means they are conditioning on the point-estimates. 
Just to be sure, I tried estimating the unconditioned log-likelihood. By using predict with re.form=NA, you can retrieve the fitted values based on fixed effects only (ignoring random effects).
    resid<-Orthodont$distance-predict(mod,newdata=Orthodont,re.form=NA)
    sum(dnorm(resid,sd=sd(resid),log=TRUE))

    > -252.7908

Interestingly, all of the above methods give roughly the same answer when using glm. So this seems to be specific to mixed effects models. 
 A: It seems that calculating log-likelihood for mixed effects models requires dealing with the covariance of error terms for the random effects. Here is a method for calculating log-likelihood by hand for both ML and REML:
data(Orthodont, package="MEMSS")

y <- Orthodont$distance
n <- nrow(Orthodont)    

mod <- lmer(distance ~ age + (1+age|Subject), data=Orthodont, REML=FALSE)
logLik(mod)

G <- diag(attr(VarCorr(mod)$Subject, "stddev")) %*% attr(VarCorr(mod)$Subject, "correlation") %*% diag(attr(VarCorr(mod)$Subject, "stddev"))
V <- lapply(split(Orthodont, Orthodont$Subject), function(x) cbind(1, x$age) %*% G %*% rbind(1, x$age) + diag(rep(sigma(mod)^2, nrow(x))))
V <- as.matrix(bdiag(V))
W <- solve(V)
X <- cbind(1, Orthodont$age)
b <- fixef(mod)

dmvnorm(y, mean = X %*% b, sigma=V, log=TRUE)
c(-n/2 * log(2*pi) - 1/2 * log(det(V)) - 1/2 * t(y - X %*% b) %*% W %*% (y - X %*% b))

mod <- lmer(distance~age+(1+age|Subject), data=Orthodont, REML=TRUE)
logLik(mod)

G <- diag(attr(VarCorr(mod)$Subject, "stddev")) %*% attr(VarCorr(mod)$Subject, "correlation") %*% diag(attr(VarCorr(mod)$Subject, "stddev"))
V <- lapply(split(Orthodont, Orthodont$Subject), function(x) cbind(1, x$age) %*% G %*% rbind(1, x$age) + diag(rep(sigma(mod)^2, nrow(x))))
V <- as.matrix(bdiag(V))
W <- solve(V)
X <- cbind(1, Orthodont$age)
b <- fixef(mod)
p <- length(b)

c(-(n-p)/2 * log(2*pi) - 1/2 * log(det(V)) - 1/2 * log(det(t(X) %*% W %*% X)) - 1/2 * t(y - X %*% b) %*% W %*% (y - X %*% b))

To calculated the likelihood for a new datapoint (or, more accurately, calculate the density for that new datapoint) using ML, given X and Y for the new subjects, compute V based on G and then calculate:
dmvnorm(y, mean = X %*% b, sigma=V, log=TRUE)
c(-n/2 * log(2*pi) - 1/2 * log(det(V)) - 1/2 * t(y - X %*% b) %*% W %*% (y - X %*% b))

It's not clear this can be done using REML, since the - 1/2 * log(det(t(X) %% W %% X)) term cannot be decomposed into the contribution of each individual subject.
Many thanks to several experts who answered questions via email.
A: For fitted lme4 models, logLik() conditions on the point estimate of (i) the fixed-effects parameters, and (ii) the random-effects covariance matrix.  It does not condition on point estimates of the random effects coefficients (e.g., subject-specific random intercepts or slopes); it integrates them out.  This is also how the log-likelihood presented by summary(mod) is computed.
Your first alternative method, resid(mod), conditions not only on (i) and (ii) above but also on point estimates of the random effects coefficients.  Your second alternative method, using predict() with the flag re.form=NA, discards the random effects coefficients altogether.  So all three of these methods will give different results.
Since vanilla generalized linear models don't have random effects, none of these issues arise and so you don't get differences across these methods with glm().
Regarding how to compute out-of-sample (and by this I think you mean observations for which the random effects factor level is not in the data used to fit the model) -- the technically correct thing to do is to integrate out the random effects for those observations too.  I'm not aware that lme4 hands you methods for doing this, but you could use the math in the documentation (see https://cran.r-project.org/web/packages/lme4/vignettes/lmer.pdf, section 3).
Finally: I'm not familiar with the method in the paper that you link to, but are you sure that it is applicable to datasets with hierarchical structure?
(Side note: it looks the documentation for logLik.merMod() may be a bit out of date as mod@resp$aic() returns an error.)
A: I can't provide a definitive answer here, but I've taken a look through the lme4 source code and it seems like, as with many areas of this package, there is an active discussion about how best to calculate the deviance of the model which, in turn, is used by logLik() to calculate the log-likelihood.
The key functions here are logLik(), devCrit(), and deviance().
You are also likely to be interested in this writeup on deviance and log-likelihood in lme4.
As for the problem at hand, I would use your last method (using residuals of unconditional estimates) for out-of-bag log-likelihood estimates, while noting that the distribution of residuals for your model fit to Orthodont data are not really normally-distributed.
