The crux of my problem (as far as I can tell) is that no examples are out there that give the proper syntax for nested groups when coding the groups argument of nlme. All examples only give a single group as an argument, if groups is even included in the function call at all.
I'm trying to fit a Non-Linear Mixed Effect model to ELISA plate data, where each plate contains sigmoidal curves for two different samples ("Ref" and "Test"). The formula is a 3-parameter logistic curve with Parameters A, B and C, and is shown here: Nonlinear mixed effect formula
"i" refers to the sample preparation, "j" refers to the plate, and "k" refers to the dilution factor (basically the "x" to the assay signal's "y"). Fixed effects are the capital letters (A, B, C) and the random effects based on plate are the lower-case ones (a, b, c).
Here is the code block to define the function (a simplified version that I understand to be equivalent to the image above), and the nlme function call that is currently producing an error. (The data argument has a data frame with columns for preparation, plate, dilution and plate signal ("y").)
# Define the 3-parameter logistic function
logistic3PL <- function(k, A, B, C, a, b, c) {
(A + a) / (1 + (k / (C + c))^ (B + b))
}
# Define the nonlinear mixed effects model
model <- nlme(y ~ logistic3PL(k, A, B, C, a, b, c),
data = framedData,
fixed = A + B + C ~ 1,
random = a + b + c ~ 1,
groups = ~ i/j,
start = c(A = 1, B = 1, C = 1, a = 0, b = 0, c = 0))
The groups argument appears to be the issue. I believe that j (the plate) should be nested within i (i.e. a collection of plate data respective to one test sample type).
The line groups = ~ i/j, as written above appears to follow the syntax from the nlme documentation. However, it produces an error reading:
Error in names(reSt) <- namGrp :
'names' attribute [2] must be the same length as the vector [1]
I've also tried groups = ~ i + j, emulating the fixed and random argument, but that seems to be even further from the mark:
Error in getGroups.data.frame(dataMix, eval(parse(text = paste("~1", deparse(groups[[2]]), :
invalid formula for groups
Any assistance will be much appreciated. Thank you in advance!
6/21/23 10:30am PST, edited to add sample data. This is two plates' worth of data - two dose-response curves for each preparation, so hopefully that should be enough to model random effects.
structure(list(i = c("Ref", "Ref", "Ref", "Ref", "Ref", "Ref",
"Ref", "Ref", "Ref", "Ref", "Ref", "Test", "Test", "Test", "Test",
"Test", "Test", "Test", "Test", "Test", "Test", "Test", "Ref",
"Ref", "Ref", "Ref", "Ref", "Ref", "Ref", "Ref", "Ref", "Ref",
"Ref", "Test", "Test", "Test", "Test", "Test", "Test", "Test",
"Test", "Test", "Test", "Test"), j = c("1", "1", "1", "1", "1",
"1", "1", "1", "1", "1", "1", "1", "1", "1", "1", "1", "1", "1",
"1", "1", "1", "1", "2", "2", "2", "2", "2", "2", "2", "2", "2",
"2", "2", "2", "2", "2", "2", "2", "2", "2", "2", "2", "2", "2"
), k = c("5", "10", "20", "40", "80", "160", "320", "640", "1280",
"2560", "5120", "5", "10", "20", "40", "80", "160", "320", "640",
"1280", "2560", "5120", "5", "10", "20", "40", "80", "160", "320",
"640", "1280", "2560", "5120", "5", "10", "20", "40", "80", "160",
"320", "640", "1280", "2560", "5120"), y = c("0.0128", "0.0269",
"0.0611", "0.1783", "0.4333", "0.6237", "0.846", "0.8952", "0.9145",
"0.9443", "0.9831", "0.00470000000000001", "0.0087", "0.0223",
"0.0557", "0.137", "0.3419", "0.555", "0.6829", "0.7723", "0.8426",
"0.8433", "0.0332", "0.042", "0.0755", "0.172", "0.4082", "0.6905",
"0.8783", "0.9687", "1.0235", "1.0798", "1.0744", "0.0114", "0.0165",
"0.0368", "0.1016", "0.1442", "0.3734", "0.6292", "0.7851", "0.8644",
"0.8932", "0.9376")), row.names = c(NA, 44L), class = "data.frame")
