Using corStruct in nlme

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

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