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.
Thanks,
Katie Grieve
Quantitative Ecology & Resource Management
University of Washington
Thanks for providing such a self-contained example by which 'nlme' crashes R. Could you please also give us 'sessionInfo()'? I don't have time to test it myself now, but perhaps if you identify your platform, you might interest someone else in checking it. I'm sorry I couldn't be more helpful. Spencer Graves grieve at 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 longeridentical 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. > > Thanks, > Katie Grieve > > Quantitative Ecology & Resource Management > University of Washington > > ______________________________________________ > R-help at 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.
On 7/18/2006 2:28 PM, grieve at 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