Generalized least squares program

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