1
$\begingroup$

I have some experimental data which I am analysing using step wise multiple regression (ANCOVA) in R using the step function. The response data (wp) is the leaf water potential of a tree which has been subject to varying degrees of water stress (Treatment). I have two covariates in the model. 1. "Road" - the presence (Y) or absence (N) of a road in the tree's vicinity, and 2. "soil" - the volumetric soil moisture data (continuous variable) for an individual tree at the time that the response data were collected.

I am investigating the effects of Treatment, Road and soil moisture on the response data as well as the interactions between Treatment:Road and Treatment:soil. The Treatment effects are the main point of interest. For some of the data, a simple anovaor lm reveals significant Treatment effects on the response, however when I include all the terms and the interactions shown above, the step function leaves me with a model (based on AIC) which includes a significant Treatment:soil interaction term and no individual Treatment effects.

Why? And, what do I need to do, adopting these statistical tests, to be able to describe the Treatment effects. I've done some reading re: Type III Anova and changing the contrasts in the analysis, but I'm not too savvy with how I should go about this. Is that the right way to go? How?

Here is the data which I'm looking at:

 Treatment<-c(6,12,3,"CONTROL",12,3,"CONTROL",6,3,12,"CONTROL",6,3,
             6,3,3,12,6,12,6,"CONTROL",12,12,6,"CONTROL",3,3,
             "CONTROL",6,"CONTROL",12)

Road<-c("N","Y","N","N","N","N","Y","N","Y","Y","Y","Y","Y","Y","Y",
        "Y","Y","Y","Y","Y","Y","N","N","N","N","N","N","N","N","N","N")

wp<-c(-0.325,-0.225,-0.375,-0.275,-0.3625,-0.375,-0.2,-0.4625,
      -0.375,-0.325,-0.25,-0.3,-0.4,-0.35,-0.55,-0.5,-0.375,-0.2,
      -0.3,-0.3,-0.25,-0.3,-0.375,-0.3,-0.35,-0.5,-0.475,-0.3,-0.5,
      -0.2,-0.35)

soil<-c(18.992299,  20.3859736, 19.4265055, 19.0402522, 19.3498457,
        18.1948846, 21.7836259, 20.3867353, 19.6153346, 21.6668146,
        17.8964699, 16.4279241, 19.1379134, 18.2698171, 18.2698171,
        18.8901119, 19.438544,  18.7045546, 17.1389654, 18.570092,
        18.8455254, 19.580172,  23.5295579, 18.6212624, 25.6860396,
        23.6555276, 21.7282271, 23.3053829, 21.9061206, 23.5122382,
        24.6748561)

wp<-data.frame(Treatment,Road,wp,soil)
wp$Treatment<-factor(wp$Treatment,levels=c("CONTROL",12,6,3))

And here is the code which I'm running to do the step wise ANCOVA

mod_1<-lm(wp ~ Treatment + Road + soil + Treatment:Road + Treatment:soil, data = wp)
summary(mod_1)

step(mod_1,direction="both")

The final model selection / last part of the step result

Step:  AIC=-164.05
wp ~ Treatment + soil + Treatment:soil

                 Df Sum of Sq      RSS     AIC
<none>                        0.093106 -164.05
- Treatment:soil  3 0.0222411 0.115347 -163.41
+ Road            1 0.0005513 0.092555 -162.23

Call:
lm(formula = wp ~ Treatment + soil + Treatment:soil, data = wp)

Coefficients:
     (Intercept)       Treatment12        Treatment6        Treatment3  
      -1.189e-01        -7.246e-02         6.285e-01        -1.326e-01  
            soil  Treatment12:soil   Treatment6:soil   Treatment3:soil  
      -6.616e-03         8.964e-05        -3.825e-02        -3.062e-03

The final model

mod_2<-lm(wp ~ Treatment + soil + Treatment:soil, data = wp)
summary(mod_2)

And the output from the summary command

Call:
lm(formula = wp ~ Treatment + soil + Treatment:soil, data = wp)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.121687 -0.035023 -0.006441  0.031295  0.129611 

Coefficients:
                   Estimate Std. Error t value Pr(>|t|)  
(Intercept)      -1.189e-01  1.926e-01  -0.617   0.5431  
Treatment12      -7.246e-02  2.807e-01  -0.258   0.7986  
Treatment6        6.285e-01  3.449e-01   1.822   0.0814 .
Treatment3       -1.326e-01  3.189e-01  -0.416   0.6814  
soil             -6.616e-03  8.912e-03  -0.742   0.4654  
Treatment12:soil  8.964e-05  1.324e-02   0.007   0.9947  
Treatment6:soil  -3.825e-02  1.747e-02  -2.190   0.0390 *
Treatment3:soil  -3.062e-03  1.555e-02  -0.197   0.8457  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.06362 on 23 degrees of freedom
Multiple R-squared:  0.6501,    Adjusted R-squared:  0.5435 
F-statistic: 6.104 on 7 and 23 DF,  p-value: 0.0004149

The simple ANOVA for Treatment only effects

mod_3<-lm(wp~Treatment, data = wp)
summary(mod_3)

And the summary output table

Call:
lm(formula = wp ~ Treatment, data = wp)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.15781 -0.04386  0.01071  0.04297  0.14219 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -0.26071    0.02678  -9.735 2.52e-10 ***
Treatment12 -0.06585    0.03667  -1.796   0.0838 .  
Treatment6  -0.08147    0.03667  -2.222   0.0349 *  
Treatment3  -0.18304    0.03667  -4.991 3.12e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.07086 on 27 degrees of freedom
Multiple R-squared:  0.4905,    Adjusted R-squared:  0.4339 
F-statistic: 8.663 on 3 and 27 DF,  p-value: 0.0003438

Biologically, the treatment effects are performing as I would expect here.

Setting the contrasts per suggestion here produces different results.

options(contrasts = c("contr.sum","contr.poly"))

mod_2b<-lm(formula = wp ~ Treatment + soil + Treatment:soil, data = wp)
summary(mod_2b)

And the summary output table

Call:
lm(formula = wp ~ Treatment + soil + Treatment:soil, data = wp)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.121687 -0.035023 -0.006441  0.031295  0.129611 

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)   
(Intercept)     -0.013022   0.118665  -0.110  0.91357   
Treatment1      -0.105856   0.180618  -0.586  0.56353   
Treatment2      -0.178311   0.186907  -0.954  0.35000   
Treatment3       0.522639   0.234561   2.228  0.03593 * 
soil            -0.016922   0.005935  -2.851  0.00904 **
Treatment1:soil  0.010306   0.008657   1.191  0.24599   
Treatment2:soil  0.010396   0.009122   1.140  0.26617   
Treatment3:soil -0.027946   0.012170  -2.296  0.03111 * 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.06362 on 23 degrees of freedom
Multiple R-squared:  0.6501,    Adjusted R-squared:  0.5435 
F-statistic: 6.104 on 7 and 23 DF,  p-value: 0.0004149

My understanding of this now, is that the Treatments shown in the summary table are not ordered as factors, as was set originally. These are now ordered as they appear in the data.frame, i.e. 6, 12, 3, Control

I've looked here, but I'm unsure if a). it's applicable to my scenario and b). if so, how to implement it correctly.

I'm probably more conversant with R, than I am with stats, so forgive me if I've explained this poorly. Any help and suggestions greatly received.

Thanks

Improve this question
$\endgroup$

2 Answers 2

Reset to default
1
$\begingroup$

I suspect that one issue you're having is that you want to conduct an anova and you are looking at the summary output, which doesn't show any anova results.

It's fine to use the default contrasts in R, and Type-II sum of squares with library(car); Anova(mod_2).

The interpretation of the soil:treatment interaction described by @a_statistician is right-on, and I suspect it is what you would expect for results. The leaf water potential is related to the soil moisture content, but this relationship varies among treatments. If this is your only significant effect, it's pretty easy to plot as wp vs. soil for each treatment, preferably on one plot, with each treatment as a different color. You can tease out the intercept and slopes for each line out of the coefficients listed in the summary of the model. I'll refer you to the chapter on ancova of the same site you cited, with the caveat that I am the author.

Improve this answer
$\endgroup$
4
  • $\begingroup$ Thanks both!! This study was carried out over some time, and I have several responses which I am anaylsing this way. It's for my Ph.D. work and the advice / instructions I am receiving from the university is that I should 'clean up the models', by removing any non-significant terms. @Sal Mangiafico, I had begun to play with the car package and run the Anova command. In doing so, from the example above, the ANOVA shows the Treatment:soil interaction as non-significant. When I drop that term, then all terms returned in the subsequent model remain significant. Is that the right way? $\endgroup$
    – A.Benson
    Oct 12, 2018 at 18:31
  • $\begingroup$ That's a valid approach. Also, because you have a classic analysis of covariance (ancova), it leads to a simple graphical explanation: if treatment and soil are significant, and the interaction is not, then the lines for wp ~ soil have different intercepts, but not-different slopes. $\endgroup$
    – Sal Mangiafico
    Oct 12, 2018 at 18:51
  • $\begingroup$ Thanks.Could you explain, in lay terms, why I use Type II and not Type III SS in this instance please? $\endgroup$
    – A.Benson
    Oct 12, 2018 at 19:24
  • $\begingroup$ No, I probably can't. I'm going to punt and let you find some discussion from someone more knowledgeable than I am about types of sums of squares. $\endgroup$
    – Sal Mangiafico
    Oct 12, 2018 at 19:45
0
$\begingroup$

Your final model includes that soil, treatment and their interaction. It means the four straight lines between soil and wp are fitted, one line for one treatment. Especially, 4 lines are not parallel. (You can try to draw 4 lines on a piece of paper based on estimated parameters.) It means the differences of wp between treatments are not constant, and they vary alone the soil values. It is possible that the difference between treatment are significant at soil = x, and not significant when soil = y. Under this situation, you need to specify the value of soil, then test if the difference of wp between treatment are significant or not at specified soil value.

The model with options(contrasts = c("contr.sum","contr.poly")) is equivalent to model without options(contrasts = c("contr.sum","contr.poly")).

Improve this answer
$\endgroup$

Your Answer

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

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