[R] lme speedup question

[Dick Beyer]

Thanks very much to the help I received on my lme speedup question.  Jake Bowers pointed out that memory fill might slow things down, so I will keep watch on that.  Renaud Lancelot pointed out using update instead of redoing lme in the inner loop (also to look at using gls as it might be faster).  That sped up things by about 9%.  Douglas Bates pointed out my unnecessary creation of structures in the inner loop, so I moved the creation of the data.frame to the outermost loop.  With these changes, I see about a 27% speedup, for which I am very grateful. 

My modified code 

for (perm in 1:216){
  vg <- vgAll[perm,]
  ge <- c(philmaanova.rloess$adjdata[1,1:24])
  dat <- data.frame(ge,vg=factor(vg),ag=factor(ag),dy=factor(dy),rp=factor(rp))
  dat$vgrp <- getGroups(dat, form = ~ 1|vg/rp, level = 2)
  ge.lme   <- lme(fixed=ge~vg+ag+dy, data=dat, random=~1|vgrp)

  for(j in 1:15000){
    dat$ge   <- c(philmaanova.rloess$adjdata[j,1:24])
    ge.lme   <- update(ge.lme, data=dat)
    cp1[j,1] <- philmaanova.rloess$cloneid[j]
    tmpInt   <- intervals(ge.lme,level=0.95,which="fixed")
    cp1[j,2] <- tmpInt$fixed[2,2]*2
    cp1[j,3] <- tmpInt$fixed[2,1]
    cp1[j,4] <- tmpInt$fixed[2,3]
  }
}

> I am hoping someone will be kind enough to have a look at the following piece
of code and tell me if there is a way to run lme() so it is a lot faster. The
inner loop, j in 1:15000, takes about 2 hrs on my 2.8GHz dual Xeon 4GB RAM
machine.  The timings I have done show the dominant execution time is in lme.
> 
> options(contrasts=c("contr.sum", "contr.sum"))
> getOption("contrasts")
> vg  <- c(1,2,1,2,2,1,2,1,1,2,1,2,2,1,2,1,1,2,1,2,2,1,2,1)
> ag  <- c(1,1,1,1,2,2,2,2,3,3,3,3,4,4,4,4,5,5,5,5,6,6,6,6)
> dy  <- c(1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2)
> rp  <- c(1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,3,3,3,3,3,3,3,3)
> 
> for (perm in 1:216){
>  vg <- vgAll[perm,]
>  for(j in 1:15000){
>   ge <- c(philmaanova.rloess$adjdata[j,1:24])
>   dat <-
data.frame(ge,vg=factor(vg),ag=factor(ag),dy=factor(dy),rp=factor(rp))
>   dat$vgrp <- getGroups(dat, form = ~ 1|vg/rp, level = 2)
>   ge.lme   <- lme(fixed=ge~vg+ag+dy, data=dat, random=~1|vgrp)
>   cp1[j,1] <- philmaanova.rloess$cloneid[j]
>   tmpInt   <- intervals(ge.lme,level=0.95,which="fixed")
>   cp1[j,2] <- tmpInt$fixed[2,2]*2
>   cp1[j,3] <- tmpInt$fixed[2,1]
>   cp1[j,4] <- tmpInt$fixed[2,3]
> }
> }

Thanks for all your help,
Dick
*******************************************************************************
Richard P. Beyer, Ph.D.	University of Washington
Tel.:(206) 616 7378	Env. & Occ. Health Sci. , Box 354695
Fax: (206) 685 4696	4225 Roosevelt Way NE, # 100
			Seattle, WA 98105-6099