[R-sig-ME] HELP on fixed effect model building with nlmer or nlme

[Lara Reichmann]

Dear R community,

I've been trying for a week to fit an non-linear mixed effect model where I want to estimate the fixed effect 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 nlme method described in "The R Book" by Crawley, p669

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 subject), 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|Subject2,data=lightr,outer = ~ Species * Fert,labels = list(x = "PAR", y = "CO2 uptake rate"),units = list(x = "(photon s-1)", y = "(umol/mˆ2 s)"))  

model.ml1<-nlme(Photo~A*(1-exp(-C*PARi/A))-B,fixed=A+B+C~1,random=A+B+C~1|Subject2,data=lightresponse,start=c(A=24.89,B=1.78,C=0.02)) #from Crawley book, page 669
summary(model.ml1)
### model with no fixed effects

model.nlme<-nlme(Photo~A*(1-exp(-C*PARi/A))-B,fixed=A+B+C~Species,random=A+B+C~1|Subject2,data=lightresponse,start=c(24.89,0,0,0,1.78,0,0,0,0.06,0,0,0))
###model with Species as fixed effects, it works

                  Value Std.Error  DF   t-value p-value
A.(Intercept)  25.030590  1.857426 281 13.475959  0.0000
A.Species1     -0.898409  3.179607 281 -0.282554  0.7777
A.Species2     -5.131793  3.179135 281 -1.614211  0.1076
A.Species3      1.964169  3.179651 281  0.617731  0.5373
B.(Intercept)   1.788700  0.133278 281 13.420773  0.0000
B.Species1     -0.265424  0.228670 281 -1.160726  0.2467
B.Species2     -0.149737  0.229510 281 -0.652419  0.5147
B.Species3      0.204079  0.228995 281  0.891197  0.3736
C.(Intercept)  0.065034  0.002574 281 25.262323  0.0000
C.Species1     0.000662  0.004402 281  0.150301  0.8806
C.Species2    -0.008721  0.004413 281 -1.976211  0.0491
C.Species3     0.009666  0.004432 281  2.180979  0.0300

QUESTION 1: How do I know which species corresponds with the Intercept, Species1, 2, and 3? Values do not seem to match the means for "levels(Species)" [1] "bb" "jg" "lb" "sw"

##Now, I want to test the effect of Species* Fert, I don't fully understand how to modify the syntax and the start, 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

Question 2: How do I compare the significance of parameter estimates between Species and Fertilizer Treatments? 

I hope someone out there has the answer!
Thanks!!!

In other desperate attempt to run non linear models, I followed the example from useR!2011, University of Warwick, and I cannot even get the example to run!
 Th.start <- c(lKe = -2.5, lKa = 0.5, lCl = -3)
> nm1 <- nlmer(conc ~ SSfol(Dose, Time,lKe, lKa, lCl) ~
+              0+lKe+lKa+lCl+(0+lKe|Subject)+(0+lKa|Subject)
+              +(0+lCl|Subject), nAGQ=0, Theoph,
+              start = Th.start, verbose=2L)
Error in FUN(X[[1L]], ...) : no loop for break/next, jumping to top level

I cannot load the lme4a package, just the lme4, is this the problem?

Any help with the nlme or nlmer syntax would be awesome. I don't really know the difference, but any of these that can give me fixed factor effects on the parameters of the nonlinear model would work for me.


Lara