[R-sig-ME] nlme repeated measures

[Ben Bolker]

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
>