[R-sig-ME] SE for predictions in lme

[Philip Harrison]

Hi Marta,

My understanding of this is that the SE*2 gives 95% confidence  
intervals on predictions. You can add a bonferroni type adjustment by  
varying the *2 (which is really 1.96)

The SE2 gives prediction intervals. Confidence intervals of  
predictions tell you how well you have predicted the mean. Prediction  
intervals tell where you could expect the next to see the next value  
sampled. The key point is that the prediction interval tells you about  
the distribution of values, not the uncertainty in determining the  
population mean.

Therefore as far as I can tell, the CIs are the most useful for making  
inferences about differences in the conditional means/BLUPs

You can also use the function predictSE.lme from the package  
AICcmodavg to produce prediction and SEs from an lme object and then  
calculate CIs accordingly. I like this method as it gives a reference  
to the method:

"predictSE.lme? computes predicted values based on fixed effects and  
associated standard errors.
Standard errors are approximated using the delta method (Oehlert 1992)"

These two methods result in identical predictions and almost identical  
SEs for my dataset

Some confusion arise for me because the predict.lme function calls the  
prediction Best linear Unbiased Predictions (BLUPs) which combine  
random and fixed effects, and the predictSE.lme function states that  
these predictions do not incorporate the random effects. Yet they give  
identical predictions for me.

Hope this helps

Phil


Quoting marta rufino <marta.m.rufino at gmail.com>:

> Hi,
>
> This topic has been discussed in several posts in this list before, but I
> could not completely understand if there was a clear answer. Also, I
> searched the books (Pinheiro, Faraway and Zuur) I could not reach a
> conclusion. I guess I am missing some point here :) so I am sorry for
> recover this topic again.
> I will start by providing an example code and then put my questions, as
> clear as possible.
>
> Following the code provided in previous posts (see for example):
> http://markmail.org/message/lhtols3t5wrleewc
> using another dataset, for example:
>
>            # Example to ask the list
>            data(iris)
>            kk=data.frame(stack(iris[,1:3]),sp=rep(iris$Species),
> subject=rep(1:dim(iris)[1],3))
>            summary(kk)
>            require(nlme)
>            kk2=lme(values~ind*sp, kk, weights = varIdent(form= ~ 1|ind),
> random=~1|subject)
>            kk2
>            anova(kk2) #all factors are significant.
>
>            # Calculate model predictions    and plot it ==== this is from
> Ben's code
>                newdat <- expand.grid(ind=levels(kk$ind), sp=levels(kk$sp))
>                newdat$pred <- predict(kk2, newdat, level = 0)
>
>                Designmat <- model.matrix(eval(eval(kk2$call$fixed)[-2]),
> newdat[-3])
>                predvar <- diag(Designmat %*% kk2$varFix %*% t(Designmat))
>                newdat$SE <- sqrt(predvar)
>                newdat$SE2 <- sqrt(predvar+kk2$sigma^2)
>
>                #install.packages("ggplot2")
>                library(ggplot2)
>                pd <- position_dodge(width=0.4)
>                ggplot(newdat,aes(x=ind, y=pred, colour=sp))+
>                geom_point(aes(x=kk$ind,y=kk$values, colour=kk$sp),
> size=.3, shape=2, position=pd)+
>                geom_point(position=pd)+
>                geom_linerange(aes(ymin=pred-2*SE2,
> ymax=pred+2*SE2),col="red", position=pd)+
>                geom_errorbar(aes(ymin=pred-2*SE, ymax=pred+2*SE),
> col="black", width=.1, position=pd)
>
> So, we can see the actual points, the SE and SE2 modelled and respective
> means, I think.
> My questions are:
> What is exactly the SE and SE2 (how do we call it in an article legend, for
> example) and how can these be interpreted?
> Can we consider that when the 'lines' do not overlap the species (or
> varieties) are different?
> The SE2 are larger than the actual data in setosa/Petal.Length and overlap
> in the remaining variaties--- how can these show significant differences
> than?
>
> I think I am getting confused or I am missing something important, because
> the results do not appear congruent in my head.
>
> Any help will be much appreciated,
>
> All the best,
> Marta
>
> 	[[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>
>



Philip Harrison MSc
PhD student (Fisheries Ecology)
Department of Biology
University of Waterloo
200 University Avenue West
Waterloo, Ontario, Canada
N2L 3G1
Cell:226-808-2309
Email:pharriso at uwaterloo.ca