[R-sig-ME] Generalized least squares program

[Gang Chen]

Dr. Bates,

Thanks a lot for the clarification!

Is there any plan to have a gls method implemented in lme4 so that I  
could take advantage of the refit function?

Also does lmer in lme4 allow for specifying covariance structure for  
the residuals? If yes, is there a way to trick lmer so that I can run  
gls by somehow zeroing the across-group random effect Z*b in model y  
= X*beta + Z*b + e?

Thanks,
Gang


On Jul 17, 2008, at 12:40 PM, Douglas Bates wrote:

> In general update applied to a fitted model object modifies the
> original call then evaluates the modified call.  Hence there will be
> no gain in speed.  The refit function in the lme4 package replaces the
> response and refits the model with the new response.  If there was
> such a function for a gls object that would be what you would want to
> use.  However, there is no refit method for gls objects.
>
> On Wed, Jul 16, 2008 at 4:23 PM, Gang Chen <gangchen at mail.nih.gov>  
> wrote:
>> 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
>>>>
>>>
>>>
>>
>> _______________________________________________
>> R-sig-mixed-models at r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>>