Posting here instead of SO because I guess this is still more of a statistical question in the sense of "am I doing it right". Also includes help needed with R.
I'm a complete newbie to nonlinear mixed-effects models and pretty new to R (and statistics in general).
What I tried
I've been googling it for days, also I've been reading Zuur et al, Pinheiro & Bates etc. but haven't so far figured it out.
I'm trying to fit a nonlinear photosynthesis-irradiance -curve to
data1 using nlme.
Response is photosynthesis (
Photo), independent variable is irradiance (
Photo is measured at several levels of
PAR, from 4
Groups, each containing 3
data1 is a longitudinal repeated measures file arranged as
groupedData as follows:
head(data1) Grouped Data: Photo ~ PAR | Groups/Genotypes Photo PAR Groups Genotypes 1 -0.56 1 BoK K1 2 0.41 25 BoK K1 3 1.07 50 Bok K1 ...
N.B.! Other groups are coded BoL, KoK and KoL. Within each group lie 3 genotypes, BUT:
BoK contains genotypes K1, K2, K3.
BoL contains genotypes L1, L2, L3.
KoK contains genotypes K1, K2, K3.
KoL contains genotypes L1, L2, L3.
i.e. BoK and KoK contain same genotypes, coded similarly in the table, and BoL and KoL contain same genotypes, also coded with the same designation.
The nonlinear model I have adopted to fit the data contains three parameters (Amax, Aqe, LCP):
Perform curve fitting using NLME, use model selection to find parameters to include in fixed and random effects and obtain parameter estimates for them, and thus make inference about possible differences between groups, while at the same time considering variation within the groups caused by genotypes.
Code so far
# This code is probably full of dirtyness and unnecessities Data1 <- read.csv2(file="data.csv", sep=";") attach(Data1) library(nlme) data1=groupedData(Photo~PAR|Groups/Genotypes) # Model picurve = function(x, Amax, Aqe, LCP) (Amax*(1-exp(-Aqe*((x-LCP))/Amax))) # Model fitting started by treating all parameters as random and without considering the effects of Groups psn1 = nlme(Photo~picurve(PAR, Amax, Aqe, LCP), fixed = (Amax + Aqe + LCP ~ 1), random = pdDiag(Amax + Aqe + LCP ~ 1), start = c(15, 0.0054, 20), data = data1) # The model selection process so far has lead me here psn5 = nlme(Photo~picurve(PAR, Amax, Aqe, LCP), fixed = (Amax + Aqe + LCP ~ 1), random = pdDiag(Amax ~ 1), start = c(15, 0.0054, 20), data = data1) # So far so good. Now if I suppose that I have determined the random structure by only leaving in Amax, I would next like to include the effects of Groups. But when I try to update my model: psn51 = update(psn5, fixed = (Amax + Aqe + LCP ~ Groups), start = rep(c(15,0.0054,20), 4)) # I always get the error: Error in MEEM(object, conLin, control$niterEM) : Singularity in backsolve at level 0, block 1
Now, I have several questions.
I have tried googling the error message, and the first thing I fear, and must ask, is it simply because of my data structure? Is it just that I'm not even doing this correctly? If it's not in the data structure itself, is it something that I did when I defined the
If it's not that, then could it be that the random effects structure doesn't work as it stands? Could that have something to do with it? The estimated variance for the parameter Amax (only one I left in) in the random effects was tiny as well, so should I actually just leave everything out of the random structure and change my approach to nls or something else at that point?
Am I even modeling things right by writing
fixed = (Amax + Aqe + LCP ~ Groups)?
What if I wanted to exclude a parameter, so that, say,
Groups affected Aqe and LCP, but not Amax, I thought I should use
fixed = (Amax ~ 1, Aqe + LCP ~ Groups) but this returns a
Error: unexpected ',' in "psn51 = update(psn5, fixed = (Amax ~ 1,". Surely it can't be right to use
fixed = (Amax ~ 1 + Aqe + LCP ~ Groups)?
(Obviously I have a lot of trouble even understanding how to use these formulae)
As an additional trouble I have a hard time figuring out how to input starting values for the iterations. When I update the fixed parameters that are affected by
Groups, the list of starting values have to be updated as well, but I seem to miss the logic behind this list? Can somebody shed some light to this seemingly simple question? Something to do with a model matrix?
As an ending note, I'm sorry to even ask these beginner questions here, it may be that I should not even be using NLME for it is something I do not fully understand, but the data I have seems to almost require it. I would be super happy to get even a little help here!