[R] nlme starting values are not the correct length

Lara Reichmann Lara.Reichmann at austin.utexas.edu
Wed Dec 5 01:52:39 CET 2012


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



More information about the R-help mailing list