I am interested in fitting both linear and nonlinear mixed-effects models. I would like to fit models to data that are grouped by two factors: yclass and d1. The variable y is a biological response measured once for every mice during the study. In the study, each mice was observed for at least two time points (i.e., 1, 3, 28, 90, 92, 162 days: yclass) with different doses (d1). The variables d2, d5, d9 are predictors.
But the nlme failed to run for nonlinear mixed model. I tried different models:
library(tidyverse)
library(nlme)
library(aomisc)
mydata <- data.frame(ID=c("Id001", "Id002", "Id003", "Id004", "Id005", "Id006", "Id007", "Id008", "Id009", "Id010", "Id011", "Id012", "Id013", "Id014", "Id015", "Id016", "Id017", "Id018", "Id019", "Id020", "Id021", "Id022", "Id023", "Id024", "Id025", "Id026", "Id027", "Id028", "Id029", "Id030", "Id031", "Id032", "Id033", "Id034", "Id035", "Id036", "Id037", "Id038", "Id039", "Id040", "Id041", "Id042", "Id043", "Id044", "Id045", "Id046", "Id047", "Id048", "Id049", "Id050", "Id051","Id052", "Id053", "Id054", "Id055", "Id056", "Id057", "Id058", "Id059", "Id060", "Id061", "Id062", "Id063", "Id064", "Id065", "Id066", "Id067"), yclass=c("Day-1", "Day-1", "Day-1", "Day-1", "Day-1", "Day-1", "Day-1", "Day-1", "Day-1", "Day-1", "Day-1", "Day-1", "Day-1", "Day-1", "Day-1", "Day-1", "Day-1", "Day-3", "Day-3", "Day-3", "Day-3", "Day-3", "Day-3", "Day-3", "Day-3", "Day-3", "Day-3", "Day-3", "Day-3", "Day-3", "Day-3", "Day-3", "Day-3", "Day-3", "Day-3", "Day-3", "Day-3", "Day-3", "Day-3", "Day-3", "Day-28", "Day-28", "Day-28", "Day-28", "Day-28", "Day-28", "Day-28", "Day-28", "Day-28", "Day-28", "Day-28", "Day-28", "Day-28", "Day-28", "Day-28", "Day-28", "Day-28", "Day-90", "Day-90", "Day-90", "Day-90", "Day-90", "Day-90", "Day-92", "Day-92", "Day-92", "Day-92"), d1=c("Dose-A", "Dose-B", "Dose-A", "Dose-A", "Dose-C", "Dose-A", "Dose-D", "Dose-A", "Dose-C", "Dose-A", "Dose-A", "Dose-C", "Dose-D", "Dose-A", "Dose-D", "Dose-D", "Dose-C", "Dose-E", "Dose-A", "Dose-E", "Dose-E", "Dose-C", "Dose-B", "Dose-B", "Dose-F", "Dose-F", "Dose-B", "Dose-A", "Dose-C", "Dose-A", "Dose-A", "Dose-C", "Dose-C", "Dose-D", "Dose-C", "Dose-C", "Dose-G", "Dose-G", "Dose-G", "Dose-G", "Dose-B", "Dose-E", "Dose-F", "Dose-B", "Dose-B", "Dose-F", "Dose-B", "Dose-F", "Dose-D", "Dose-D", "Dose-C", "Dose-D", "Dose-C", "Dose-C", "Dose-D", "Dose-A", "Dose-D", "Dose-A", "Dose-A", "Dose-A", "Dose-G", "Dose-G", "Dose-D", "Dose-A", "Dose-A", "Dose-G", "Dose-A"), d2=c(1, 22.5, 40, 12, 1, 64.2, 61, 22.1, 26.9, 55.6, 40, 30.2, 31, 1.69, 10, 32.7, 1.69, 20, 40, 40, 25, 40, 40, 25, 85, 40, 85, 40, 40, 61, 40, 61, 40, 61, 14, 14, 1, 74, 1, 7.8, 22.5, 22.5, 20, 85, 40, 40, 25, 25, 61, 1, 30.2, 22.1, 26.9, 61, 12, 12, 26.9, 25, 30, 14, 95, 31, 25, 26.9, 30.2, 14, 22.1), d5=c(453.1, 87.7, 50, 135, 453.1, 18, 25.6, 150, 152, 82, 50, 119, 84.6, 411, 99, 74, 411, 104, 50, 27.7, 86.9, 136, 27.7, 86.9, 27.4, 27.7, 27.4, 136, 100, 25.6, 50, 25.6, 50, 25.6, 338, 338, 24.7, 26, 2.74, 73.59, 87.7, 87.7, 104, 27.4, 27.7, 27.7, 86.9, 86.9, 25.6, 453.1, 141, 150, 152, 25.6, 135, 135, 152, 10.1, 11.42, 338, 23.8, 84.6, 10.1, 152, 141, 338, 150), d9=c(2, 3, 3, 2, 2, 2, 2, 2, 2, 2, 3, 2, 3, 1, 3, 2, 1, 3, 3, 3, 3, 3, 3, 3, 2, 3, 2, 3, 3, 2, 3, 2, 3, 2, 3, 3, 1, 2, 1, 1, 3, 3, 3, 2, 3, 3, 3, 3, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 2, 2, 3, 2), y=c(1.43, 1.45, 1.5, 1.56, 1.57, 2.9, 2.93, 2.94, 3.01, 3.06, 3.08, 3.14, 4.17, 4.23, 4.28, 4.29, 4.6, 1.73, 1.95, 1.97, 2.08, 2.1, 2.13, 2.2, 2.27, 2.28, 2.38, 2.4, 2.42, 2.43, 2.45, 2.53, 2.87, 2.88, 6.46, 6.48, 7.04, 8.94, 9.32, 11.7, 1.58, 1.6, 1.7, 1.77, 1.8, 1.82, 1.84, 1.88, 1.92, 1.95, 2.0, 2.0, 2.03, 2.08, 2.13, 2.18, 2.19, 3.76, 3.83, 3.94, 4.12, 4.17, 7.91, 2.24, 2.5, 2.76, 2.76))
mydata2 <- mydata %>% mutate( grp = str_c(yclass, "_", d1) )
mydata2_grp <- groupedData( y ~ d2 + d5 + d9 | grp, data = mydata2)
# Linear mixed model fitting
m1 <- lme(y ~ d2 + d5 + d9, random = ~ 1 | yclass/d1, correlation = corAR1(), method = "REML", data = mydata2_grp)
# Nonlinear mixed model fitting
m2 <- nlsList(y ~ SSlogis(d2 + d5 + d9, Asym, xmid, scal), data = mydata2_grp)
m3 <- nlme(y ~ NLS.YL(d2 + d5 + d9, i, A), data=mydata2_grp,
fixed = list(i ~ 1, A ~ 1),
random = pdSymm(list(i ~ 1, A ~ 1)))
How do I get the above model to run? I am also open to any other suggestions about the most effective way to model my data. Thank you for your time and help!
