I have gene expression data of one gene at 6 different time-points from 3 samples (biological replicates). I am trying to determine if the expression follows a 24h cycle using a non-linear regression model. I have successfully been able to fit the data to a harmonic curve using the mean of the expression with NLS as follows:
# Expression data, the first 3 are replicates of time-point 1,
# the next 3 are replicates of time-point 5 and so on.
data <- c(21.974280, 28.056115, 26.143261, 5.471526, 6.095989, 7.513615
, 12.386668, 6.849121, 7.730032, 27.745546, 26.861919, 24.709190
, 52.788829, 40.106706, 44.755854, 45.154904, 54.638577, 43.061662)
# Time-points for each measurement
groups <-c(1,1,1,5,5,5,9,9,9,13,13,13,17,17,17,21,21,21)
# Mean of expression in each time-point
y <- c(25.391219, 6.360377, 8.988607, 26.438885, 45.883796, 47.618381)
# time-points
t <- c(1, 5, 9, 13, 17, 21)
# Model
fit <- nls(y ~ amp*sin((2*pi/per*t) + phase) + C
, start = list(amp = diff(range(y))/2
, per = 24
, phase = 0
, C = mean(y))
, control = list(maxiter = 500, warnOnly = TRUE)
)
I used this as a first approximation to harmonic curves and regression, just to make sure I understood how everything works. However now I'm trying to incorporate info about the replicates as random effects for each time-point; that is, use all the information instead of calculating the mean. NLS can't handle that, so I've been advised to try NLME instead. I don't quite understand how to incorporate the random and fixed effects, though. I've looked for examples or tutorials on several forums and web pages unsuccessfully, and the R package documentation does not contain an example I can follow. The closest thing to what I want is this post, but it uses LME which doesn't seem to have the same syntax as NLME. From that post I've come up with this attempt:
samples <- paste("sample", c(1:18), sep="_")
nlme( y ~ amp*sin((2*pi/per*t)+phase)+C
, start = list(amp = diff(range(y))/2
, per = 24
, phase = 0, C = 0)
, random = ~ 1 | sample/groups
, data = data)
Which throws the following error:
Error in nlme.formula(y ~ amp * sin((2 * pi/per * t) + phase) + C, start = list(amp = diff(range(y))/2, : argument "fixed" is missing, with no default
because I'm not sure how to introduce the random and fixed effects of my model, or if indeed I am specifying correctly the random effect of individuals.
Can anyone please explain how to use NLME in this particular case?
