Re: [R] nlme predict with se?

From: Berton Gunter <gunter.berton_at_gene.com>
Date: Sat 18 Mar 2006 - 10:08:48 EST


This won't be much help, I'm afraid but ...

The problem is, of course, that the prediction is a **nonlinear** function of the parameters (including possibly the variance components, depending on what level of the variance hierarchy you are trying to predict). So you need to somehow estimate how to propogate the error of a nonlinear function of parameters which, themselves, have an approximate variance-covariance matrix estimated from an approximation to the likelihood surface. Asymptotic approaches (Taylor series/delta method; Laplace approximations) can do this, but they are of questionable accuracy, as Doug Bates has repeatedly emphasized on this list.

A better approach might well be to bootstrap: repeatedly fit and predict your desired values via bootstrapping your original data "appropriately". This, too, could be somewhat tricky, as you have to bootstrap from the variance components appropriately (preserving group identities, for example).

As I have no real experience with this, take this with more than just a grain of salt. But it might give you some insight into the difficulty and why you were unable to find "good" answers.

Cheers,

"The business of the statistician is to catalyze the scientific learning process." - George E. P. Box    

> -----Original Message-----
> From: r-help-bounces@stat.math.ethz.ch
> [mailto:r-help-bounces@stat.math.ethz.ch] On Behalf Of Dr.
> Emilio A. Laca
> Sent: Friday, March 17, 2006 1:55 PM
> To: R-help@stat.math.ethz.ch
> Subject: [R] nlme predict with se?
>
> I am trying to make predictions with se's using a nlme (kew11.nlme
> below). I get an error indicating levels for a factor are
> not allowed.
> I have searched and read Rnews, MEMSS, MASS, R-Help, and other lists
> in Spanish where I found questions similar to mine but not solution.
> I do not really care about the method used. Any suggestions
> to obtain
> predictions with se's from an nlme would be appreciated.
> eal
>
>
> R Version 2.2.1 (2005-12-20 r36812)
> mac osx 10.4.5 on G5 DP
>
> > predict(kew11.nlme,x.for.lw.pred, se=T)
> Error in predict.nlme(kew11.nlme, x.for.lw.pred, se = T) :
> Levels 1. Fall-Winter,2. Spring,3. Summer,4. Fall not
> allowed for
> Season
>
> > x.for.lw.pred[1:10,]
> Season ecoreg MBreed ageyr2 sheep farm
> 1 1. Fall-Winter desert Meat 0.60 1 AksToKa
> 2 2. Spring desert Meat 1.00 1 AksToKa
> 3 3. Summer desert Meat 1.25 1 AksToKa
> 4 4. Fall desert Meat 1.50 1 AksToKa
> 5 1. Fall-Winter foothills Half 0.60 1 AksToKa
> 6 2. Spring foothills Half 1.00 1 AksToKa
> 7 3. Summer foothills Half 1.25 1 AksToKa
> 8 4. Fall foothills Half 1.50 1 AksToKa
> 9 1. Fall-Winter foothills Meat 0.60 1 AksToKa
> 10 2. Spring foothills Meat 1.00 1 AksToKa
>
>
> > kew11.nlme$call
> nlme.formula(model = lw ~ SSasympOff(ageyr2, mw, lgr, age0),
> data = kew, fixed = list(mw + lgr + age0 ~ Season * MBreed +
> ecoreg + ecoreg:Season), random = list(farm = list(mw ~
> 1, lgr ~ 1), sheep = list(mw ~ 1)), start = c(fixef
> (kew8.nlme)[1:15],
> 0, 0, 0, 0, 0, 0, 0, 0, 0, fixef(kew8.nlme)[16:30], 0,
> 0, 0, 0, 0, 0, 0, 0, 0, fixef(kew8.nlme)[31:45], 0, 0,
> 0, 0, 0, 0, 0, 0, 0), correlation = corSymm())
>
>
> > levels(kew$Season)==levels(x.for.lw.pred$Season)
> [1] TRUE TRUE TRUE TRUE
>
>
> > kew[1:10,]
> Grouped Data: lw ~ ageyr2 | farm
> X sheep doe ageyr2 MBreEco Season
> ecoreg
> farm village Breed MBreed date doy lw ebw
> 1 1 1 1 2.500000 Meatsemidesert 1. Fall-Winter semidesert
> AksToKa Aksenger KZP Meat 11/23/02 327 63.8 55.63211
> 2 2 1 139 2.883333 Meatsemidesert 2. Spring semidesert
> AksToKa Aksenger KZP Meat 4/10/03 100 51.7 44.53119
> 3 3 1 250 3.191667 Meatsemidesert 3. Summer semidesert
> AksToKa Aksenger KZP Meat 7/30/03 211 58.3 50.58624
> 4 4 1 330 3.413889 Meatsemidesert 4. Fall semidesert
> AksToKa Aksenger KZP Meat 10/18/03 291 59.9 52.05413
> 5 5 2 1 6.000000 Meatsemidesert 1. Fall-Winter semidesert
> AksToKa Aksenger KZP Meat 11/23/02 327 58.0 50.31101
> 6 6 2 139 6.383333 Meatsemidesert 2. Spring semidesert
> AksToKa Aksenger KZP Meat 4/10/03 100 41.2 34.89817
> 7 7 2 250 6.691667 Meatsemidesert 3. Summer semidesert
> AksToKa Aksenger KZP Meat 7/30/03 211 53.3 45.99908
> 8 8 2 330 6.913889 Meatsemidesert 4. Fall semidesert
> AksToKa Aksenger KZP Meat 10/18/03 291 63.7 55.54037
> 9 9 3 1 4.000000 Meatsemidesert 1. Fall-Winter semidesert
> AksToKa Aksenger KZP Meat 11/23/02 327 62.3 54.25596
> 10 10 3 139 4.383333 Meatsemidesert 2. Spring semidesert
> AksToKa Aksenger KZP Meat 4/10/03 100 47.3 40.49450
> bcs prcfat cageyrs
> 1 2.75 26.533 -0.46162659
> 2 2.00 20.192 -0.07829326
> 3 2.50 24.478 0.23004008
> 4 2.50 24.414 0.45226230
> 5 2.50 24.490 3.03837341
> 6 1.75 18.337 3.42170674
> 7 2.25 22.403 3.73004008
> 8 3.25 31.087 3.95226230
> 9 2.50 24.318 1.03837341
> 10 2.00 20.368 1.42170675
>
> ______________________________________________
> R-help@stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide!
> http://www.R-project.org/posting-guide.html
>



R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html Received on Sat Mar 18 10:29:19 2006

This archive was generated by hypermail 2.1.8 : Sun 19 Mar 2006 - 18:09:12 EST