Re: [R] Using corStruct in nlme

From: <grieve_at_u.washington.edu>
Date: Fri 21 Jul 2006 - 07:23:49 EST


Duncan, i could not find one seed which caused both corGaus and corExp to crash in R-patched. but, i found a seed that caused each to fail individually. Thanks for your help.

# For corExp:

set.seed(26)
for(i in 1:length(CO2$conc)){

        CO2$conc[i]<-(CO2$conc[i]+rnorm(1))
}

fm1CO2.lis <- nlsList(SSasympOff, CO2)
fm1CO2.nlme <- nlme(fm1CO2.lis, control = list(tolerance = 1e-2))
fm2CO2.nlme <- update(fm1CO2.nlme, random = Asym + lrc ~ 1)
CO2.nlme.var <- update(fm2CO2.nlme,
  fixed = list(Asym ~ Type * Treatment, lrc + c0 ~ 1),   start = c(32.412, 0, 0, 0, -4.5603, 49.344), weights=varConstPower(fixed=list(const=0.1, power=1)), verbose=T)   

CO2.nlme.exp<-update(CO2.nlme.var, correlation=corExp(form=~as.numeric(conc)|Plant,nugget=F), data=CO2)

# For corGaus:

set.seed(4)
for(i in 1:length(CO2$conc)){

        CO2$conc[i]<-(CO2$conc[i]+rnorm(1))
}

fm1CO2.lis <- nlsList(SSasympOff, CO2)
fm1CO2.nlme <- nlme(fm1CO2.lis, control = list(tolerance = 1e-2))
fm2CO2.nlme <- update(fm1CO2.nlme, random = Asym + lrc ~ 1)
CO2.nlme.var <- update(fm2CO2.nlme,
  fixed = list(Asym ~ Type * Treatment, lrc + c0 ~ 1),   start = c(32.412, 0, 0, 0, -4.5603, 49.344), weights=varConstPower(fixed=list(const=0.1, power=1)), verbose=T)   

CO2.nlme.gauss<-update(CO2.nlme.var, correlation=corGaus(form=~as.numeric(conc)|Plant,nugget=F), data=CO2)

On Thu, 20 Jul 2006, Duncan Murdoch wrote:

> On 7/18/2006 2:28 PM, grieve@u.washington.edu wrote:

>> I am having trouble fitting correlation structures within nlme. I would
>> like to fit corCAR1, corGaus and corExp correlation structures to my data.
>> I either get the error "step halving reduced below minimum in pnls step"
>> or alternatively R crashes.
>>
>> My dataset is similar to the CO2 example in the nlme package. The one
>> major difference is that in my case the 'conc' steps are not the same for
>> each 'Plant'. I have replicated the problem using the CO2 data in nlme
>> (based off of the Ch08.R script).
>>
>> This works (when 'conc' is the same for each 'Plant':
>>
>> (fm1CO2.lis <- nlsList(SSasympOff, CO2))
>> (fm1CO2.nlme <- nlme(fm1CO2.lis, control = list(tolerance = 1e-2)))
>> (fm2CO2.nlme <- update(fm1CO2.nlme, random = Asym + lrc ~ 1))
>> CO2.nlme.var <- update(fm2CO2.nlme,
>> fixed = list(Asym ~ Type * Treatment, lrc + c0 ~ 1),
>> start = c(32.412, 0, 0, 0, -4.5603, 49.344),
>> weights=varConstPower(fixed=list(const=0.1, power=1)), verbose=T)
>>
>> CO2.nlme.CAR<-update(CO2.nlme.var, corr=corCAR1())
>>
>> CO2.nlme.gauss<-update(CO2.nlme.var,
>> correlation=corGaus(form=~as.numeric(conc)|Plant,nugget=F), data=CO2)
>>
>> CO2.nlme.exp<-update(CO2.nlme.var,
>> correlation=corExp(form=~as.numeric(conc)|Plant,nugget=F), data=CO2) But,
>> if i change each of the 'conc' numbers slightly so that they are no longer
>> identical between subjects i can only get the corCAR1 correlation to work
>> while R crashes for both corExp and corGaus:
>>
>> for(i in 1:length(CO2$conc)){
>> CO2$conc[i]<-(CO2$conc[i]+rnorm(1))
>> }
>>
>> (fm1CO2.lis <- nlsList(SSasympOff, CO2))
>> (fm1CO2.nlme <- nlme(fm1CO2.lis, control = list(tolerance = 1e-2)))
>> (fm2CO2.nlme <- update(fm1CO2.nlme, random = Asym + lrc ~ 1))
>> CO2.nlme.var <- update(fm2CO2.nlme,
>> fixed = list(Asym ~ Type * Treatment, lrc + c0 ~ 1),
>> start = c(32.412, 0, 0, 0, -4.5603, 49.344),
>> weights=varConstPower(fixed=list(const=0.1, power=1)), verbose=T)
>>
>> CO2.nlme.CAR<-update(CO2.nlme.var, corr=corCAR1())
>>
>> CO2.nlme.gauss<-update(CO2.nlme.var,
>> correlation=corGaus(form=~as.numeric(conc)|Plant,nugget=F), data=CO2)
>>
>> CO2.nlme.exp<-update(CO2.nlme.var,
>> correlation=corExp(form=~as.numeric(conc)|Plant,nugget=F), data=CO2) I
>> have read Pinheiro & Bates (2000) and i think that it should be possible
>> to fit these correlation structures to my data, but maybe i am mistaken.
>>
>> I am running R 2.3.1 and have recently updated all packages.
> 
> I reproduced this once in R-patched, but since then have been unable to do 
> so. I can reproduce it reliably with "set.seed(1)" at the start in R 2.3.1. 
> So it looks to me as though we've probably done something to make the error 
> less likely, but not completely fixed it.
> 
> If you can find a script (including set.seed() to some value) that reliably 
> causes a crash in R-patched, could you let me know?
> 
> You can get R-patched from CRAN in the bin/windows/base directory.
> 
> Duncan Murdoch
>

______________________________________________
R-help@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 and provide commented, minimal, self-contained, reproducible code. Received on Fri Jul 21 07:31:37 2006

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.1.8, at Fri 21 Jul 2006 - 10:18:27 EST.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help. Please read the posting guide before posting to the list.