[R-sig-ME] Generalized least squares program

[Gang Chen]

No, update() doesn't make much difference in my case either. I've  
been reluctant to dive into the source code, but I may have to do so  
if no better choices turn up.

Thanks again,
Gang


On Jul 16, 2008, at 4:32 PM, Steven McKinney wrote:

> Hi Gang,
>
> Sorry, I should have said "update()" instead
> of "predict()".
>
> However, a timing test on a small example shows no
> considerable gain to updating an existing object
> compared to refitting a new object.  This may be different
> for your scenario.  If not, then it's an exercise of
> diving into the code base to figure out how to adapt
> the update.gls() function or related fitting
> functions so you can run your extra fits without
> recomputing everything.
>
>
>> require("nlme")
>>
>> Ovary <- Ovary
>> set.seed(123)
>> Ovary$jitterFollicles <- Ovary$follicles + sample(((-2):2), nrow 
>> (Ovary), replace = TRUE)
>>
>> fm1 <- gls(follicles ~ sin(2*pi*Time) + cos(2*pi*Time), Ovary,
> +            correlation = corAR1(form = ~ 1 | Mare))
>> fm2 <- update(fm1, model. = jitterFollicles ~ .)
>>
>> fm3 <-  gls(jitterFollicles ~ sin(2*pi*Time) + cos(2*pi*Time), Ovary,
> +            correlation = corAR1(form = ~ 1 | Mare))
>>
>> all.equal(fitted(fm2), fitted(fm3))
> [1] TRUE
>>
>> system.time(
> +             for(i in 1:100) {
> +               fm2 <- update(fm1, model. = jitterFollicles ~ .)
> +             }
> +             )
>    user  system elapsed
>  12.022   5.238  17.264
>> system.time(
> +             for(i in 1:100) {
> +               fm3 <-  gls(jitterFollicles ~ sin(2*pi*Time) + cos 
> (2*pi*Time), Ovary,
> +                           correlation = corAR1(form = ~ 1 | Mare))
> +             }
> +             )
>    user  system elapsed
>  11.898   5.227  17.126
>>
>
> Steve McKinney
>
>
> -----Original Message-----
> From: Gang Chen [mailto:gangchen at mail.nih.gov]
> Sent: Wed 7/16/2008 11:34 AM
> To: Steven McKinney
> Cc: r-sig-mixed-models at r-project.org
> Subject: Re: [R-sig-ME] Generalized least squares program
>
> Thanks for the suggestion!
>
> It seems to me predict() is the opposite of what I'm looking for.
> With an extended linear model
>
> Y = X*b + e, e ~ N(0, R s^2), R s^2 is the  covariance matrix of the
> residual term e
>
> both X and the structure of R remain the same across all the loops in
> my case, and Y is different and so are the estimates of b and R s^2
> from one iteration to another.
>
> Gang
>
>
> On Jul 16, 2008, at 1:29 PM, Steven McKinney wrote:
>
>> Hi Gang,
>>
>> Have you tried the predict() function?
>> You should be able to hand off new data
>> to your saved fitted model object.
>> See
>> ?predict.gls
>> for an example.
>>
>> HTH
>>
>> Steve McKinney
>>
>> -----Original Message-----
>> From: r-sig-mixed-models-bounces at r-project.org on behalf of Gang Chen
>> Sent: Wed 7/16/2008 8:21 AM
>> To: r-sig-mixed-models at r-project.org
>> Subject: [R-sig-ME] Generalized least squares program
>>
>> I'm using function gls in nlme to run some extended linear model with
>> generalized least squares many many times in a loop with the same
>> design matrix X and same residual covariance structure but different
>> values for the response variable. A lot of time is wasted in
>> reformulating the same model structures and overhead matrix
>> calculations. Such a waste makes the runtime range from one day to
>> one month. Are there some slots available in gls that would allow me
>> to avoid the unnecessary waste? If not, is there a counterpart of gls
>> for extended linear model with generalized least squares available in
>> lme4 that would give such slots?
>>
>> Thanks,
>> Gang
>>
>> _______________________________________________
>> R-sig-mixed-models at r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>>
>
>