Sorry, there was an error in my code. I have corrected the code below, now it exactly yields the same results as lme4.
There is only the last question remaining (why the constant column).

SL



________________________________

À : "r-sig-mixed-models@r-project.org" <r-sig-mixed-models@r-project.org> 
Envoyé le : Vendredi 29 juin 2012 10h28
Objet : [R-sig-ME] profile likelihood confidence interval

Dear all,

I'm trying to reproduce with my bare hands the confidence interval given by confint(profile(fit)) with the latest development version of lme4 (version 0.99999911-0 with R-2.15.1 on Windows XP 32 bits).

I have done this job for the simplest mixed model: the balanced random one-way ANOVA model. I have found a little discrepancy between my results and the lme4 results. I would like to have an explanation for that. My code is given below. Moreover I ask a question at the end of the code: the profile() function returns a constant column, I do not understand why.

Best regards,
Stéphane


library(lme4)
library(mvtnorm)
library(optimx)


## function for simulating data
sims <- function(I, J, sigmab0, sigmaw0){
      Mu <- rnorm(I, mean=0, sd=sigmab0)
      y <- c(sapply(Mu, function(mu) rnorm(J, mu, sigmaw0)))
data.frame(y=y, group=gl(I,J))    
}


## set parameters 
I <- 3 # number of groups
J <- 8 # number of repeats per group
sigmab0 <- 0.15 # between standard deviation
sigmaw0 <- 0.15 # within standard deviation 


## likelihood function 
LHD <- function(dd, mu, sigma2b, sigma2w){
    ag <- aggregate(y~group, data=dd, FUN=length)
    J <- length(ag$group)
    K <- ag$y
    L <- 1
    for(j in 1:J){
        M1 <- cbind(rep(sigma2b, K[j]), diag(sigma2w, K[j]))
        M2 <- rbind(rep(1, K[j]), diag(K[j]))
        Sigma <- M1 %*%M2
        y <- subset(dd, subset=group==as.character(j))$y
        L <- L*dmvnorm(y, rep(mu, K[j]), Sigma)
    }
L
}

# simulates data
set.seed(666)
dd <- sims(I, J, sigmab0, sigmaw0)


####################################
### confidence interval about mu ###
####################################

## firstly, determine a range for the confidence bounds 
( fit <- lmer(y~(1|group), data=dd, REML=FALSE) )
mu.hat <- getME(fit, "beta")
se <- coef(summary(fit))[1, "Std. Error"]
low <- mu.hat-4*se
up <- mu.hat+4*se
mu <- seq(low, up, le=51)


## determine the profile likelihood for mu 
l <- length(mu)
sigmas.hat <- array(NA, dim=c(l,2))
for(i in 1:l){
    f <- function(sigmas){
        sigma2b <- sigmas[1]^2
        sigma2w <- sigmas[2]^2
        LHD(dd, mu[i], sigma2b, sigma2w)
    }
    sigmas.hat[i,] <- optimx(c(0.4,0.2), f, method="Nelder-Mead", control=list(fnscale=-1))$par$par
}


## plot the profile likelihood 
PL <- sapply(1:l, function(i) 
    LHD(dd, mu[i], sigmas.hat[i,1]^2, sigmas.hat[i,2]^2)  
)
plot(mu,PL, type="l")
abline(v= mu.hat, lty="dotted")
abline(h= exp(as.numeric(logLik(fit))) , lty="dotted")


## determine the profile likelihood confidence interval 
F <- splinefun(mu,log(PL))
curve(F(x), from=low, to=up, xlab="mu (Intercept)", 
    main="profile likelihood")
abline(v= mu.hat, lty="dotted", col="blue")
abline(h= as.numeric(logLik(fit)), lty="dotted", col="green")
abline(h= as.numeric(logLik(fit))- qchisq(95/100,1)/2, lty="dotted", col="yellow") 
FF <- function(x) F(x) - as.numeric(logLik(fit))+ qchisq(95/100,1)/2
(upb <- uniroot(FF, c(mu.hat, up))$root) # upper confidence bound
(lowb <- uniroot(FF, c(low, mu.hat))$root) # lower confidence bound
abline(v=c(lowb,upb), lty=2, col="red")


## compare my confidence interval with the one given through the profile() function 
pr <- profile(fit)
confint(pr)


## compare profile likelihood
subset(pr, subset= .par=="(Intercept)") # why is the column .sigma constant ?
cbind(sigmas.hat, mu)
    [[alternative HTML version deleted]]


_______________________________________________
R-sig-mixed-models@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
	[[alternative HTML version deleted]]

