[R] nlme predict with se?

[Berton Gunter]

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,

-- Bert Gunter
Genentech Non-Clinical Statistics
South San Francisco, CA
 
"The business of the statistician is to catalyze the scientific learning
process."  - George E. P. Box
 
 

> -----Original Message-----
> From: r-help-bounces at stat.math.ethz.ch 
> [mailto:r-help-bounces at stat.math.ethz.ch] On Behalf Of Dr. 
> Emilio A. Laca
> Sent: Friday, March 17, 2006 1:55 PM
> To: R-help at 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 at 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
>