Dear R community,
I am trying to fit an nlme model where I want to estimate the fixed effects of
two treatments on the parameters on the following equation
Photo~(a*(1-exp(-c*PARi/a)))-b
I was able to fit a simple model without covariates following the method
described in Mixed-Effects Methods and Classes for S and S-PLUS, version 3.0,
but when I add the covariates, I get the error " starting values for the fixed
component are not the correct length"
My data has the following structure "Subject" "Species" "Fert" "Photo" "PARi" ,
where several "Photo" measurements where taken on the same subject by changing
"PARi", 4 Species levels and 2 Fert levels, there are 31 Subjects (one missing
value), and 323 observations
DATA extract
Subject Species Fert Photo PARi
1 bb 1 22.5389 1499.3307
1 bb 1 21.881 1248.913
1 bb 1 21.2862 999.3387
1 bb 1 20.5836 799.9308
1 bb 1 19.3758 601.1412
1 bb 1 15.5915 399.815
1 bb 1 8.7978 200.1087
1 bb 1 4.4347 99.686
1 bb 1 2.0387 49.7842
1 bb 1 -1.4854 0.0576
2 sw 0 6.782 1500.5337
2 sw 0 7.1432 1249.2749
2 sw 0 7.3319 1000.9891
2 sw 0 7.5848 799.1752
2 sw 0 7.1882 599.5544
2 sw 0 6.809 399.988
2 sw 0 5.3877 198.7574
2 sw 0 3.5104 100.7115
2 sw 0 0.8856 50.7015
2 sw 0 -1.121 0.0569
3 jg 1 16.0827 2000.4941
3 jg 1 16.0236 1501.1957
3 jg 1 16.3818 1248.9551
3 jg 1 16.7815 1499.6414
3 jg 1 17.175 2000.6851
3 jg 1 16.6529 1000.2707
3 jg 1 15.7987 799.676
3 jg 1 15.5437 598.9409
3 jg 1 11.7683 400.7715
3 jg 1 4.89 200.7468
3 jg 1 4.1294 100.9664
3 jg 1 1.6008 50.9254
3 jg 1 -0.89 0.5347
4 sw 1 25.2889 2000.1454
4 sw 1 24.7284 1499.6191
4 sw 1 24.3637 1249.7523
4 sw 1 23.3523 1000.0944
4 sw 1 21.6057 800.2209
4 sw 1 18.8926 599.7022
4 sw 1 14.6598 398.9366
4 sw 1 7.7182 201.5697
4 sw 1 3.4775 100.5139
4 sw 1 1.169 49.7045
4 sw 1 -1.3558 1.6914
5 jg 0 6.1626 2000.9351
5 jg 0 7.5573 1499.6581
5 jg 0 7.7129 1249.5073
5 jg 0 7.442 1000.7276
5 jg 0 7.5135 799.1286
5 jg 0 7.1559 599.5568
5 jg 0 6.8161 400.3576
5 jg 0 4.0097 199.7442
5 jg 0 2.7202 101.1253
5 jg 0 1.0746 51.1787
5 jg 0 -0.5913 0.975
This works so far:
lightresponse<-groupedData(Photo~PARi|Subject,data=lightr,outer = ~ Species *
Fert,labels = list(x = "PAR", y = "CO2 uptake rate"),units = list(x = "(photon
s-1)", y = "(umol/mˆ2 s)"))
Photo.resp<-function(PARi,A,B,C)A*(1-exp(-C*PARi/A))-B
Photo.resp<-deriv ((~A
*(1-exp(-C*PARi/A))-B),c("A","B","C"),function(PARi,A,B,C){})
> lightresp.fit1<-nlme(model=Photo~Photo.resp(PARi,A,B,C),fixed=A+B+C~1,data=lightresponse,start=c(30,-5,0.1))#fitting
> nlme without any covariates
>lightresp.fit1
OUTPUT
> lightresp.fit1
Nonlinear mixed-effects model fit by maximum likelihood
Model: Photo ~ Photo.resp(PARi, A, B, C)
Data: lightresponse
Log-likelihood: -494.5926
Fixed: A + B + C ~ 1
A B C
24.89334793 1.77983637 0.06499634
Random effects:
Formula: list(A ~ 1, B ~ 1, C ~ 1)
Level: Subject2
Structure: General positive-definite, Log-Cholesky parametrization
StdDev Corr
A 10.67382785 A B
B 0.52572012 1.000
C 0.01433605 0.371 0.384
Residual 0.71900020
Number of Observations: 323
Number of Groups: 31
##Now, I want to test the effect of Species and Fert, I don't fully understand
how to modify the start vector, as I tried several options and no one seems to
be correct. Do the number of levels in each factor matter? In that case 4
Species and 2 Fert levels, I would need 6 initial parameters x 3? This didn't
work either
>lightresp.fit2<-nlme(model=Photo~Photo.resp(PARi,A,B,C),fixed=A+B+C ~
>Species*Fert,random=A+B+C~1,data=lightresponse,
>start=c(24.89,0,0,0,1.78,0,0,0,0.065,0,0,0))
Error in nlme.formula(model = Photo ~ Photo.resp(PARi, A, B, C), fixed = A + :
starting values for the fixed component are not the correct length
I hope someone out there has the answer!
Thanks!!!
Lara
Lara G. Reichmann
Postdoctoral Fellow
USDA-ARS
808 E Blackland Rd
Temple, TX 76502
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
