nlme repeated measures
On 13-11-25 09:32 PM, Elizabeth Webb wrote:
Hi All-
I am trying to fit a light response curve (model photosynthesis based on
photosynthetically active radiation) with measurements taken on the same
plots over many days. I can get the model to converge in nls but to take
into account the repeated measures, I need to use nlme. My code is here
(where the starting values were taken from the nls output):
ww.light_grouped<-groupedData(gC~par|block/fence/plot,data=ww.light)
mixedWW<-nlme(model=gC~(alpha*par*Fmax)/(alpha*par+Fmax)+Rd,fixed=alpha+Fmax~1,start=c(alpha=-0.0001795,Fmax=-0.0092121),
data=ww.light_grouped,random=Rd~1|block/fence/plot)
When I run this, I get the following error:
Error in solve.default(estimates[dimE[1] - (p:1), dimE[2] - (p:1), drop =
FALSE]) :
system is computationally singular: reciprocal condition number =
3.55879e-37
I know that I cannot be the first person to try to fit a repeated measures
light response curve (most of the literature seems to be on coding in SAS
though). Any ideas on how to solve this error or other ways to fit this
type of curve in R?
Thanks,
Elizabeth
Your main problem is that you also need 'Rd' in the model as a fixed
effect -- the random effects are defined to have mean zero. (Note that
these data are sufficiently noisy that the block variance = residual
variance/100 (i.e. ratio of 10 in the standard deviations) and that the
other two variance terms are essentially zero ...)
ww.light_grouped<-groupedData(gC~par|block/fence/plot,data=ww.light)
nls0 <- nls(gC~alpha*par*Fmax/(alpha*par+Fmax)+Rd,
start=c(alpha=-0.0001795,Fmax=-0.0092121,Rd=0),
data=ww.light)
pframe <- data.frame(par=1:700)
pframe$gC <- predict(nls0,newdata=pframe)
ww.light <- transform(ww.light,
bfp=interaction(block,fence,plot))
library(ggplot2)
ggplot(ww.light,aes(x=par,y=gC)) +geom_point(aes(colour=bfp))+
geom_line(aes(group=bfp,
colour=bfp))+
geom_line(data=pframe,lwd=2,alpha=0.4)
library(mgcv)
ggplot(ww.light,aes(x=par,y=gC)) +geom_point(aes(colour=bfp))+
geom_smooth(aes(group=interaction(block,fence,plot),
colour=interaction(block,fence,plot)),
se=FALSE,method="gam")+
geom_line(data=pframe,lwd=2,alpha=0.4)
ggplot(ww.light,aes(x=par,y=gC)) +geom_point(aes(colour=bfp))+
geom_smooth(aes(group=bfp,
colour=bfp),se=FALSE)+
coord_cartesian(ylim=c(-0.03,0.04))
mixedWW<-nlme(model=gC~(alpha*par*Fmax)/(alpha*par+Fmax)+Rd,
fixed=alpha+Fmax+Rd~1,start=c(alpha=-0.0001795,Fmax=-0.0092121,
Rd=0.0139),
data=ww.light,random=Rd~1|bfp)
mixedWW2<-nlme(model=gC~(alpha*par*Fmax)/(alpha*par+Fmax)+Rd,
fixed=alpha+Fmax+Rd~1,start=c(alpha=-0.0001795,Fmax=-0.0092121,
Rd=0.0139),
data=ww.light,random=Rd~1|block/fence/plot)
P.S. My data:
dput(ww.light)
structure(list(par = c(212, 208, 503, 89, 61.9, 155, 49, 35,
641.5, 532.6, 345, 559, 80, 33, 21, 173, 750, 210, 151, 169,
334, 34, 17, 121, 689.5, 689.5, 140, 152, 116, 142, 183, 366,
24, 68, 119, 626, 603.9, 364, 194, 118, 185, 149, 45, 383, 83,
81, 141, 308.6, 300, 297, 500, 27, 25, 315, 97, 76, 165, 487.7,
453.2), gC = c(0.00823490377389, 0.001915093741032, -0.022289772410622,
0.0227122610049, 0.02557586905398, 0.02456920537407, 0.002374206642954,
0.013696390897218, 0.040337047839618, 0.009981229463982, -0.00674478119502,
-0.001055959594434, 0.030835240675164, 0.006106703096412, 0.01421021063421,
-0.000526876101036, 0.00202679213112, -0.000838672558272, 0.00509918585139,
0.006495516099864, 0.001626921909792, 0.006247252588266, 0.004622651392392,
0.00063338751144, 0.002997109554102, 0.010496218151088, 0.002770538321952,
-0.003660578166654, -0.00229353987582, 0.004299458228208,
0.003570011155386,
-0.004095960060402, 0.00200178972513, 0.007373300376096, 0.001331055108216,
0.009435928103748, 0.026637682283982, -0.006434379185256,
-0.001959871385466,
0.010275562922436, -0.000594160634394, 0.015110760372264, 0.01182487186701,
0.004280099006556, 0.00569932660827, 0.00988832461734, 0.002471339105676,
0.020481226398234, 0.004056278470344, -0.00501195668094, -0.0063009487386,
0.00757298331279, 0.021282565362282, 0.006687425246994, 0.00579015733752,
0.010894561817346, 0.007070991129486, 0.020669872048848, 0.005616670768956
), block = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L,
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L,
3L), .Label = c("A", "B", "C"), class = "factor"), fence = c(1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 3L, 3L,
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 4L, 4L, 4L,
4L, 4L, 4L, 4L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 6L,
6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L), plot = structure(c(6L,
8L, 8L, 6L, 8L, 8L, 8L, 8L, 6L, 8L, 8L, 8L, 8L, 8L, 8L, 6L, 8L,
8L, 6L, 8L, 8L, 8L, 8L, 8L, 6L, 8L, 6L, 8L, 8L, 6L, 8L, 8L, 8L,
8L, 8L, 6L, 8L, 6L, 8L, 6L, 8L, 6L, 8L, 8L, 8L, 8L, 8L, 6L, 6L,
8L, 8L, 6L, 8L, 8L, 8L, 8L, 8L, 6L, 8L), .Label = c("1", "2",
"3", "4", "5", "6", "7", "8", "A", "B", "C", "D"), class = "factor")),
.Names = c("par",
"gC", "block", "fence", "plot"), row.names = c(630L, 632L, 634L,
638L, 640L, 642L, 646L, 648L, 649L, 651L, 654L, 656L, 659L, 663L,
665L, 669L, 671L, 673L, 675L, 677L, 679L, 681L, 683L, 685L, 687L,
689L, 691L, 693L, 695L, 697L, 699L, 701L, 703L, 705L, 707L, 709L,
711L, 713L, 715L, 717L, 719L, 721L, 723L, 725L, 727L, 729L, 731L,
733L, 737L, 739L, 741L, 743L, 745L, 747L, 749L, 751L, 753L, 755L,
757L), class = "data.frame")
[[alternative HTML version deleted]]
_______________________________________________ R-sig-mixed-models at r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models