1
$\begingroup$

I am investigating the effects of leaf toughness, leaf's Nitrogen content and fire on the proportion of herbivory of plants.

My fire treatment is ordinal (fire severity): I have 75 plants located in unburned (control), burned once, and burned twice plots (25 in each fire treatment). Dataset: https://www.dropbox.com/s/kk4l0fmbwfyz821/df.csv?dl=0

I adjusted the following model:

data <- read.csv("df.csv", h=T)

#converting the 'treatment' variable to ordinal
data$treatment <- factor(data$treatment, order = TRUE, 
                                    levels = c("Control", "one_fire", "two_fires"))

class(data$treatment)
[1] "ordered" "factor"

#Adjusting betareg model
m1 <- betareg(prop.herbivory~treatment+nitrogen_content+LeafToughness, data=data)

And got the following result:

summary(m1)
Call:
betareg(formula = prop.herbivory ~ treatment + nitrogen_content + LeafToughness, 
    data = data)

Standardized weighted residuals 2:
    Min      1Q  Median      3Q     Max 
-1.9636 -0.6267 -0.0581  0.7116  2.3585 

Coefficients (mean model with logit link):
                  Estimate Std. Error z value Pr(>|z|)    
(Intercept)      -4.410491   0.438525 -10.058  < 2e-16 ***
treatment.L       0.362331   0.153044   2.367  0.01791 *  
treatment.Q      -0.043466   0.148272  -0.293  0.76941    
nitrogen_content  0.015457   0.005408   2.858  0.00426 ** 
LeafToughness     0.706605   1.615106   0.437  0.66175    

Phi coefficients (precision model with identity link):
      Estimate Std. Error z value Pr(>|z|)    
(phi)   43.931      7.828   5.612    2e-08 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Type of estimator: ML (maximum likelihood)
Log-likelihood: 197.5 on 6 Df
Pseudo R-squared: 0.2011
Number of iterations: 27 (BFGS) + 5 (Fisher scoring) 

My questions are:

1- How do I know which original fire treatment corresponds to output's treatment.L and treatment.Q?

2- How do I know the significance (p-value) of the third fire treatment?

3- From this output, can I conclude that nitrogen_content, treatment.L and the third fire treatment mentioned in question 2 have significant effects on the proportion of herbivory?

Many thanks for your thoughts.

$\endgroup$
2
  • $\begingroup$ Your questions seem to arise from your (uncommon) choice to code treatment as an ordered factor. An unordered factor might be more helpful to answer your research questions? $\endgroup$
    – Michael M
    Commented Jul 26, 2019 at 19:56
  • $\begingroup$ Thank you for your reply, @MichaelM. I think my fire treatments need to be ordered because there are a clear ordering of the levels (i.e., fire frequency – unburnt, burnt once, and burnt twice). Hence, I think this variable is not categorical, but indeed ordinal. Any thought on the questions considering I need them as ordered factors? Thanks!! $\endgroup$
    – Lucnp
    Commented Jul 28, 2019 at 16:02

1 Answer 1

1
$\begingroup$

None of the questions is specific to beta regression but is the same for any other model based on a linear predictor, i.e., linear regression, ANOVA, logistic regression etc. Thus, I encourage you to read more about factors and the corresponding contrasts in (generalized) linear models.

  1. Your factor has three levels, thus two coefficients can be freely estimated, the third one needs to be restricted. For an unordered factor the default in R is to use treatment contrasts, i.e., restrict the first coefficient to be zero and then estimate the other two in relation to that. However, for an ordered factor the default is to use a polynomial contrast. This fits a quadratic polynomial through the three points/levels and reports the linear (.L) and quadratic (.Q) coefficient. Here, you see that the linear term is significant while the quadratic is relatively small and not significant. Thus, approximately the effect for two fires is twice that of one fire (albeit a somewhat smaller, as indicated by the negative coefficient in the .Q term). I guess this has a natural interpretation that is easy to summarize.
  2. By itself "significance of the third treatment" is not yet precise enough. What contrast do you want to test? Significantly higher effect for two fires compared to one fire? If so, successive difference contrasts (MASS::contr.sdif) would be better suited than the default polynomial contrasts.
  3. Clearly, nitrogen_content is significant and the effect of the treatment appears to be captured well by a linear effect.

Note also that no matter which conding of the contrasts you choose, the fitted values and hence the maximized log-likelihood remain unchanged. Thus, the argument for or against a certain coding is not goodness of fit but just ease of interpretation (and of coding the contrasts you acually want to test).

$\endgroup$
3
  • 1
    $\begingroup$ I like this a lot. My only nit-pick is the conclusion "approximately the effect for two fires is twice that of one fire" - that depends on the .Q effect being small, not just non-significant! ("the effect of number of fires is not clearly nonlinear ...") $\endgroup$
    – Ben Bolker
    Commented Jul 28, 2019 at 23:35
  • $\begingroup$ Thanks, @BenBolker, I have refined my response by pointing out that the quadratic term is negative and how this can be interpreted. $\endgroup$ Commented Jul 28, 2019 at 23:49
  • $\begingroup$ Dear @AchimZeileis, thank you very much for your time. I assumed it was related to betareg because this was the first time I used both betareg and an ordinal variable– sorry about that. Now I see the output compares linear and quadratic, but not my treatment levels. I discovered that levels in an ordinal variable are assumed to be equally spaced, but I can't assume this for my data. So I decided to run the same model but using treatment as a factor, not ordered factor. As the Control level is already the baseline (alphabetical ordering by default) I got my contrasts with other treatment levels $\endgroup$
    – Lucnp
    Commented Jul 30, 2019 at 21:53

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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