I have only a very basic background in statistics, and I have a possibly simple question, but I'm having a bit of trouble with my model. I suppose this is also an R question, but also statistical!

I'm looking for some advice concerning the interpretation my nlme results.

I am modelling photosynthesis-irradiance -relationship and have fitted a nonlinear mixed model in nlme such as this with 3 parameters to estimate:

nlme(Foto~picurve(PAR, Amax, Aqe, LCP), fixed = list(Amax ~ Place, Aqe + LCP ~ 1), random = pdDiag(Amax ~ 1), weight = varPower(), correlation = corARMA(p = 1, q = 1), method = "ML",  start = rep(c(15, 0.0054, 20), 2), data = data1)

The summary is as follows


    Fixed effects: list(Amax ~ Place, Aqe + LCP ~ 1) 
                                           Value Std.Error  DF  t-value p-value
    Amax.(Intercept)                   16.427779 0.6567488 959 25.01380  0.0000
    Amax.Place 2                       -1.328056 0.9169505 959 -1.44834  0.1478
    Amax.Place 3                       -1.063690 0.8996467 959 -1.18234  0.2374
    Amax.Place 4                       -3.207345 0.9171032 959 -3.49726  0.0005
    Aqe                                 0.057579 0.0015047 959 38.26518  0.0000
    LCP                                21.388703 0.7486608 959 28.56928  0.0000


I need to make inference on the fixed effects at the 4 places. Each place contains 3 genotypes, and each genotype contain several individuals, but they're not modelled here. Some genotypes contain more individuals than others (unbalanced?), and some individuals may have a few more measurements in them than the others. But I suppose that's beyond the point.

Now my questions is, what is the current recommended way of making inferences about the fixed effects parameter estimates? If I understood correctly, MCMCs etc. are the "most correct" approach, but they are probably beyond my comprehension at the moment and probably not what my superiors would want anyway. I know I can obtain the conditional F-test statistics and compare the places with anova(), even if that probably is a bad approach due to the denominator df problems I've been hearing about. But other than that, I don't know any way of doing this (recommendations are welcome).

Next, and more importantly, how does one make pairwise/multiple comparisons between the fixed effects? I need to compare all places to each others, as so:

    Place 1 - Place 2
    Place 1 - Place 3
    Place 1 - Place 4
    Place 2 - Place 3


Now, again, MCMCs might not be my cup of tea here. The conditional t-test values from the summary can probably be used, at least in some (balanced?) instances, but would still leave me wanting to compare other groups besides the first.

What I would have to do, if I am correct (am I?), is set my contrasts so that I can compare the other groups as well. I tried approaching this with the multcomp-package:

summary(glht(psn7, linfct =  matrix(c(0, 0, 1, -1, 0, 0), 1)))

Linear Hypotheses:
       Estimate Std. Error z value Pr(>|z|)  
1 == 0   2.1437     0.8862   2.419   0.0156 *
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Adjusted p values reported -- single-step method)

This is where the troubles begin.

1) Is glht unsuited for nlme? It reports z-values (df problems?) - I thought it always reported t-values...

2) I still cannot compare the first group and the others, do I need to manipulate the order of my data or is there a more elegant way of getting all the comparisons?

I know this a very simple question, and I apologize. I'm quite new to this all...

p.s. also, bonus question, how about likelihood ratio tests? Can they be used for inference about fixed effects?


You can get the needed comparisons as follows:

con <- matrix(c(0, 1, 0, 0, 0, 0,
                0, 0, 1, 0, 0, 0,
                0, 0, 0, 1, 0, 0,
                0,-1, 1, 0, 0, 0,
                0,-1, 0, 1, 0, 0,
                0, 0,-1, 1, 0, 0),
               nrow = 6, byrow = TRUE)
dimnames(con)[[1]] <- c("2-1", "3-1", "4-1", "3-2", "4-2", "4-3")
glht(your.model, con)

The reason is because the way it is parameterized, the intercept for Amax estimates the value for the omitted indicator, which is Place 1, and the other coefficients estimates differences between their respective places and Place 1.

It is important to include all the comparisons in one matrix so that glht performs the correct adjustment for multiple comparisons.

It is not straightforward to get degrees of freedom from nlme objects, so glht instead outputs the asymptotic statistics ($z$ rather than $t$). It appears you have a lot of data so this won't make much difference.

I hope this is still useful, even though it's been several days...

| cite | improve this answer | |
  • $\begingroup$ That's excellent! However, during the few days that have passed, I actually managed to do this with the estimable() -function from package gmodels, although I did it only a step at a time like this: estimable(model, cm=c(0,-1,1,0,0,0)) The estimates would be identical to your suggestion, but estimables would also report std. errors and t-statistics. From your code using glht I would only be able to get the estimates. Is there a way to get the test statistics as well, and will they indeed be different then? $\endgroup$ – tuhinokkaeläin Feb 19 '15 at 22:31
  • $\begingroup$ Yes - it's subtle, but if you use summary on the results of the glht call, you get that stuff! $\endgroup$ – Russ Lenth Feb 19 '15 at 22:43
  • $\begingroup$ Well I guess I should have seen that one coming (: Your answer was very useful, I would have ended up doing it wrong otherwise, without imputing the whole matrix at once, and using a test that estimates the df. I'm seeing some big changes in p-values now, but like you said, not much difference in the estimates. Indeed, apart from trusting a z- or a t-statistic, and p-values therein, could I still ask you to help point me the way to how inference in NLME-models could be done otherwise? I would be eager to learn a better way. Some entry-level references are welcome, if such exist. $\endgroup$ – tuhinokkaeläin Feb 20 '15 at 0:04

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.