[R] Using corStruct in nlme

[Thomas Lumley]

On Thu, 20 Jul 2006, grieve at u.washington.edu wrote:

> 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.


Running these under Valgrind they both show the same problem

==10878== Invalid read of size 4
==10878==    at 0x1CCDB9EB: mult_mat (matrix.c:84)
==10878==    by 0x1CCD81B9: corStruct_recalc (corStruct.c:72)
==10878==    by 0x1CCDC06E: nlme_wtCorrAdj (nlme.c:152)
==10878==    by 0x1CCDC4F7: fit_nlme (nlme.c:347)
==10878==    by 0x80A518E: do_dotCode (dotcode.c:1777)
==10878==    by 0x80C45A9: Rf_eval (eval.c:444)
==10878==    by 0x80C600A: do_set (eval.c:1340)
==10878==    by 0x80C44A2: Rf_eval (eval.c:427)
==10878==    by 0x80C6097: do_begin (eval.c:1104)
==10878==    by 0x80C44A2: Rf_eval (eval.c:427)
==10878==    by 0x80C61AD: do_repeat (eval.c:1066)
==10878==    by 0x80C44A2: Rf_eval (eval.c:427)
==10878==  Address 0x1D4A60DC is 4 bytes after a block of size 6048 alloc'd
==10878==    at 0x1B909B71: calloc (vg_replace_malloc.c:175)
==10878==    by 0x80E09D5: R_chk_calloc (memory.c:2315)
==10878==    by 0x1CCDC415: fit_nlme (nlme.c:113)
==10878==    by 0x80A518E: do_dotCode (dotcode.c:1777)
==10878==    by 0x80C45A9: Rf_eval (eval.c:444)
==10878==    by 0x80C600A: do_set (eval.c:1340)
==10878==    by 0x80C44A2: Rf_eval (eval.c:427)
==10878==    by 0x80C6097: do_begin (eval.c:1104)
==10878==    by 0x80C44A2: Rf_eval (eval.c:427)
==10878==    by 0x80C61AD: do_repeat (eval.c:1066)
==10878==    by 0x80C44A2: Rf_eval (eval.c:427)
==10878==    by 0x80C6097: do_begin (eval.c:1104)


They also both go on to crash so hard that Valgrind crashes and says

--9895-- INTERNAL ERROR: Valgrind received a signal 11 (SIGSEGV) - exiting
--9895-- si_code=1 Fault EIP: 0xB00313AD; Faulting address: 0x3C439F95
--9895--   esp=0xB0653E48
valgrind: the `impossible' happened:
    Killed by fatal signal


Ouch.

 	-thomas



> # 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 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
>>
>
> ______________________________________________
> 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.
>

Thomas Lumley			Assoc. Professor, Biostatistics
tlumley at u.washington.edu	University of Washington, Seattle