[R-sig-ME] lme with heteroscedacity- how to plot within-person variability

[Ben Bolker]

Rapsomaniki, Eleni <e.rapsomaniki at ...> writes:

> 
> Dear list members,
> 
> I wasn't sure whether to post this in the sig-mixed-models list, as it
concerns the nlme package.
> 
> Suppose the following model:
> model.lme=lme(distance ~ Sex*age, data=Orthodont,random=~age|Subject,
> weights=varIdent(form=~1|Sex),correlation=corExp(form=~age, nugget=T))
 
> I would like to plot within-subject variability (e.g. the standard
> deviation of the fitted 'distance') over age. I expect that the
> vector of this variability (one per person per observed age?) is
> returned in model.lme, but I'm a bit confused with what to extract
> from the output.

  There may be a good better to do this, but here's what I came up
with:


library(nlme)
model.lme=lme(distance ~ Sex*age, data=Orthodont,random=~age|Subject,
 weights=varIdent(form=~1|Sex),
correlation=corExp(form=~age, nugget=TRUE))


## extract relative sd terms per sex
sdw <- 1/unique(varWeights(model.lme$modelStruct))
## fixed-effect variance-covariance matrix
fv <- vcov(model.lme)
## data for prediction
dd <- expand.grid(Sex=factor(levels(Orthodont$Sex)),
      age=seq(min(Orthodont$age),max(Orthodont$age),length.out=41))
## design matrix
X <- model.matrix(~Sex*age,data=dd)
## sample variances of predictions of a linear model
dd$cvar <- diag(X %*% fv %*% t(X))
## residual variances, scaled
##    by sex-specific weighting (make sure the order is right!)
dd$pvar <- summary(model.lme)$sigma^2*sdw[as.numeric(dd$Sex)]^2
dd$sd <- sqrt(dd$cvar+dd$pvar)

## one way to plot the results
library(lattice)
xyplot(sd~age,group=Sex,data=dd,type="l",
   ylim=c(0,2))

Note that these are prediction intervals (including the
residual variance, which is what the heteroscedasticity
term applies to). The estimate of the model standard deviations
(i.e. based on the fixed effect sampling variances) does *not*
take the uncertainty of the variances themselves into account.
 
> Also regarding the nugget, should I initialise it to some value? Or
> just assume the value chosen by the model fitted are the actual
> 'measurement error' (please correct me if this is not what the
> nugget represents)?

I would let the nugget be estimated automatically unless you
have reason to believe the estimate is wonky.
 
> Finally, am I right to think that heteroscedastic variances are
>  still not possible with lme4?

  Correct.