R – Error in contr.treatment(n = 0L): Not Enough Degrees of Freedom to Define Contrasts in nlme

degrees of freedomlme4-nlmelogisticr

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?

Best Answer

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! :)

Related Solutions

Solved – Non-linear mixed effects regression in R

I wanted to share some of the things I learned since asking this question. nlme seems a reasonably way to model non-linear mixed effects in R. Start with a simple base model:

library(nlme)
data <- groupedData(y ~ t | UID, data=data) ## not strictly necessary
initVals <- getInitial(y ~ SSlogis(t, Asym, xmid, scal), data = data)
baseModel<- nlme(y ~ SSlogis(t, Asym, xmid, scal),
    data = data,
    fixed = list(Asym ~ 1, xmid ~ 1, scal ~ 1),
    random = Asym + xmid + scal ~ 1|UID,
    start = initVals
)

Then use update to increase model complexity. The start parameter is slightly tricky to work with, it may take some tinkering to figure out the order. Note how the new fixed effect for var1 on Asym follows the regular fixed effect for Asym.

 nestedModel <- update(baseModel, fixed=list(Asym ~ var1, xmid ~ 1, scal ~ 1), start = c(fixef(baseModel)[1], 0, fixef(baseModel)[2], fixef(baseModel)[3]))

lme4 seemed more robust against the outliers in my dataset and seemed to offer more reliable convergence for the more complex models. However, it seems the downside is that the relevant likelihood functions need to be specified manually. The following is the logistic growth model with a fixed effect of var1 (binary) on Asym. You can add fixed effects on xmid and scal in a similar fashion. Note the strange way of specifying the model using a double formula as outcome ~ fixed effects ~ random effects.

library(lme4) ## careful loading nlme and lme4 concurrently
customLogitModel <- function(t, Asym, AsymVar1, xmid, scal) {
    (Asym+AsymVar1*var1)/(1+exp((xmid-t)/scal))
}

customLogitModelGradient <- deriv(
    body(customLogitModel)[[2]], 
    namevec = c("Asym", "AsymVar1", "xmid", "scal"), 
    function.arg=customLogitModel
)

## find starting parameters
initVals <- getInitial(y ~ SSlogis(t, Asym, xmid, scal), data = data)

# Fit the model
model <- nlmer(
    y ~ customLogitModelGradient(t=t, Asym, AsymVar1, xmid, scal, var1=var) ~ 
    # Random effects with a second ~
    (Asym | UID) + (xmid | UID) + (scal | UID), 
    data = data, 
    start = c(Asym=initVals[1], AsymVar1=0, xmid=initVals[2], scal=initVals[3])
)

Solved – group fixed-effects, not individual-fixed effects using plm in R

I have worked on similar projects and am confronting one right now. The way that we handle this is to put in a fixed effect for each village and then to cluster the standard errors by village. This is not a perfect solution, but is fairly standard practice.

The plm package in R and xtreg ..., fe command in Stata, and the traditional fixed effect (within) estimator are designed to follow individuals. I believe one of the names for the method that you want is called a hierarchical linear model.

The simplest implementation in R would be something like

myLM <- lm(y ~ x + v v.t*t, data=df)

where y is the outcome of interest, x is some set of controls, v is a factor variable for the villages, v.t is a binary (factor) variable indicating whether a village was treated, and t is an indicator for pre-post treatment.

For standard inference, it is typical and recommended to produce clustered standard errors use either the multiwayvcov package or clusterSEs package.

Another method for inference, and the preferred method in Bertrand, Duflo & Mullainathan, 2004 is to perform a placebo test, where you vary "treatment" across all villages, form an empirical CDF, and see where the effect of treatment for the truly treated village sits in that distribution. Note that this is roughly the same method recommended for inference with synthetic controls of Abadie, Diamond, and Hainmueller, and has ties back to Fisher's 1935 text.

Best Answer

Related Solutions

Solved – Non-linear mixed effects regression in R

Solved – group fixed-effects, not individual-fixed effects using plm in R

Related Question