3

I have a treatment (2 groups) for which I want to compare parameter estimate of a sigmoid curve (i.e., do any of the three parameters differ significantly between my two treatment group).

I am following this example (Compare non-linear model parameter estimates between conditions), which is also explained in the book "Nonlinear Regression R provides" p.115.

Therefore, I am trying to fit two models: one where all parameters are assumed to be identical for group1 and group2 and another one where all parameters vary between group1 and group2. Then I planned to compared them with a F-test (anova, again as suggested in the book). Without random effect it works fine, but I am struggling to do so when there is the random effect.

Basically, I struggle to fit the second model, where I allow my treatment group to vary in all parameters, while assigning the random effect on all three parameters.

Here is a reproducible example:

library(nlme)
#add fake treatment
Loblolly$mygroup = c(rep(1,84/2), rep(2,84/2))
head(Loblolly)

#works
fm1 = nlme(height ~ SSasymp(age, Asym, R0, lrc), data = Loblolly,
fixed = Asym + R0 + lrc ~ 1,
random = Asym ~ 1 | Seed,
start = c(Asym = 103, R0 = -8.5, lrc = -3.3))


#trying to allow the group to vary in all parameters, but it does not work: 
fm2 <- nlme(height ~ SSasymp(age, Asym[mygroup], R0[mygroup], lrc[mygroup]), data = Loblolly,
fixed = Asym + R0 + lrc ~ 1,
random = Asym ~ 1 | Seed,
start = list(fixed = as.data.frame(rbind(fixef(fm1), fixef(fm1)))))

I am working with SSlogis and not SSasymp, but I don't this this is important. Here is a small reproducible example of what I am aiming at but which does not contain the random effect:

#first model(m0)
m0 = nls(height ~ SSlogis(age, Asym, xmid, scal), Loblolly)
#second model (m12)
m1 = nls(height ~ SSlogis(age, Asym, xmid, scal), Loblolly, subset = mygroup == 1) 
m2 = nls(height ~ SSlogis(age, Asym, xmid, scal), Loblolly, subset = mygroup == 2)
m12 = nls(height ~ SSlogis(age, Asym[mygroup], xmid[mygroup], scal[mygroup]),
          data=Loblolly, 
          start = as.data.frame(rbind(coef(m1), coef(m2)))) #give starting points
#model comparison
anova(m0, m12)

Any help would be really greatly appreciated.

miki
  • 294

1 Answers1

1

Obviously, you don't understand the nlme documentation.

fm2 <- update(fm1, 
              #linear models for parameters
              fixed = Asym + R0 + lrc ~ mygroup,
              #starting values with intercept and slope 
              #based on treatment contrasts
              start = as.numeric(rbind(c(Asym = 103, R0 = -8.5, lrc = -3.3), 0)))

summary(fm2)


anova(fm1, fm2)


Model df      AIC      BIC    logLik   Test  L.Ratio p-value


#fm1     1  5 239.4856 251.6397 -114.7428

#fm2     2  8 240.3754 259.8219 -112.1877 1 vs 2 5.110246  0.1639

I assume you intend to only have a random effect for the asymptote and not for the other parameters.

Roland
  • 6,611