UVR consequences on polyp reproduction

Importing and tidying asexual reproduction data

Reproduction Results 1. One-way ANOVA (parametric); Tukey Post-Hoc 2. Two-way ANOVA (parametric); Tukey Post-Hoc 3. Kruskall-Wallace (non-parametric); Pairwise Wilcox Test

Detachment Results 1. Chi-squaredtest of independence

  1. Importing data.
# data where all the controls are combined into one group 
df.combinedcontrol <- read.csv(here::here("data_raw/polypclonedata_combinedcontrol.csv"), stringsAsFactors = F)
  1. Add column with final polyp numbers
df.combinedcontrol$FinalPolyp <- df.combinedcontrol$Day_27
  1. Calculate daily and weekly budding rates
# daily budding rate for experiment
df.combinedcontrol$DBR <- ((df.combinedcontrol$FinalPolyp - df.combinedcontrol$Day_00)/27)
# weekly budding rate for experiment
df.combinedcontrol$WBR <- (((df.combinedcontrol$FinalPolyp - df.combinedcontrol$Day_00)/27)*7)
  1. Reorganize/Reformat the data frame
# long format 
df.combinedcontrol.gather <- df.combinedcontrol %>% gather(key = "Time", value = "Polyps", Day_00:Day_27)

# separate the days of experiment out and only keep the number of the day
df <- separate(data = df.combinedcontrol.gather, col = Time, into = c(NA,'Day'), sep = '_')

# reformat number of day to integer
df$Day <- as.integer(df$Day)

knitr::kable(head(df))
Treatment JarID DishID FinalPolyp DBR WBR Day Polyps
Control C1 D1 11 0.3703704 2.5925926 0 1
Control C1 D2 4 0.1111111 0.7777778 0 1
Control C1 D3 8 0.2592593 1.8148148 0 1
Control C2 D4 9 0.2962963 2.0740741 0 1
Control C2 D5 10 0.3333333 2.3333333 0 1
Control C2 D6 7 0.2222222 1.5555556 0 1

Data analysis

We compared differences in weekly budding rate among the UVR treatment groups using a one-way ANOVA with treatment as the main effect.

wbraov <- aov(WBR~Treatment, data = df.combinedcontrol)

summary(wbraov)
##             Df Sum Sq Mean Sq F value   Pr(>F)    
## Treatment    3  23.43   7.809   10.94 1.23e-05 ***
## Residuals   50  35.69   0.714                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

A post-hoc Tukey test was used when significant differences were detected.

TukeyHSD(wbraov)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = WBR ~ Treatment, data = df.combinedcontrol)
## 
## $Treatment
##                  diff          lwr         upr     p adj
## AB-AA      -1.8580247 -2.774636818 -0.94141256 0.0000113
## BB-AA      -0.9506173 -1.867229411 -0.03400516 0.0393725
## Control-AA -0.4609053 -1.297653914  0.37584321 0.4665907
## BB-AB       0.9074074 -0.009204719  1.82401953 0.0532755
## Control-AB  1.3971193  0.560370778  2.23386791 0.0002857
## Control-BB  0.4897119 -0.347036630  1.32646050 0.4130403

UVR consequences on polyp health

Creating table of detachment data

Polyp detachment was classified binomially as either yes (polyps in dish detached at > 1 timepoints) or no (polyps in dish detached at ≤ 1 timepoints).

D <- as.table(rbind(c(11,7), c(11,1), c(4,8), c(0,12)))

dimnames(D) <- list(treatment = c("Control", "AA", "BB", "AB"), status = c("Attached", "Detached"))

knitr::kable(D)
Attached Detached
Control 11 7
AA 11 1
BB 4 8
AB 0 12

Data analysis

Chi-squared test of independence to test whether there were non-random differences in detachment probabilities across the treatment groups.

library(chisq.posthoc.test)

chisq.test(D, correct = F)
## 
##  Pearson's Chi-squared test
## 
## data:  D
## X-squared = 22.512, df = 3, p-value = 5.103e-05

Post-Hoc Comparisons

6 comparisons so alpha is 0.05/6 = 0.0083
Comparison is bolded if p value is significant

Control vs AA

chisq.test(D[c(1,2),], correct = F)
## 
##  Pearson's Chi-squared test
## 
## data:  D[c(1, 2), ]
## X-squared = 3.4375, df = 1, p-value = 0.06373

Control vs BB

chisq.test(D[c(1,3),], correct = F)
## 
##  Pearson's Chi-squared test
## 
## data:  D[c(1, 3), ]
## X-squared = 2.2222, df = 1, p-value = 0.136

Control vs AB

chisq.test(D[c(1,4),], correct = F)
## 
##  Pearson's Chi-squared test
## 
## data:  D[c(1, 4), ]
## X-squared = 11.579, df = 1, p-value = 0.000667

AA vs BB

chisq.test(D[c(2,3),], correct = F)
## 
##  Pearson's Chi-squared test
## 
## data:  D[c(2, 3), ]
## X-squared = 8.7111, df = 1, p-value = 0.003163

AA vs AB

chisq.test(D[c(2,4),], correct = F)
## 
##  Pearson's Chi-squared test
## 
## data:  D[c(2, 4), ]
## X-squared = 20.308, df = 1, p-value = 6.593e-06

BB vs AB

chisq.test(D[c(3,4),], correct = F)
## 
##  Pearson's Chi-squared test
## 
## data:  D[c(3, 4), ]
## X-squared = 4.8, df = 1, p-value = 0.02846

Figures

Fig 1A. Total number of polyps observed (mean +/- s.e.) under each UV treatment over the 27-day experiment.

Fig 1B. Weekly budding rates (polyps/week) for each UV treatment. Lowercase letters (a, b, c) indicate signficant differences (p < 0.5) between groups.

Data Analysis in Response to Reviewers

  1. Importing data. Added UVA and UVB column independently. ‘1’ indicates presence of that wavelength and ‘0’ indicates absence.
# data where all the controls are combined into one group 
df <- read.csv(here::here("data_raw/polypclonedata_combinedcontrol_v2.csv"), stringsAsFactors = F)
  1. Add column with final polyp numbers
df$FinalPolyp <- df$Day_27
  1. Calculate daily and weekly budding rates
# daily budding rate for experiment
df$DBR <- ((df$FinalPolyp - df$Day_00)/27)
# weekly budding rate for experiment
df$WBR <- (((df$FinalPolyp - df$Day_00)/27)*7)
  1. Reorganize/Reformat the data frame
# Make UVA and UVB factors
df$UVA <- as.factor(df$UVA)
df$UVB <- as.factor(df$UVB)

knitr::kable(head(df))
Treatment UVA UVB JarID DishID Day_00 Day_01 Day_02 Day_03 Day_04 Day_05 Day_06 Day_07 Day_09 Day_11 Day_13 Day_15 Day_17 Day_19 Day_21 Day_23 Day_25 Day_27 FinalPolyp DBR WBR
Control 0 0 C1 D1 1 1 1 1 3 3 3 3 3 4 4 5 5 8 8 9 9 11 11 0.3703704 2.5925926
Control 0 0 C1 D2 1 1 1 1 2 2 2 2 2 2 3 3 3 4 4 4 3 4 4 0.1111111 0.7777778
Control 0 0 C1 D3 1 1 1 1 2 2 2 2 2 3 3 3 4 5 7 7 8 8 8 0.2592593 1.8148148
Control 0 0 C2 D4 1 1 1 1 3 3 3 3 3 3 3 3 4 5 5 6 7 9 9 0.2962963 2.0740741
Control 0 0 C2 D5 1 1 1 1 2 2 3 3 3 4 4 4 6 6 9 10 10 10 10 0.3333333 2.3333333
Control 0 0 C2 D6 1 1 1 1 2 2 2 2 2 2 2 2 2 3 3 3 4 7 7 0.2222222 1.5555556

Data analysis

interaction <- aov(WBR ~ UVA * UVB, data = df)

plot(interaction)

Shapiro-Wilk Test of Normality. WBR departs significantly from parametric assumption of normality.

plot(interaction, 2)

rstatix::shapiro_test(residuals(interaction))
## # A tibble: 1 × 3
##   variable               statistic p.value
##   <chr>                      <dbl>   <dbl>
## 1 residuals(interaction)     0.925 0.00240

Levene’s of homogeneity of variances. WBR departs significantly from parametric assumption.

library(rstatix)
## 
## Attaching package: 'rstatix'
## The following object is masked from 'package:stats':
## 
##     filter
plot(interaction, 1)

levene_test(WBR~UVA * UVB, data = df)
## # A tibble: 1 × 4
##     df1   df2 statistic       p
##   <int> <int>     <dbl>   <dbl>
## 1     3    50      4.97 0.00429

Because the data failed parametric assumption, we decided to run an non-paramentric Kruskal-Wallis Rank Sum Test. Given that there was a signficant effect, we ran Pairwise Wilcox test with a Bonferroni correction to correct for multiple comparisons.

kruskal.test(WBR ~ Treatment, data = df)
## 
##  Kruskal-Wallis rank sum test
## 
## data:  WBR by Treatment
## Kruskal-Wallis chi-squared = 31.488, df = 3, p-value = 6.71e-07
pairwise.wilcox.test(df$WBR, df$Treatment, p.adjust.method = "bonferroni")
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties
## 
##  Pairwise comparisons using Wilcoxon rank sum test with continuity correction 
## 
## data:  df$WBR and df$Treatment 
## 
##         AA      AB      BB     
## AB      6.0e-05 -       -      
## BB      0.04955 0.00021 -      
## Control 1.00000 1.5e-05 0.64963
## 
## P value adjustment method: bonferroni

We compared the relative effects of UVA and UVB on weekly budding rate using a two-way ANOVA, that included an interaction between UVA and UVB. We also created a model that considered JarID since it was a grouping factor.

interaction <- aov(WBR~UVA * UVB, data = df)
two.way <-aov(WBR~UVA + UVB, data = df)
blocking <- aov(WBR ~ UVA * UVB + JarID, data = df)

We compared the AICs of the models to see which is of best fit.

library(AICcmodavg)
## Warning: package 'AICcmodavg' was built under R version 4.2.2
model.set <- list(two.way, interaction, blocking)
model.names <- c("two.way", "interaction", "blocking")

aictab(model.set, modnames = model.names)
## 
## Model selection based on AICc:
## 
##              K   AICc Delta_AICc AICcWt Cum.Wt     LL
## interaction  5 142.13       0.00   0.95   0.95 -65.44
## two.way      4 148.25       6.12   0.04   1.00 -69.72
## blocking    19 155.51      13.38   0.00   1.00 -47.58
summary(interaction)
##             Df Sum Sq Mean Sq F value   Pr(>F)    
## UVA          1   0.99   0.988   1.384   0.2449    
## UVB          1  16.31  16.313  22.855 1.58e-05 ***
## UVA:UVB      1   6.13   6.127   8.585   0.0051 ** 
## Residuals   50  35.69   0.714                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

A post-hoc Tukey test was used when significant differences were detected.

TukeyHSD(interaction)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = WBR ~ UVA * UVB, data = df)
## 
## $UVA
##           diff        lwr       upr     p adj
## 1-0 -0.2722222 -0.7369398 0.1924953 0.2449369
## 
## $UVB
##          diff       lwr       upr    p adj
## 1-0 -1.100556 -1.565273 -0.635838 1.71e-05
## 
## $`UVA:UVB`
##               diff        lwr          upr     p adj
## 1:0-0:0  0.4609053 -0.3758432  1.297653914 0.4665907
## 0:1-0:0 -0.4897119 -1.3264605  0.347036630 0.4130403
## 1:1-0:0 -1.3971193 -2.2338679 -0.560370778 0.0002857
## 0:1-1:0 -0.9506173 -1.8672294 -0.034005157 0.0393725
## 1:1-1:0 -1.8580247 -2.7746368 -0.941412565 0.0000113
## 1:1-0:1 -0.9074074 -1.8240195  0.009204719 0.0532755
summary(blocking)
##             Df Sum Sq Mean Sq F value   Pr(>F)    
## UVA          1  0.988   0.988   1.931  0.17314    
## UVB          1 16.313  16.313  31.887 2.06e-06 ***
## JarID       15 23.398   1.560   3.049  0.00308 ** 
## Residuals   36 18.417   0.512                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1