I am getting an error message in a call to nlme and cannot understand what
is happening. I explain the steps below in the hope that someone can
explain the error and how to correct it.
STEP 1: Data set: name: marouane.data. This is a data frame whose first few
lines are as follows:
> marouane.data[1:13,]
species plant leaf irradiance photosynthesis chlorophyll
1 Asclepias.incarnata 1 1 0 -2.091359 0.02619
2 Asclepias.incarnata 1 1 50 1.153241 0.02619
3 Asclepias.incarnata 1 1 100 2.241963 0.02619
4 Asclepias.incarnata 1 1 200 3.541258 0.02619
5 Asclepias.incarnata 1 1 300 5.012547 0.02619
6 Asclepias.incarnata 1 1 400 5.710689 0.02619
7 Asclepias.incarnata 1 1 600 6.632571 0.02619
8 Asclepias.incarnata 1 1 800 7.314284 0.02619
9 Asclepias.incarnata 1 1 1000 8.092402 0.02619
10 Asclepias.incarnata 1 1 1310 8.110368 0.02619
11 Asclepias.incarnata 1 2 0 -2.051934 0.02707
12 Asclepias.incarnata 1 2 50 1.213854 0.02707
13 Asclepias.incarnata 1 2 100 2.271453 0.02707
absorbance nitrogen sla
1 0.8816 18.35913 77.53
2 0.8816 18.35913 77.53
3 0.8816 18.35913 77.53
4 0.8816 18.35913 77.53
5 0.8816 18.35913 77.53
6 0.8816 18.35913 77.53
7 0.8816 18.35913 77.53
8 0.8816 18.35913 77.53
9 0.8816 18.35913 77.53
10 0.8816 18.35913 77.53
11 0.8813 18.35913 77.53
12 0.8813 18.35913 77.53
13 0.8813 18.35913 77.53
The nesting structure is species (25), plants within species (4 each
species) and leaves within plants (2 each plant). There is only 1 missing
value in the data set.
STEP 2: I constructed a self starting function (photo) that gives the
nonlinear function and provides initial estimates of the three parameters.
This self starting function works. I use this to call the nlsList function,
which fits separate nonlinear functions to each species (I am only working
on a 2 level model so far: between and within species):
fit<-nlsList(photosynthesis~photo(irradiance,Am,Q,LCP)|species,
+ data=marouane.data,na.action=na.omit)
This works and I get the three parameter estimates per species.
STEP 3: use nlsList to call nlme. I use the nlsList object (fit) to fit a
variance components model:
fit.nlme<-nlme(fit)
This works (at least it runs without an error message). To see what syntax
is used, I extract the call:
> fit.nlme$call
nlme.formula(model = photosynthesis ~ photo(irradiance, Am, Q,
LCP), data = marouane.data, fixed = list(Am ~ 1, Q ~ 1, LCP ~
1), random = list(species = c(2.41887429868478, -1.53351210977838,
2.47282942435836, 0, 0, 0)), groups = ~species, start = list(
fixed = c(11.7384877502532, 0.109091641284836, 9.91905527125462
)), na.action = na.omit)
Question 1: the syntax to the random argument seems wrong! It should be
something like: random=(list(Am~1,Q~1,LCP~1)). I have no idea where the 6
numbers (2.41887...) are coming from in the random argument. Can someone
explain how the nlme function obtains such an argument for random? The
values in fixed are the estimated fixed effects from nlsList.
If I follow the instructions in Pinheiro & Bates, the call should be (with a
self-starting function):
nlme(model = photosynthesis ~ photo(irradiance, Am, Q,
LCP), data = marouane.data, fixed = list(Am ~ 1, Q ~ 1, LCP ~
1), random = list(Am~1,Q~1,LCP~1),groups=~species,
na.action=na.omit).
When I use this, I get an error message:
> nlme(model=photosynthesis~photo(irradiance,Am,Q,LCP),data=marouane.data,
+ fixed=list(Am~1,Q~1,LCP~1),groups=~species,
+ random=list(Am~1,Q~1,LCP~1),
+ na.action = na.omit)
Error in nlmeCall[[i]] <- NULL : subscript out of bounds
Bill Shipley
[[alternative HTML version deleted]]
