QK Box 7.5

Long and Porturas (2014) examined the effect of multiple stressors on the performance of a saltmarsh plant that is important for ecological restoration. Their focus was the potential for salinity stress to modify the herbivory coming from scale insects, and they experimentally removed scale insects or left them intact (Factor: Scale), on plots with salinity at ambient levels or elevated (Factor: Salinity). The experiment was repeated at two sites chosen to be very different in overall elevation (Factor: Site, a fixed effect in this context) within a marsh in southern California. These three factors form a three-way factorial design, and from each replicate experimental plot (of which there were 7 and 8 at the two sites). Their response variable was the time to senescence of a single stem of the cordgrass Spartina foliosa.

Spartina foliosa. Pacific Southwest Region U.S. Fish and Wildlife Service, Public domain, via Wikimedia Commons

The paper is here

Long, J. D. & Porturas, L. D. (2014). Herbivore impacts on marsh production depend upon a compensatory continuum mediated by salinity stress. PLoS One, 9, e110419.

Preliminaries

First, load the required packages (car, sjstats, afex, pwr); apaTables added for cleaner output

Import longport data file (longport.csv)

longport <- read.csv("../data/longport.csv")
head(longport,10)

Fit model to untransformed data and check residuals

longport.aov <- aov(days~site*salinity*scale, data=longport)
plot(longport.aov)

Residuals fine so examine analysis with untransformed data

summary(longport.aov)

                    Df Sum Sq Mean Sq F value Pr(>F)   
site                 1   6801    6801    7.58 0.0081 **
salinity             1    491     491    0.55 0.4629   
scale                1    450     450    0.50 0.4818   
site:salinity        1    563     563    0.63 0.4319   
site:scale           1    195     195    0.22 0.6427   
salinity:scale       1   1530    1530    1.71 0.1974   
site:salinity:scale  1   5182    5182    5.78 0.0199 * 
Residuals           51  45740     897                  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Get Type III SS

longport2 <- Anova(lm(longport.aov), type=3)
longport2

Anova Table (Type III tests)

Response: days
                    Sum Sq Df F value  Pr(>F)    
(Intercept)          80143  1   89.36 8.5e-13 ***
site                     1  1    0.00   0.974    
salinity              6174  1    6.88   0.011 *  
scale                 3590  1    4.00   0.051 .  
site:salinity         4482  1    5.00   0.030 *  
site:scale            1783  1    1.99   0.165    
salinity:scale        6320  1    7.05   0.011 *  
site:salinity:scale   5182  1    5.78   0.020 *  
Residuals            45740 51                    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Get effect size measures (sjstats package) - point estimates match SPSS

effectsize(longport2)

Type 3 ANOVAs only give sensible and informative results when covariates are mean-centered and factors are coded with
  orthogonal contrasts (such as those produced by `contr.sum`, `contr.poly`, or `contr.helmert`, but *not* by the default
  `contr.treatment`).

# Effect Size for ANOVA (Type III)

Parameter           | Eta2 (partial) |       95% CI
---------------------------------------------------
site                |       2.04e-05 | [0.00, 1.00]
salinity            |           0.12 | [0.02, 1.00]
scale               |           0.07 | [0.00, 1.00]
site:salinity       |           0.09 | [0.00, 1.00]
site:scale          |           0.04 | [0.00, 1.00]
salinity:scale      |           0.12 | [0.02, 1.00]
site:salinity:scale |           0.10 | [0.01, 1.00]

- One-sided CIs: upper bound fixed at [1.00].

Simple interaction effects - salinity x scale at North

One way to approach this 3-factor interaction is to look at the interaction between two predictors at each level of the third. For this experiment, salinity and scale are biological effects of interest, and site measures whether those effects are stable across space, so it makes sense to examine each site separately.

One way to do this is to analyze each site separately using, e.g. aov:

longportN.aov <- aov(days~salinity*scale, data=longport, subset=c(site=='North'))
summary(longportN.aov)

               Df Sum Sq Mean Sq F value Pr(>F)  
salinity        1   1073    1073    1.11  0.303  
scale           1     14      14    0.01  0.905  
salinity:scale  1   6320    6320    6.53  0.018 *
Residuals      23  22269     968                 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

This approach creates a residual MS specific to this data subset. An alternative approach is to use the the original error term from the 3 factor model. Anova() function from car package tests simple main effect against whole model residual.

longportN1.aov <- Anova(lm(longportN.aov), lm(longport.aov), type=3)
longportN1.aov

Anova Table (Type III tests)

Response: days
               Sum Sq Df F value  Pr(>F)    
(Intercept)     80143  1   89.36 8.5e-13 ***
salinity         6174  1    6.88   0.011 *  
scale            3590  1    4.00   0.051 .  
salinity:scale   6320  1    7.05   0.011 *  
Residuals       45740 51                    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

effectsize(longportN1.aov)

Type 3 ANOVAs only give sensible and informative results when covariates are mean-centered and factors are coded with
  orthogonal contrasts (such as those produced by `contr.sum`, `contr.poly`, or `contr.helmert`, but *not* by the default
  `contr.treatment`).

# Effect Size for ANOVA (Type III)

Parameter      | Eta2 (partial) |       95% CI
----------------------------------------------
salinity       |           0.12 | [0.02, 1.00]
scale          |           0.07 | [0.00, 1.00]
salinity:scale |           0.12 | [0.02, 1.00]

- One-sided CIs: upper bound fixed at [1.00].

Simple interaction effects - salinity x scale at South

longportS.aov <- aov(days~salinity*scale, data=longport, subset=c(site=='South'))
longportS1.aov <- Anova(lm(longportS.aov), lm(longport.aov), type=3)
longportS1.aov

Anova Table (Type III tests)

Response: days
               Sum Sq Df F value  Pr(>F)    
(Intercept)     92450  1  103.08 7.7e-14 ***
salinity          196  1    0.22    0.64    
scale              12  1    0.01    0.91    
salinity:scale    392  1    0.44    0.51    
Residuals       45740 51                    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

effectsize(longportS1.aov)

Type 3 ANOVAs only give sensible and informative results when covariates are mean-centered and factors are coded with
  orthogonal contrasts (such as those produced by `contr.sum`, `contr.poly`, or `contr.helmert`, but *not* by the default
  `contr.treatment`).

# Effect Size for ANOVA (Type III)

Parameter      | Eta2 (partial) |       95% CI
----------------------------------------------
salinity       |       4.27e-03 | [0.00, 1.00]
scale          |       2.68e-04 | [0.00, 1.00]
salinity:scale |       8.50e-03 | [0.00, 1.00]

- One-sided CIs: upper bound fixed at [1.00].

Generate figure

Residual plot

p1<-ggplot(longport.aov, aes(x = longport.aov$fitted.values, y = longport.aov$residuals)) +
geom_point(color=sc, alpha=0.5) +
theme_classic(base_size = 10)+
theme(
axis.text = element_text(colour = ac),
axis.line = element_line(color = ac),
axis.ticks = element_line(color = ac),
)+labs(x = "Predicted", y = "Residuals"
)

p1

Interaction plots

Use emmeans to make file of means and se

emm1<-emmeans(longport.aov, ~site|salinity|scale)
emm2<-as.data.frame(emm1)
emm2

salinity = Ambient, scale = NoScale:
 site  emmean   SE df lower.CL upper.CL
 North  107.0 11.3 51     84.3    129.7
 South  107.5 10.6 51     86.2    128.8

salinity = Enhanced, scale = NoScale:
 site  emmean   SE df lower.CL upper.CL
 North   65.0 11.3 51     42.3     87.7
 South  114.5 10.6 51     93.2    135.8

salinity = Ambient, scale = Scale:
 site  emmean   SE df lower.CL upper.CL
 North   73.7 12.2 51     49.1     98.2
 South  105.7 10.6 51     84.5    127.0

salinity = Enhanced, scale = Scale:
 site  emmean   SE df lower.CL upper.CL
 North   93.0 11.3 51     70.3    115.7
 South   98.7 10.6 51     77.5    120.0

Confidence level used: 0.95 

Separate interaction plots for N and S Means and error bars Use filter to subset data for ggplot

emm3<-subset(emm2,site=="North")
pd=position_dodge(width=0)
p2<-ggplot(emm3,aes(x=salinity,y=emmean,shape=scale, group=scale, color=scale))+
  geom_point(position=pd,aes(shape=scale), size=3,show.legend = FALSE)+
  geom_errorbar(aes(ymin = emmean-SE, ymax = emmean+SE), width=0, position = pd,show.legend = FALSE)+
  geom_line(aes(color=scale), position=pd, size=1.5, show.legend = FALSE)+
  scale_color_uchicago(labels = c("No scale", "Scale"))+
    scale_linetype_manual(values=c("solid", "solid"))+
  labs(x = "Salinity", y = "Days to senescence", title="North Site"
       )+
  ylim(50,150)+
  scale_x_discrete(labels=c("Ambient", "Enhanced"))+
  theme_classic(base_size = 10)+
  theme(
    axis.text.x = element_text(color="black"),
    axis.text.y= element_text(color=ac),
    axis.line = element_line(color = ac),
    axis.ticks = element_line(color = ac),
        )+
  theme(
  legend.position = c(.9, .50),
  legend.justification = c("right", "top"),
  legend.box.just = "right",
  legend.margin = margin(6, 6, 6, 6),
  plot.title = element_text(hjust=0.5),
  legend.title=element_blank()
)
  
p2

Now S site

emm4<-subset(emm2, site=="South")
pd=position_dodge(width=0)
p3<-ggplot(emm4, aes(x=salinity,y=emmean,shape=scale, group=scale, color=scale))+
  geom_point(position=pd,aes(shape=scale), size=3,show.legend = FALSE)+
  geom_errorbar(aes(ymin = emmean-SE, ymax = emmean+SE), width=0, position = pd,show.legend = FALSE)+
  geom_line(aes(color=scale), position=pd, size=1.5)+
  scale_color_uchicago(labels = c("No scale", "Scale"))+
    scale_linetype_manual(values=c("solid", "solid"))+
  labs(x = "Salinity", y = NULL, title="South Site"
       )+
  ylim(50,150)+
  scale_x_discrete(labels=c("Ambient", "Enhanced"))+
  theme_classic(base_size = 10)+
  theme(
    axis.text.x = element_text(color="black"),
    axis.text.y= element_text(color=ac),
    axis.line = element_line(color = ac),
    axis.ticks = element_line(color = ac),
        )+
  theme(
  legend.position = c(.8, .4),
  legend.justification = c("right", "top"),
  legend.box.just = "right",
  legend.margin = margin(6, 6, 6, 6),
  plot.title = element_text(hjust=0.5),
  legend.title=element_blank()
)
  
p3

p1+p2+p3