[R-sig-ME] Multiple covariates in SSbiexp model and an error message

[Jesus Frias]

Hi Andrew,

I think you just needed a little change. I've changed slightly your 
code to show you a plotting method in nlsList and nlme that allows to 
inspect graphically the covariates (is really nice, at least for me). 
You should be able to use the coef() or ranef() to extract the random 
effects and test for covariate dependence before you build your final 
model (two-stage). This is the main reason why I love the suite of 
nlList, gnls and nlme in R.

I've modified the call to update so that it includes one dependence of 
A1 with the temperature, but it is easily extended to include others. 
It works in my computer, but send me an email it it doesn't in yours.


library(nlme)
temp<-rnorm(66, mean = 3, sd = 0.5)
temp<-data.frame(temp)
temp
Indometh1<-cbind(Indometh, temp)

Indometh1$temp1 <- 10^-(Indometh1$temp)

Indometh1<-groupedData(conc ~ time + temp1 | Subject, data = Indometh1)

fm1Indom.lis <- nlsList(conc ~ SSbiexp(time, A1, lrc1, A2, lrc2), data 
= Indometh1)

plot(coef(fm1Indom.lis,augFrame=T),A1~temp)
plot(coef(fm1Indom.lis,augFrame=T),A2~temp)


fm1Indom.nlme <- nlme(fm1Indom.lis, random = pdDiag(A1 + lrc1 + A2 + 
lrc2~1))
summary(fm1Indom.nlme)
plot(ranef(fm1Indom.nlme,augFrame=T),A1~temp)
plot(ranef(fm1Indom.nlme,augFrame=T),A2~temp)


fm2Indom.nlme <- update(fm1Indom.nlme,
                        fixed = list(A1~temp,
                          lrc1~1,
                          A2~1,
                          lrc2~1),
                        start=c(2.8328430, 0,0.7731807, 0.4625497, 
-1.3422416))


regards,

Jesus


Andrew Close wrote:


>Dear all,
>
>I wanted to experiment with the SSbiexp model but include multiple 
covariates. To this end I have been following the presented in P&B for 
the CO2 uptake model
>
>I was presented with the error message:
>
>### Error in eval(expr, envir, enclos) : object 'temp1' not found ###
>
>The following code illustrates my modelling process:
>I attached a new column to the Indometh data from randomly generated 
numbers to represent temperature.
>library(nlme)
>temp<-rnorm(66, mean = 3, sd = 0.5)
>temp<-data.frame(temp)
>temp
>Indometh1<-cbind(Indometh, temp)
>
>I converted the temperature values to the same order of magnitude using 
the following line of code.
>> Indometh1$temp1 <- 10^-(Indometh1$temp)
>
>converted the data to a grouped object and began the model fitting 
process:
>> Indometh1<-groupedData(conc ~ time + temp1 | Subject, data = Indometh1)
>
>> fm1Indom.lis <- nlsList(conc ~ SSbiexp(time + temp, A1, lrc1, A2, 
lrc2), data = Indometh1)
>
>> fm1Indom.nlme <- nlme(fm1Indom.lis, random = pdDiag(A1 + lrc1 + A2 + 
lrc2~1))
>
>> summary(fm1Indom.nlme)
>
>Nonlinear mixed-effects model fit by maximum likelihood
>  Model: conc ~ SSbiexp(time + H, A1, lrc1, A2, lrc2)
> Data: Indometh1.g
>        AIC       BIC   logLik
>  -91.35715 -71.65026 54.67858
>Random effects:
> Formula: list(A1 ~ 1, lrc1 ~ 1, A2 ~ 1, lrc2 ~ 1)
> Level: Subject
> Structure: Diagonal
>               A1      lrc1        A2         lrc2  Residual
>StdDev: 0.5640351 0.1556456 0.1122092 9.439076e-06 0.0815653
>Fixed effects: list(A1 ~ 1, lrc1 ~ 1, A2 ~ 1, lrc2 ~ 1)
>          Value Std.Error DF   t-value p-value
>A1    2.8328430 0.2613774 57 10.838132   0e+00
>lrc1  0.7731807 0.1095290 57  7.059140   0e+00
>A2    0.4625497 0.1132419 57  4.084615   1e-04
>lrc2 -1.3422416 0.2309100 57 -5.812833   0e+00
> Correlation:
>     A1     lrc1   A2
>lrc1  0.057
>A2   -0.102  0.634
>lrc2 -0.140  0.581  0.834
>Standardized Within-Group Residuals:
>       Min         Q1        Med         Q3        Max
>-3.1672960 -0.3592548 -0.1301288  0.3445645  2.9932219
>Number of Observations: 66
>Number of Groups: 6
>
>##
>Using the estimates from the above model I tried to generate new 
starting values for the fixed effects. My aim to see how the decay may 
vary according to time and the newly generated "temp1" variable.
>
>> fm2Indom.nlme <- update(fm1Indom.nlme, fixed = list(A1~time+temp, 
lrc1~1 + A2~time+temp + lrc2~1),
>+ start=c(2.8328430, 0, 0, 0.7731807, 0, 0.4625497, 0, 0, -1.3422416, 
0))
>
>I was presented with the the error message:
>
>### Error in eval(expr, envir, enclos) : object 'temp1' not found ###
>
>Can anyone clarify if it possible to generate biexp models with 
multiple covariates and if so what is it that I am doing wrong?
>
>Thank you for your patience.
>
>Best wishes
>
>Andrew Close
><mailto:a.j.close at ncl.ac.uk>
>
>_______________________________________________
>R-sig-mixed-models at r-project.org mailing list
>https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>

--------------------------------------------------
Jesús María Frías Celayeta
Head of the Dept of Food Science 
School of Food Sci. and Env. Health 
DIT 
t + 353 1 402 4459 
f + 353 1 402 4495
http://fseh.dit.ie/o4/StaffListing/JesusFrias.html

This message has been scanned for content and viruses by the DIT Information Services E-Mail Scanning Service, and is believed to be clean. http://www.dit.ie