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?