1
$\begingroup$

We are attempting to model and compare logistic growth over time for 6 different treatments using nlme. So far, we have successfully added random effects of individuals. However, when we try to add any combination or number of fixed effects of treatments we get "Error in contr.treatment(n = 0L) : not enough degrees of freedom to define contrasts". We get the same error when adding fixed effects to a model containing any combination or number of random effects. Below is a sample of our data and the code that is producing the error. In our full data set each treatment has 12 individuals, measured at 12 time steps, the sample contains measurements of one individual per treatment. Searching the error on google only revealed a few hits that I wasn't able to gain any clarity from.

pacman::p_load(
  tidyverse,
  nlme)
df <- structure(list(Individual = c("SO_BF_05_MF_G", "SO_BF_05_MF_G", "SO_BF_05_MF_G", "SO_BF_05_MF_G", "SO_BF_05_MF_G", "SO_BF_05_MF_G", "SO_BF_05_MF_G", "SO_BF_05_MF_G", "SO_BF_05_MF_G", "SO_BF_05_MF_G", "SO_BF_05_MF_G", "PP_BF_05_MF_D", "PP_BF_05_MF_D", "PP_BF_05_MF_D", "PP_BF_05_MF_D", "PP_BF_05_MF_D", "PP_BF_05_MF_D", "PP_BF_05_MF_D", "PP_BF_05_MF_D", "PP_BF_05_MF_D", "PP_BF_05_MF_D", "PP_BF_05_MF_D", "PP_BF_05_MF_D", "BN_BF_02_MF_A", "BN_BF_02_MF_A", "BN_BF_02_MF_A", "BN_BF_02_MF_A", "BN_BF_02_MF_A", "BN_BF_02_MF_A", "BN_BF_02_MF_A", "BN_BF_02_MF_A", "BN_BF_02_MF_A", "BN_BF_02_MF_A", "BN_BF_02_MF_A", "BN_BF_02_MF_A", "TD_BF_02_MF_H", "TD_BF_02_MF_H", "TD_BF_02_MF_H", "TD_BF_02_MF_H", "TD_BF_02_MF_H", "TD_BF_02_MF_H", "TD_BF_02_MF_H", "TD_BF_02_MF_H", "TD_BF_02_MF_H", "TD_BF_02_MF_H", "TD_BF_02_MF_H", "TD_BF_02_MF_H", "E6_BF_13_MF_L", "E6_BF_13_MF_L", "E6_BF_13_MF_L", "E6_BF_13_MF_L", "E6_BF_13_MF_L", "E6_BF_13_MF_L", "E6_BF_13_MF_L", "E6_BF_13_MF_L", "E6_BF_13_MF_L", "E6_BF_13_MF_L", "E6_BF_13_MF_L", "E6_BF_13_MF_L", "PO_BF_08_MF_L", "PO_BF_08_MF_L", "PO_BF_08_MF_L", "PO_BF_08_MF_L", "PO_BF_08_MF_L", "PO_BF_08_MF_L", "PO_BF_08_MF_L", "PO_BF_08_MF_L", "PO_BF_08_MF_L", "PO_BF_08_MF_L", "PO_BF_08_MF_L", "PO_BF_08_MF_L"), Treatment = c("O2", "O2", "O2", "O2", "O2", "O2", "O2", "O2", "O2", "O2", "O2", "PP", "PP", "PP", "PP", "PP", "PP", "PP", "PP", "PP", "PP", "PP", "PP", "PV", "PV", "PV", "PV", "PV", "PV", "PV", "PV", 
"PV", "PV", "PV", "PV", "S1", "S1", "S1", "S1", "S1", "S1", "S1", 
"S1", "S1", "S1", "S1", "S1", "S2", "S2", "S2", "S2", "S2", "S2", 
"S2", "S2", "S2", "S2", "S2", "S2", "SO", "SO", "SO", "SO", "SO", 
"SO", "SO", "SO", "SO", "SO", "SO", "SO"), Time.days = c(0, 
5, 12, 19, 26, 40, 47, 58, 65, 72, 99, 0, 5, 12, 19, 26, 33, 
40, 47, 58, 65, 72, 99, 0, 5, 12, 19, 26, 33, 40, 47, 58, 65, 
72, 99, 0, 5, 12, 19, 26, 33, 40, 47, 58, 65, 72, 99, 0, 5, 12, 
19, 26, 33, 40, 47, 58, 65, 72, 99, 0, 5, 12, 19, 26, 33, 40, 
47, 58, 65, 72, 99), Size = c(0, 0, 0, 0, 0, 0.0004444, 
0.256, 0.423, 2.511, 51.08, 55.935, 0, 0, 0, 0, 0.058, 0.149, 
1.711, 2.396, 14.747, 39.905, 51.462, 71.469, 0, 0, 1.111e-05, 
0, 1.094, 8.351, 15.96, 18.982, 52.926, 67.767, 70.759, 74.131, 
0, 0, 0.028, 0.971, 0.625, 2.302, 6.67, 15.891, 19.676, 51.914, 
72.486, 33.264, 0, 0, 0, 0, 0.078, 0.271, 0.638, 2.165, 24.953, 
39.935, 56.932, 64.002, 0, 0, 0, 0, 0, 0.271, 0.042, 0.621, 0.338, 
18.682, 59.92, 0)), class = c("grouped_df", "tbl_df", "tbl", 
"data.frame"), row.names = c(NA, -71L))
# get start values
df_nls <- nls(Size ~ SSlogis(Time.days, Asymp, xmid, scal), data = df)
# make growth model
growth.model = function(Time.days, Asymp, xmid, scal){Asymp/(1 + exp(-(Time.days - xmid)/scal))}
# nlme with fixed effect (f) of treatment on Asymptote (A), random effect of individual on asymptote
## x is nls data, y is all data
nlme_fA <- function(x,y) {
  nlme_start <- summary(x)$parameters
  Asymp_st <- nlme_start[,"Estimate"][1]
  xmid_st <- nlme_start[,"Estimate"][2]
  scal_st <- nlme_start[,"Estimate"][3]
  start_param <- c(Asymp_st, xmid_st, scal_st)
  nlme(Size ~ growth.model(Time.days, Asymp, xmid, scal), fixed = list(Asymp ~ Treatment, xmid + scal ~ 1), random = Asymp ~ 1|Individual, data = y, start = start_param)
}

df_fA <- nlme_fA(df_nls, df)

I have tried:

  • various combinations of individual and multiple fixed effects and random effects
  • adding duplicates of the data with slight variation to see if the issue was number of replicates
  • using grouped data frames df_grp_Trt_code <- groupedData(Size ~ Time.days|Treatment/Individual, data = df)
  • using self start function SSlogis for nlme nlme(Size ~ SSlogis(Time.days, Asymp, xmid, scal), fixed = list(Asymp ~ Treatment, xmid + scal ~ 1), random = Asymp ~ 1|Individual, data = y)
  • various combinations of the above attempts

All gave the same error

Is the error a result of too few replicates?

If so, is there a way for us to compare logistic growth between treatments in nlme while taking into account random effects?

If it isn't a result of too few replicates, what does the error mean and how can we solve it?

$\endgroup$

1 Answer 1

2
$\begingroup$

The main problem is that you did not define Treatment as a factor.

Note: I removed your last entry since I assume it was an error

df2 <- df[-71,]

Define Treatment as a factor

df2$Treatment <- as.factor(df2$Treatment)

Create groupedData

df2G <- groupedData(Size ~ Time.days | Individual, data = df2)

Visualize

plot(df2G)

Fit logistic to each Individual

fitL <- nlsList(Size ~ SSlogis(Time.days, Asym, xmid, scal), data = df2G)

Fit NLME

fnm1 <- nlme(fitL, random = pdDiag(Asym + xmid + scal ~ 1))

Extract estimate for fixed effects

fxf <- fixef(fnm1)

Update model incorporating the effect of treatment

fnm2 <- update(fnm1, fixed = Asym + xmid + scal ~ Treatment, start = c(fxf[1], rep(0, 5), fxf[2], rep(0, 5), fxf[3], rep(0, 5)))

Test effect of each treatment

anova(fnm2)

Check model assumptions

plot(fnm2)

Hope this helps! :)

$\endgroup$
1
  • $\begingroup$ that helps A LOT! this has pretty much entirely cleared things up $\endgroup$ Commented Jun 2, 2020 at 19:43

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.