1
$\begingroup$

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.

Background

I'm trying to fit a nonlinear photosynthesis-irradiance -curve to data1 using nlme.

Response is photosynthesis (Photo), independent variable is irradiance (PAR).

Photo is measured at several levels of PAR, from 4 Groups, each containing 3 Genotypes.

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 model

The nonlinear model I have adopted to fit the data contains three parameters (Amax, Aqe, LCP):

(Amax*(1-exp(-Aqe*((x-LCP))/Amax)))

Research goals

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.

1)

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 groupedData?

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)?

2)

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)

3)

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!

$\endgroup$
  • $\begingroup$ If you could give us ful access to your data maybe somebody wil give it a try ... $\endgroup$ – kjetil b halvorsen Jan 26 '15 at 11:49

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.