QK Box 7.2

Breitwieser et al. (2016) surveyed variegated scallops, Mimachlamys varia, living at several sites along the French Atlantic coast. Four sites were chosen (factor site); three sites were chosen because they were potentially contaminated, with different sources of contaminants, and a fourth site was considered relatively clean (and considered a reference site). The authors sampled scallops at two times of the year (March and September), chosen to correspond to before and at the end of the scallop’s reproductive season. These two times were levels of the factor season (although strictly these levels represent two different sampling times that may or may not reflect seasonal differences). Breitwieser and colleagues measured several variables to assess the condition of scallops, and here we use their data using a biomarker, Malondialdehyde (MDA). MDA is a stress marker and was measured (μM/g fresh tissue) in ten scallops from each combination of site and season. A two-factor linear model (7.3) including the fixed main effects of site and season and their interaction was fitted to these data.

Dugornay Olivier (2020). Ifremer. CC BY-SA 4.0

Breitwieser, M., et al. (2016). Short-term and long-term biological effects of chronic chemical contamination on natural populations of a marine bivalve. PLoS One, 11, e0150184. The paper is here

Preliminaries

First, load the required packages (car, effectsizes, afex) + ggplot2, patchwork

Import breit data file (breit.csv)

breit <- read.csv("../data/breit.csv")
#For convenience, the original data file has one site as Port-Neuf. R dislikes arithmetic operators in characters, so recode site
breit$site <- case_match(
    breit$site,
    "Port-Neuf" ~ "Port Neuf",
    .default = breit$site
    )
head(breit,10)

Initial look at raw data

Start with boxplot

boxplot(mda~group,data=breit)

Fit model to untransformed data and check residuals

breit.aov<-aov(mda~site*season, data=breit)
plot(breit.aov)

Transform response to logs & redo diagnostics

Boxplots

breit$lmda<- log10(breit$mda) 
boxplot(lmda~group,data=breit)

Refit model with log data and check residuals

lbreit.aov<-aov(lmda~site*season, data=breit)
plot(lbreit.aov)

Analysis with untransformed data

Some oddities in the boxplots, but no real improvement when log-transforming, so go with original data

summary(breit.aov)

            Df Sum Sq Mean Sq F value  Pr(>F)   
site         3  13946    4649   3.088 0.03244 * 
season       1  14018   14018   9.313 0.00319 **
site:season  3  13131    4377   2.908 0.04038 * 
Residuals   72 108378    1505                   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Get effect size measures (sjstats package)

effectsize(breit.aov)

# Effect Size for ANOVA (Type I)

Parameter   | Eta2 (partial) |       95% CI
-------------------------------------------
site        |           0.11 | [0.01, 1.00]
season      |           0.11 | [0.02, 1.00]
site:season |           0.11 | [0.00, 1.00]

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

Interaction plot

interaction.plot(breit$site,breit$season,breit$mda)

afex_plot(breit.aov, "site", "season", dodge=0.1, data_plot=FALSE)+theme_light()

dv column detected: mda
No id column passed. Assuming all rows are independent samples.

Simple main effects - effect of site in March

breitmar.aov <- aov(mda~site, data = breit, subset = c(season == 'march')) 
# Anova() function from car package tests simple main effect against original error term from 2 factor model 
Anova(lm(breitmar.aov), lm(breit.aov), type=3)

Anova Table (Type III tests)

Response: mda
            Sum Sq Df  F value    Pr(>F)    
(Intercept) 157314  1 104.5100 1.153e-15 ***
site          5524  3   1.2232    0.3075    
Residuals   108378 72                       
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Effect of site in September

breitsep.aov <- aov(mda~site, data = breit, subset=c(season == 'sept')) 
Anova(lm(breitsep.aov), lm(breit.aov), type = 3)

Anova Table (Type III tests)

Response: mda
            Sum Sq Df F value    Pr(>F)    
(Intercept)  99509  1 66.1074 8.744e-12 ***
site         21553  3  4.7727  0.004328 ** 
Residuals   108378 72                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Now include contrast between Loix and average of rest

Note that the book has an error and the original code we provided treated Les Palles, rather than Loix, as the reference cite.

breit$site <- factor(breit$site)
contrasts(breit$site)<- c(1,-3,1,1)    #This is a correction from the original listed contrast of -3/1/1/1
contrasts(breit$site)

           [,1]          [,2]          [,3]
Les Palles    1 -5.773503e-01 -5.773503e-01
Loix         -3 -2.775558e-17  2.775558e-17
Minimes       1  7.886751e-01 -2.113249e-01
Port Neuf     1 -2.113249e-01  7.886751e-01

breit2.aov <- aov(mda~site*season, data = breit)

summary.lm(breit2.aov)

Call:
aov(formula = mda ~ site * season, data = breit)

Residuals:
    Min      1Q  Median      3Q     Max 
-60.828 -24.762  -3.249  14.991 171.602 

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)      133.7895     6.1344  21.810  < 2e-16 ***
site1              4.1408     3.5417   1.169  0.24620    
site2             18.3149    12.2689   1.493  0.13986    
site3              3.3449    12.2689   0.273  0.78591    
seasonsept       -26.4748     8.6754  -3.052  0.00319 ** 
site1:seasonsept   0.6624     5.0088   0.132  0.89515    
site2:seasonsept -36.3224    17.3508  -2.093  0.03984 *  
site3:seasonsept  36.0776    17.3508   2.079  0.04115 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 38.8 on 72 degrees of freedom
Multiple R-squared:  0.2749,    Adjusted R-squared:  0.2044 
F-statistic:   3.9 on 7 and 72 DF,  p-value: 0.001176

Figure

Residual plot

p1<-ggplot(breit.aov, aes(x = breit.aov$fitted.values, y = breit.aov$residuals)) +
  geom_point(color=sc) +
  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<-p1+ annotate ("text", x=95, y=170, label="Untransformed")
p1

p3<-ggplot(lbreit.aov, aes(x = lbreit.aov$fitted.values, y = lbreit.aov$residuals)) +
  geom_point(color=sc) +
  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"
       )
p3<-p3 + annotate ("text", x=1.96, y=0.35, label="Log-transformed")
p3

Interaction plot

breit2<-summarySE(data=breit,measurevar="mda", groupvars=c("site", "season"))
p5<-ggplot(breit2, aes(x = site, y = mda, fill = season)) +
  geom_errorbar(aes(ymin = mda-se, ymax = mda+se), width=0, position = position_dodge(width=0.6)) +
  geom_bar(stat = "identity", width=.6, position = position_dodge(width=0.6), color=lc) +
  scale_fill_uchicago(na.value = "red",labels = c("March", "September"))+
  scale_y_continuous(expand = c(0,0))+
  scale_x_discrete(labels = str_wrap(c("Les Palles", "Loix", "Minimes", "Port Neuf"), width = 7))+
  labs(x = NULL, y = "MDA")+
    theme_qk()+
    theme(
    axis.text.x = element_text(colour = lc, size=8),
    legend.key.size = unit(0.35, 'cm'), #change legend key size
#        legend.key.height = unit(1, 'cm'), #change legend key height
#        legend.key.width = unit(1, 'cm'), #change legend key width
 #       legend.title = element_text(size=14), #change legend title font size
#        legend.text = element_text(size=10)) #change legend text font size
  legend.position = "right",
  legend.title = element_blank(),
  legend.text=element_text(size=7),
  )+
theme(
  legend.position = c(.5, .95),
  legend.justification = c("right", "top"),
  legend.box.just = "right",
  legend.margin = margin(6, 6, 6, 6),
  legend.title = element_blank()
)
p5

Combine figures

p1+p3+p5+plot_layout(widths = c(1,1,2))