R help - Oct 2004 - gnls or nlme : how to obtain confidence intervals of fitted values

Hi

I use gnls to fit non linear models of the form y = alpha * x**beta 
(alpha and beta being linear functions of a 2nd regressor z i.e. 
alpha=a1+a2*z and beta=b1+b2*z) with variance function 
varPower(fitted(.)) which sounds correct for the data set I use.

My purpose is to use the fitted models for predictions with other sets 
of regressors x, z than those used in fitting. I therefore need to 
estimate y with (95%) confidence intervals.

Does any body knows how to do this with R ?

Thanks

Pinhiero and Bates (2000) Mixed-Effects Models in S and S-Plus 
(Springer) describe the use of "intervals" for that.  In
library(nlme),
R 1.9.1 for Windows, I found documentation for intervals, intervals.gls, 
intervals.lme, intervals.lmList, and gnls.  To an example in the 
documentation for gnls, I added intervals: 

 >      data(Soybean)
 >      # variance increases with a power of the absolute fitted values
 >      fm1 <- gnls(weight ~ SSlogis(Time, Asym, xmid, scal), Soybean,
+                  weights = varPower())
 >      summary(fm1)
Generalized nonlinear least squares fit
  Model: weight ~ SSlogis(Time, Asym, xmid, scal)
  Data: Soybean
       AIC      BIC    logLik
  983.7947 1003.900 -486.8974

Variance function:
 Structure: Power of variance covariate
 Formula: ~fitted(.)
 Parameter estimates:
    power
0.8815437

Coefficients:
        Value Std.Error  t-value p-value
Asym 17.35682 0.5226885 33.20682       0
xmid 51.87232 0.5916820 87.66925       0
scal  7.62053 0.1390958 54.78617       0

 Correlation:
     Asym  xmid
xmid 0.787     
scal 0.485 0.842

Standardized residuals:
         Min           Q1          Med           Q3          Max
-2.309670367 -0.646844555 -0.004897024  0.498606088  4.986727281

Residual standard error: 0.3662752
Degrees of freedom: 412 total; 409 residual
 > intervals(fm1)
Approximate 95% confidence intervals

 Coefficients:
         lower      est.     upper
Asym 16.329331 17.356822 18.384313
xmid 50.709200 51.872317 53.035434
scal  7.347094  7.620525  7.893957
attr(,"label")
[1] "Coefficients:"

 Variance function:
          lower      est.    upper
power 0.8369085 0.8815437 0.926179
attr(,"label")
[1] "Variance function:"

 Residual standard error:
    lower      est.     upper
0.3373374 0.3649392 0.3947995

      Is this what you want? 
      hope this helps. 
      spencer graves

Pierre MONTPIED wrote:
> Hi
>
> I use gnls to fit non linear models of the form y = alpha * x**beta 
> (alpha and beta being linear functions of a 2nd regressor z i.e. 
> alpha=a1+a2*z and beta=b1+b2*z) with variance function 
> varPower(fitted(.)) which sounds correct for the data set I use.
>
> My purpose is to use the fitted models for predictions with other sets 
> of regressors x, z than those used in fitting. I therefore need to 
> estimate y with (95%) confidence intervals.
>
> Does any body knows how to do this with R ?
>
> Thanks
>
> ______________________________________________
> 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

-- 
Spencer Graves, PhD, Senior Development Engineer
O:  (408)938-4420;  mobile:  (408)655-4567

Thanks Spencer but the intervals function gives confidence intervals of 
the parameters of the model not the predicted values. In the Soybean 
example it would be the CI of predicted weight for a given time, knowing 
all the parameters (Asym, xmid, scal, variance function and residual) 
and their distributions.

And what I need to calculate is precisely this CI.

My question is therefore is there an analytical way to calculate such 
CI, whatever the model, or could I try some randomizing techniques such 
as bootstrap or other ?

Thanks.


Spencer Graves wrote:
>      Pinhiero and Bates (2000) Mixed-Effects Models in S and S-Plus 
> (Springer) describe the use of "intervals" for that.  In 
> library(nlme), R 1.9.1 for Windows, I found documentation for 
> intervals, intervals.gls, intervals.lme, intervals.lmList, and gnls.  
> To an example in the documentation for gnls, I added intervals:
> >      data(Soybean)
> >      # variance increases with a power of the absolute fitted values
> >      fm1 <- gnls(weight ~ SSlogis(Time, Asym, xmid, scal), Soybean,
> +                  weights = varPower())
> >      summary(fm1)
> Generalized nonlinear least squares fit
>  Model: weight ~ SSlogis(Time, Asym, xmid, scal)
>  Data: Soybean
>       AIC      BIC    logLik
>  983.7947 1003.900 -486.8974
>
> Variance function:
> Structure: Power of variance covariate
> Formula: ~fitted(.)
> Parameter estimates:
>    power
> 0.8815437
>
> Coefficients:
>        Value Std.Error  t-value p-value
> Asym 17.35682 0.5226885 33.20682       0
> xmid 51.87232 0.5916820 87.66925       0
> scal  7.62053 0.1390958 54.78617       0
>
> Correlation:
>     Asym  xmid
> xmid 0.787     scal 0.485 0.842
>
> Standardized residuals:
>         Min           Q1          Med           Q3          Max
> -2.309670367 -0.646844555 -0.004897024  0.498606088  4.986727281
>
> Residual standard error: 0.3662752
> Degrees of freedom: 412 total; 409 residual
> > intervals(fm1)
> Approximate 95% confidence intervals
>
> Coefficients:
>         lower      est.     upper
> Asym 16.329331 17.356822 18.384313
> xmid 50.709200 51.872317 53.035434
> scal  7.347094  7.620525  7.893957
> attr(,"label")
> [1] "Coefficients:"
>
> Variance function:
>          lower      est.    upper
> power 0.8369085 0.8815437 0.926179
> attr(,"label")
> [1] "Variance function:"
>
> Residual standard error:
>    lower      est.     upper
> 0.3373374 0.3649392 0.3947995
>
>      Is this what you want?      hope this helps.      spencer graves
>
> Pierre MONTPIED wrote:
>
>> Hi
>>
>> I use gnls to fit non linear models of the form y = alpha * x**beta 
>> (alpha and beta being linear functions of a 2nd regressor z i.e. 
>> alpha=a1+a2*z and beta=b1+b2*z) with variance function 
>> varPower(fitted(.)) which sounds correct for the data set I use.
>>
>> My purpose is to use the fitted models for predictions with other 
>> sets of regressors x, z than those used in fitting. I therefore need 
>> to estimate y with (95%) confidence intervals.
>>
>> Does any body knows how to do this with R ?
>>
>> Thanks
>>
>> ______________________________________________
>> 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
>
>
>

On Mon, 4 Oct 2004, Pierre MONTPIED wrote:
> Thanks Spencer but the intervals function gives confidence intervals of the
> parameters of the model not the predicted values. In the Soybean example it
> would be the CI of predicted weight for a given time, knowing all the 
> parameters (Asym, xmid, scal, variance function and residual) and their 
> distributions.
>
> And what I need to calculate is precisely this CI.
>
> My question is therefore is there an analytical way to calculate such CI, 
> whatever the model, or could I try some randomizing techniques such as 
> bootstrap or other ?
>
I find I have a need for this too but for the fixed effects in a mixed 
model fitted with lme. I have fitted a polynomial to some longitudinal 
data, plus some other fixed effects and some random effects. I can use 
predict to get the fitted polynomial for particular values of the other 
predictors, but I would really like to put some confidence bounds on it. I 
am not wild about bootstrapping---just fitting the model can take 30 
minutes.

Any suggestions will be gratefully received (although on second thoughts 
if someone tells me to use SAS I will be less than grateful).

David Scott



_________________________________________________________________
David Scott	Department of Statistics, Tamaki Campus
 		The University of Auckland, PB 92019
 		Auckland	NEW ZEALAND
Phone: +64 9 373 7599 ext 86830		Fax: +64 9 373 7000
Email:	d.scott at auckland.ac.nz


Graduate Officer, Department of Statistics

I believe what Pierre means to be asking for are "prediction
intervals", or
possibly "tolerance intervals", which are different from
"confidence
intervals". A great reference is:

Hahn G and Meeker WQ: Statistical Intervals. John Wiley and Sons, New York,
1991

For example, in a one-sample Normal problem, \bar{X} \pm 2
\hat{sigma}/\sqrt{n} is an approximate 95% confidence interval, while
\bar{X} \pm 2 \hat{sigma} is an approximate 95% prediction interval
(assuming \bar{X} and \hat{sigma} are high precision estimates; if they are
not you probably want to consider tolerance intervals, which are discussed
in the above reference). 

Both intervals.lme{nlme} and estimable{gmodels} will give you confidence
intervals, but neither will give you prediction intervals. I am not aware of
any published R function that gives you prediction intervals or tolerance
intervals for lme models. It is not easy to write such a function for the
general case, but it may be relatively easy to write your own for special
cases of lme models.  

Jim 
> Message: 9
> Date: Mon, 4 Oct 2004 21:42:20 +1000
> From: Andrew Robinson <andrewr at uidaho.edu>
> Subject: Re: [R] gnls or nlme : how to obtain confidence intervals of
> 	fitted	values
> To: David Scott <d.scott at auckland.ac.nz>
> Cc: r-help at stat.math.ethz.ch, Spencer Graves <spencer.graves at
pdf.com>
> Message-ID: <20041004114220.GJ67507 at uidaho.edu>
> Content-Type: text/plain; charset=us-ascii
> 
> The function estimable() from the gmodels part of the gregmisc package
> will do this, if applied appropriately.
> 
> It allows for the estimation of arbitrary linear combinations of the
> parameter estimates of a model object, and calcualtes confidence
> intervals.
> 
> I hope that this helps,
> 
> Andrew
> 
> 
> 
> On Tue, Oct 05, 2004 at 12:31:48AM +1300, David Scott wrote:
> > On Mon, 4 Oct 2004, Pierre MONTPIED wrote:
> > 
> > >Thanks Spencer but the intervals function gives confidence
intervals of
> > >the parameters of the model not the predicted values. In the
Soybean
> > >example it would be the CI of predicted weight for a given time,
knowing 
> > >all the parameters (Asym, xmid, scal, variance function and
residual)
and 
> > >their distributions.
> > >
> > >And what I need to calculate is precisely this CI.
> > >
> > >My question is therefore is there an analytical way to calculate
such
CI, 
> > >whatever the model, or could I try some randomizing techniques
such as
> > >bootstrap or other ?
> > >

LEGAL NOTICE\ Unless expressly stated otherwise, this messag...{{dropped}}