[R] nlme starting values are not the correct length
Lara Reichmann
Lara.Reichmann at austin.utexas.edu
Thu Dec 6 16:45:41 CET 2012
Hi Petr,
Thanks for the email, the number of parameters is not 3 because I am trying to fit a parameter for each treatment level, therefore the start vector changes depending on the fixed statement (see below).
I finally got some of these to work....but for some reason I cannot run a complete model with both fixed effects and interactions (Species*Fert), I still couldn't get the syntax for the start vector!
Thanks again,
Lara
All these work:
<unfert <- subset(lightresponse, Fert==0)
<model.cov4<-nlme(Photo~A*(1-exp(-lC*PARi/A))-B,fixed=A+B+lC~Species,random=A+B+lC~1,data=unfert,start=c(24.89,0,0,0,1.78,0,0,0,0.06,0,0,0))
<plot(augPred(model.cov4,outer=T))
#Here is the summary of the non-linear mixed model:
<summary(model.cov4)
<fert <- subset(lightresponse, Fert==1)
<model.cov5<-nlme(Photo~A*(1-exp(-lC*PARi/A))-B,fixed=A+B+lC~Species,random=A+B+lC~1,data=fert,start=c(24.89,0,0,0,1.78,0,0,0,0.06,0,0,0))
<plot(augPred(model.cov5))
#Here is the summary of the non-linear mixed model:
<summary(model.cov5)
<model.cov6<-nlme(Photo~A*(1-exp(-lC*PARi/A))-B,fixed=A+B+lC~Fert,random=A+B+lC~1,data=lightresponse,start=c(24.89,0,1.78,0,0.06,0))
<plot(augPred(model.cov6))
#Here is the summary of the non-linear mixed model:
<summary(model.cov6)
This one I am interested in but cannot figure out the start vector!
<model.nlme<-nlme(Photo~A*(1-exp(-lC*PARi/A))-B,fixed=A+B+lC~Species*Fert,random=A+B+lC~1,data=lightresponse,start=c(?????)
Also, once I have the summary of effects, how do I do contrasts for specific species differences?
I understand from this output that the p value is comparing each species with the baseline species, but how do I compare jg with sw for example? Is there a simple "Tukey" or family contrast for nlme method?
Fixed effects: A + B + lC ~ Species
Value Std.Error DF t-value p-value
A.(Intercept) 24.132181 3.649618 281 6.612249 0.0000
A.Speciesjg -4.233378 5.160757 281 -0.820302 0.4127
A.Specieslb 2.862573 5.161393 281 0.554612 0.5796
A.Speciessw 4.964391 5.344752 281 0.928835 0.3538
B.(Intercept) 1.523276 0.262774 281 5.796899 0.0000
B.Speciesjg 0.115693 0.372653 281 0.310457 0.7564
B.Specieslb 0.469510 0.372018 281 1.262062 0.2080
B.Speciessw 0.476496 0.380826 281 1.251220 0.2119
lC.(Intercept) 0.065695 0.005050 281 13.009079 0.0000
lC.Speciesjg -0.009383 0.007155 281 -1.311311 0.1908
lC.Specieslb 0.009004 0.007178 281 1.254371 0.2107
lC.Speciessw -0.002268 0.007369 281 -0.307801 0.7585
On Dec 6, 2012, at 12:52 AM, PIKAL Petr <petr.pikal at precheza.cz> wrote:
> Hi
>
> see inline
>
>> -----Original Message-----
>> From: r-help-bounces at r-project.org [mailto:r-help-bounces at r-
>> project.org] On Behalf Of Lara Reichmann
>> Sent: Wednesday, December 05, 2012 1:53 AM
>> To: R-help at lists.R-project.org
>> Subject: [R] nlme starting values are not the correct length
>>
>> 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,d
>>> ata=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
>>
>
> A aqm not an expert in nlme models but it seems that you have only three fixed parameters A,B,C and want to feed them with 12 starting values. This is probably the reason for error.
>
> Maybe you shall also search in R-sig-mixed-models
>
> Regards
> Petr
>
>
>> 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 at 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.
More information about the R-help
mailing list