[R] Using corStruct in nlme

[Duncan Murdoch]

On 7/20/2006 7:16 PM, Thomas Lumley wrote:
> 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.

So that looks like an nlme problem, but there's nothing obviously wrong 
in the code.  Doug, can you spot what the problem might be?

Duncan Murdoch
> 
>  	-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