I'm running a self-starting 3-parameter logistic model using nlme and would like clarification in how starting values should be chosen when there are covariates in the model.
Create an example data frame:
sampledf <- tibble(
age = rep(c(0, 3, 6, 9, 12, 15, 18, 21), 4),
mass = c(1, 4, 8, 12, 16, 20, 22, 23,
1, 5, 8, 13, 17, 21, 24.5, 24.5,
1, 2.5, 6.5, 10, 14, 18.5, 20.5, 21,
1, 5, 8, 13, 17, 21, 24.5, 24.5),
treatment = c(rep(0:1, each=8, 2)),
year = c(rep(0:1, each=16))) %>%
slice(rep(1:n(), each = 3))
id = c("rae", "raphael", "reshmi", "rudy", "ruth", "roland", "remy", "roxy", "red", "randy", "raj", "roy")
sampledf <- sampledf %>%
arrange(age) %>%
mutate(id = rep(id, 8),
id = as.factor(id), year = as.factor(year), treatment=as.factor(treatment))
Plot to visualize sigmoidal growth curves:
ggplot(sampledf, aes(age, mass, col=treatment, linetype=year)) +
geom_point() +
geom_smooth(method="gam", formula = y ~ s(x, k = 4))
Because the model takes 3 starting values for each parameter {e.g., Asym=c(23.0, 26.1, 24.5)}, does this mean that the first two starting values in the list should be the estimated values from treatment group, while the third should be the estimated value from both years averaged together? Or, should the third value be the estimate for 'treatment=0, year=0'? Further, in the model summary, how is the intercept determined, given that there are 3 fixed effect estimates per parameter (e.g., Intercept.Asym, Xmid.Asym, Scal.Asym)? Is the 'treatment=0, year=0' group the reference point, even though there is no estimate which represents 'treatment=0, year=1'?
model1 <- nlme(mass ~ SSlogis(age, Asym, xmid, scal),
fixed=list(Asym + xmid + scal ~ treatment + year),
random = Asym + xmid ~ 1 | id,
start=list(fixed=c(Asym=c(23.0, 26.1, 24.5),
xmid=c(9.33, 9.21, 9.1),
scal=c(3.63, 3.85, 3.7))),
data=sampledf)
summary(model1)
How I determined which starting values to use:
sampledf %>%
group_by(treatment) %>%
do(fit = nls(mass ~ SSlogis(age, Asym, xmid, scal),
data = .)) %>%
tidy(fit) %>%
select(treatment, term, estimate) %>%
spread(term, estimate)