Load packages.
Session info.
sessionInfo()## R version 4.2.3 (2023-03-15)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Big Sur ... 10.16
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] huxtable_5.5.2 ggpubr_0.6.0 nlme_3.1-162 performance_0.10.3
## [5] report_0.5.7 lme4_1.1-32 Matrix_1.5-3 ggbeeswarm_0.7.2
## [9] gridExtra_2.3 knitr_1.42 readxl_1.4.2 pliman_1.1.0
## [13] lubridate_1.9.2 forcats_1.0.0 stringr_1.5.0 dplyr_1.1.1
## [17] purrr_1.0.1 readr_2.1.4 tidyr_1.3.0 tibble_3.2.1
## [21] ggplot2_3.4.2 tidyverse_2.0.0 here_1.0.1
##
## loaded via a namespace (and not attached):
## [1] sass_0.4.5 jsonlite_1.8.4 splines_4.2.3
## [4] carData_3.0-5 bslib_0.4.2 assertthat_0.2.1
## [7] vipor_0.4.5 cellranger_1.1.0 tiff_0.1-11
## [10] yaml_2.3.7 backports_1.4.1 pillar_1.9.0
## [13] lattice_0.20-45 glue_1.6.2 digest_0.6.31
## [16] ggsignif_0.6.4 minqa_1.2.5 colorspace_2.1-0
## [19] htmltools_0.5.5 pkgconfig_2.0.3 broom_1.0.4
## [22] fftwtools_0.9-11 scales_1.2.1 jpeg_0.1-10
## [25] tzdb_0.3.0 timechange_0.2.0 car_3.1-2
## [28] EBImage_4.40.1 generics_0.1.3 cachem_1.0.7
## [31] withr_2.5.0 BiocGenerics_0.44.0 cli_3.6.1
## [34] crayon_1.5.2 magrittr_2.0.3 evaluate_0.20
## [37] fansi_1.0.4 MASS_7.3-58.3 rstatix_0.7.2
## [40] beeswarm_0.4.0 tools_4.2.3 hms_1.1.3
## [43] lifecycle_1.0.3 munsell_0.5.0 locfit_1.5-9.7
## [46] compiler_4.2.3 jquerylib_0.1.4 rlang_1.1.0
## [49] grid_4.2.3 RCurl_1.98-1.12 nloptr_2.0.3
## [52] rstudioapi_0.14 htmlwidgets_1.6.2 bitops_1.0-7
## [55] rmarkdown_2.21 boot_1.3-28.1 gtable_0.3.3
## [58] abind_1.4-5 R6_2.5.1 fastmap_1.1.1
## [61] utf8_1.2.3 rprojroot_2.0.3 insight_0.19.1
## [64] stringi_1.7.12 parallel_4.2.3 Rcpp_1.0.10
## [67] vctrs_0.6.2 png_0.1-8 tidyselect_1.2.0
## [70] xfun_0.38Generate random sequence for assigning cups to ALAN and control groups.
light_day <- c(470,390,560.480,350,760,550,370)
mean(light_day)## [1] 492.9257sd(light_day)## [1] 144.803light_night_ALAN <- c(30,46,47,39)
mean(light_night_ALAN)## [1] 40.5sd(light_night_ALAN)## [1] 7.852813light_night_control <- c(0.2,0.2,0.2,0.2)
mean(light_night_control)## [1] 0.2sd(light_night_control)## [1] 0wilcox.test(light_night_ALAN, light_night_control)## Warning in wilcox.test.default(light_night_ALAN, light_night_control): cannot
## compute exact p-value with ties##
## Wilcoxon rank sum test with continuity correction
##
## data: light_night_ALAN and light_night_control
## W = 16, p-value = 0.02107
## alternative hypothesis: true location shift is not equal to 0Growth data was measured as weekly leaf counts (Day 0 - Day 49) and leaf area was measured from photos taken on Day 28.
Load and prepare main growth data.
dat <- read_csv(here("data", "duckweed_data_main.csv"))## Rows: 161 Columns: 15
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (3): sample_code, group, file_name
## dbl (12): sample_nr, day0_leaves, day7_leaves, day14_leaves, day21_leaves, d...
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.#dim(dat) #161
#names(dat)
dat$group <- as.factor(dat$group) #convert to factor
#levels(dat$group)
levels(dat$group) <- c("ALAN", "Control") #use full group names
dat$group01 <- ifelse(dat$group == "Control", 0, 1) #make a new group variable coded 0 (Control) and 1 (ALAN)
dat <- rename(dat, day28_leaf_area = "day28_TotalLeafArea_cm2)8") #rename variable
dat <- rename(dat, sample = "sample_code") #rename variable
#table(is.na(dat$area)) # 1 missing value
dat$obs <- 1:dim(dat)[1] #add new column with observation numbers
dat$day49_leaves_all <- dat$day49_leaves + dat$day49_deadleaves #calculate total number of leaves (alive + dead) on day 42
dat$day42_leaves_all <- dat$day42_leaves + dat$day42_deadleaves #calculate total number of leaves (alive + dead) on day 42
#sapply(dat, function(x) sum(is.na(x))) #check fo numbers of missing values in each columnLoad and merge processed leaf pignmentation data with count data.
Note: code chunks for actual processing of the images are at the end of this document - they are not run or knitted, as this take long time and the analyses below just use processed and saved output.
leaf_colorcat <- read.csv2(here("data", "leafcolours_severity.csv")) #load processed data
#dim(leaf_colorcat) #159 values
dat_image <- left_join(dat, leaf_colorcat) #join (merge) to main dataframe with experimental group info## Joining with `by = join_by(file_name)`#names(dat_image)
dat_image$day49_prop_dark_area <- dat_image$symptomatic / (dat_image$symptomatic + dat_image$healthy) #calculate proportion of leaves that are dark
#dat_image <- dat_image[!is.na(dat_image$day49_prop_dark_area), ] #remove rows with NA in dat_image$prop_dark_area
#dim(dat_image)
#sapply(dat_image, function(x) sum(is.na(x))) #check fo numbers of missing values in each column
#table(is.na(dat_image$img), dat_image$group) #79 for ALAN and 80 for ControlPrepare weekly leaf count data in a long format.
#names(dat)
## scatterplot plot of counts across the weeks by group using package ggbeeswarm
dat %>% pivot_longer(c(4:10,12), names_to = "days", values_to = "counts") %>%
select(sample, group = group, area = day28_leaf_area, days, counts) ->
dat2 #transform selected columns to long data format
dat2$day <- parse_number(dat2$days) #convert from character string to day number
dat2$zday <- scale(dat2$day) #scale day numbers to mean of 0
dat2$obs <- 1:dim(dat2)[1] #add observation numbers
dat2$group01 <- ifelse(dat2$group == "Control", 0, 1) #recode groups as 0 and 1 values
dat2 <- dat2[!is.na(dat2$counts), ] #remove rows with missing values
#dim(dat2)Prepare data for plotting weekly leaf counts by group .
#names(dat)
#scatterplot plot of counts across the weeks by group using package ggbeeswarm
dat %>% pivot_longer(c(4:10,12), names_to = "days", values_to = "counts") %>%
select(sample, group = group, area = day28_leaf_area, days, counts) ->
dat2 #transform selected columns to long data format
dat2$day <- parse_number(dat2$days) #convert from character string to day number
dat2$zday <- scale(dat2$day) #scale day numbers to mean of 0
dat2$obs <- 1:dim(dat2)[1] #add observation numbers
dat2$group01 <- ifelse(dat2$group == "Control", 0, 1) #recode groups as 0 and 1 values
dat2 <- dat2[!is.na(dat2$counts), ] #remove rows with missing values
#dim(dat2)Code for Figure 3 - plot of weekly leaf counts by group - for the manuscript (plot not shown here - saving into a file instead).
dat2 %>%
ggplot(aes(day, counts, col = group, shape = group)) +
geom_quasirandom(alpha = 0.2, size = 2) + # plot raw data using quasirandom method to distribute the points
stat_summary(fun = mean, geom = 'point', size = 2) + # compute mean points and plot
stat_summary(fun.data = mean_cl_normal, geom = 'errorbar', width = 0) + # compute confidence intervals
scale_color_manual(values = c("sienna1", "deepskyblue3")) +
theme_classic() +
theme(legend.position = "top") +
theme(plot.title = element_text(hjust = 0.5)) +
scale_x_continuous(limits = c(-5, 54), expand = c(0, 0), breaks = c(0,7,14, 21, 28, 35, 42, 49)) +
labs(title = "Duckweed growth as weekly leaf counts") +
labs(y = "Leaf count")
ggsave(here("plots", "fig_weekly_growth.png"), width = 18, height = 12, unit = "cm")Use mixed models for weekly leaf counts.
## using a generalised mixed model with poisson link function and repeated measurements (nesting in sample identity)
model <- glmer(counts ~ 1 + group01 + zday + (1|sample) , family = "poisson", data = dat2) #control is coded as group 0 and ALAN is coded as group 1
summary(model)
report(model)
# report_parameters(model)
#report_statistics(model)
r2_nakagawa(model)Use non-parametric Wilcoxon test to compare medians (without assuming that variances are equal).
# wilcox.test(day0_leaves ~ group, data = dat_image)
# wilcox.test(day7_leaves ~ group, data = dat_image)
# wilcox.test(day14_leaves ~ group, data = dat_image)
# wilcox.test(day21_leaves ~ group, data = dat_image)
# wilcox.test(day28_leaves ~ group, data = dat_image)
# wilcox.test(day35_leaves ~ group, data = dat_image)
# wilcox.test(day42_leaves ~ group, data = dat_image)
# wilcox.test(day49_leaves ~ group, data = dat_image)
# wilcox.test(day42_deadleaves ~ group, data = dat_image)
# wilcox.test(day49_deadleaves ~ group, data = dat_image)
# wilcox.test(day42_leaves_all ~ group, data = dat_image)
# wilcox.test(day49_leaves_all ~ group, data = dat_image)
# wilcox.test(day28_leaf_area ~ group, data = dat_image)
# wilcox.test(day49_prop_dark_area ~ group, data = dat_image)
tab_01 <- t(sapply(dat_image[ , c("day0_leaves","day7_leaves","day14_leaves","day21_leaves","day28_leaves","day35_leaves","day42_leaves","day49_leaves","day42_deadleaves","day49_deadleaves", "day42_leaves_all", "day49_leaves_all","day28_leaf_area", "day49_prop_dark_area")],
function(x) unlist(wilcox.test(x ~ dat$group) [c("statistic","p.value")])))
str(tab_01)## num [1:14, 1:2] 3169 3814 3351 3726 3750 ...
## - attr(*, "dimnames")=List of 2
## ..$ : chr [1:14] "day0_leaves" "day7_leaves" "day14_leaves" "day21_leaves" ...
## ..$ : chr [1:2] "statistic.W" "p.value"tab_01 <- round(as.data.frame(tab_01), 3)Assemble all test results into a table.
kable(
tab_01,
col.names = c("W statistics", "p-value"),
caption = "S1. Results of Wilcoxon runk sum tests to compare medians of two groups."
)| W statistics | p-value | |
|---|---|---|
| day0_leaves | 3169.0 | 0.906 |
| day7_leaves | 3814.5 | 0.031 |
| day14_leaves | 3351.0 | 0.602 |
| day21_leaves | 3726.0 | 0.070 |
| day28_leaves | 3750.5 | 0.059 |
| day35_leaves | 3462.0 | 0.369 |
| day42_leaves | 3459.0 | 0.375 |
| day49_leaves | 3355.0 | 0.596 |
| day42_deadleaves | 4923.0 | 0.000 |
| day49_deadleaves | 4830.5 | 0.000 |
| day42_leaves_all | 3808.0 | 0.037 |
| day49_leaves_all | 3745.0 | 0.062 |
| day28_leaf_area | 4020.5 | 0.005 |
| day49_prop_dark_area | 4738.0 | 0.000 |
Use F-test to test for homogeneity of variances.
#names(dat_image)
# var.test(day0_leaves ~ group, data = dat_image)
# var.test(day7_leaves ~ group, data = dat_image)
# var.test(day14_leaves ~ group, data = dat_image)
# var.test(day21_leaves ~ group, data = dat_image)
# var.test(day28_leaves ~ group, data = dat_image)
# var.test(day35_leaves ~ group, data = dat_image)
# var.test(day42_leaves ~ group, data = dat_image)
# var.test(day49_leaves ~ group, data = dat_image)
# var.test(day42_deadleaves ~ group, data = dat_image)
# var.test(day49_deadleaves ~ group, data = dat_image)
# var.test(day42_leaves_all ~ group, data = dat_image)
# var.test(day49_leaves_all ~ group, data = dat_image)
# var.test(day28_leaf_area ~ group, data = dat_image)
# var.test(day49_prop_dark_area ~ group, data = dat_image)
tab_02 <- t(sapply(dat_image[ , c("day0_leaves","day7_leaves","day14_leaves","day21_leaves","day28_leaves","day35_leaves","day42_leaves","day49_leaves","day42_deadleaves","day49_deadleaves", "day42_leaves_all", "day49_leaves_all","day28_leaf_area", "day49_prop_dark_area")],
function(x) unlist(var.test(x ~ dat$group) [c("statistic","p.value")])))
str(tab_02)## num [1:14, 1:2] 1.04 0.85 1.46 1.46 1.9 ...
## - attr(*, "dimnames")=List of 2
## ..$ : chr [1:14] "day0_leaves" "day7_leaves" "day14_leaves" "day21_leaves" ...
## ..$ : chr [1:2] "statistic.F" "p.value"tab_02 <- round(as.data.frame(tab_02), 3)Assemble all test results into a table.
kable(
tab_02,
col.names = c("F statistics", "p-value"),
caption = "S2. Results of F tests to compare variances of two groups."
)| F statistics | p-value | |
|---|---|---|
| day0_leaves | 1.035 | 0.878 |
| day7_leaves | 0.850 | 0.473 |
| day14_leaves | 1.459 | 0.095 |
| day21_leaves | 1.458 | 0.096 |
| day28_leaves | 1.900 | 0.005 |
| day35_leaves | 1.953 | 0.003 |
| day42_leaves | 1.890 | 0.005 |
| day49_leaves | 2.007 | 0.002 |
| day42_deadleaves | 1.421 | 0.121 |
| day49_deadleaves | 2.158 | 0.001 |
| day42_leaves_all | 1.886 | 0.005 |
| day49_leaves_all | 2.193 | 0.001 |
| day28_leaf_area | 1.814 | 0.009 |
| day49_prop_dark_area | 0.628 | 0.041 |
Results of two tests as one table for the main text (not shown here - can be exported into a file).
Additional model for leaf pigmentation using glm for proportion data.
## fit linear model for proportions
glm_model <- glm(day49_prop_dark_area ~ group, dat_image, family = quasibinomial())
summary(glm_model) #signif difference between groups##
## Call:
## glm(formula = day49_prop_dark_area ~ group, family = quasibinomial(),
## data = dat_image)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.96553 -0.35331 0.01325 0.37514 1.01374
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.0224 0.1069 9.561 < 0.0000000000000002 ***
## groupControl -0.8512 0.1424 -5.976 0.0000000148 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasibinomial family taken to be 0.1757555)
##
## Null deviance: 36.058 on 158 degrees of freedom
## Residual deviance: 29.599 on 157 degrees of freedom
## (2 observations deleted due to missingness)
## AIC: NA
##
## Number of Fisher Scoring iterations: 4Code for Figure 4 for the manuscript - multiple panels with different measures for specific experimental days - (not including the plot here - saving into a file).
## plot for Day 42
dead42 <- dat_image %>%
filter(!is.na(day42_deadleaves)) %>%
ggplot(aes(group, day42_deadleaves, col = group, shape = group)) +
geom_boxplot() +
geom_quasirandom(alpha = 0.2, size = 3) +
theme_classic() +
theme(legend.position = "none") +
theme(plot.title = element_text(hjust = 0.5)) +
ylim(0, 5) +
labs(title = "A. Counts of dead leaves on Day 42") +
scale_fill_manual(values = c("sienna1", "deepskyblue3")) +
labs(y = "Leaf count")
## plot for Day 49
dead49 <- dat_image %>%
filter(!is.na(day49_deadleaves)) %>%
ggplot(aes(group, day49_deadleaves, col = group, shape = group)) +
geom_boxplot() +
geom_quasirandom(alpha = 0.2, size = 3) +
theme_classic() +
theme(legend.position = "none") +
theme(plot.title = element_text(hjust = 0.5)) +
ylim(0, 5) +
labs(title = "B. Counts of dead leaves on Day 49") +
scale_fill_manual(values = c("sienna1", "deepskyblue3")) +
labs(y = "Leaf count")
## plot total leaves for Day 42
all42 <- dat_image %>%
filter(!is.na(day42_leaves_all)) %>%
ggplot(aes(group, day42_leaves_all, col = group, shape = group)) +
geom_boxplot() +
geom_quasirandom(alpha = 0.2, size = 3) +
theme_classic() +
theme(legend.position = "none") +
theme(plot.title = element_text(hjust = 0.5)) +
ylim(0, 40) +
labs(title = "C. Counts of live and dead leaves on Day 42") +
scale_fill_manual(values = c("sienna1", "deepskyblue3")) +
labs(y = "Leaf count")
## plot total leaves for Day 49
all49 <- dat_image %>%
filter(!is.na(day49_leaves_all)) %>%
ggplot(aes(group, day49_leaves_all, col = group, shape = group)) +
geom_boxplot() +
geom_quasirandom(alpha = 0.2, size = 3) +
theme_classic() +
theme(legend.position = "none") +
theme(plot.title = element_text(hjust = 0.5)) +
ylim(0, 40) +
labs(title = "D. Counts of live and dead leaves on Day 49") +
scale_fill_manual(values = c("sienna1", "deepskyblue3")) +
labs(y = "Leaf count")
## plot area for Day 28
area28 <- dat_image %>%
filter(!is.na(day28_leaf_area)) %>%
ggplot(aes(group, day28_leaf_area, col = group, shape = group)) +
geom_boxplot() +
geom_quasirandom(alpha = 0.2, size = 3) +
theme_classic() +
theme(legend.position = "none") +
theme(plot.title = element_text(hjust = 0.5)) +
ylim(0, 1.5) +
labs(title = "E. Total leaf area on Day 28") +
scale_fill_manual(values = c("sienna1", "deepskyblue3")) +
labs(y = "Leaf area [cm^2]")
## plot pigmentation for Day 49
pigm49 <- dat_image %>%
filter(!is.na(day49_prop_dark_area)) %>%
ggplot(aes(group, day49_prop_dark_area, col = group, shape = group)) +
geom_boxplot() +
geom_quasirandom(alpha = 0.2, size = 3) +
theme_classic() +
theme(legend.position = "none") +
theme(plot.title = element_text(hjust = 0.5)) +
ylim(0, 1) +
labs(title = "F. Underside pigmentation on Day 49") +
scale_fill_manual(values = c("sienna1", "deepskyblue3")) +
labs(y = "Proportion")
## save the figure
ggsave(here("plots", "plot.groupdifferences.png"), grid.arrange(dead42, dead49, all42, all49, area28, pigm49, ncol = 2, nrow = 3), width = 16, height = 18, units = "cm", dpi = 300, scale = 1.4) Code for Figure 5 - relationships between leaf counts and area or pigmentation - for the manuscript (not including here - saving into a file).
Correlation between leaf count and leaf area on Day 28, by group.
## Day28 leaf count vs. total leaf area
dat_image %>%
filter(!is.na(day28_leaves)) %>%
filter(!is.na(day28_leaf_area)) %>%
group_by(group) %>%
reframe(est = cor.test(day28_leaves, day28_leaf_area)$estimate,
l.ci = cor.test(day28_leaves, day28_leaf_area)$conf.int[1],
u.ci = cor.test(day28_leaves, day28_leaf_area)$conf.int[2],
p = cor.test(day28_leaves, day28_leaf_area)$p.value)| group | est | l.ci | u.ci | p |
|---|---|---|---|---|
| ALAN | 0.902 | 0.851 | 0.936 | 3.81e-30 |
| Control | 0.855 | 0.782 | 0.905 | 6.09e-24 |
Correlation between total leaf count and proportion of dark leaf area on Day 49, by group.
## Day49 leaf count vs. pigmentation
dat_image %>%
filter(!is.na(day49_leaves)) %>%
filter(!is.na(day49_prop_dark_area)) %>%
group_by(group) %>%
reframe(est = cor.test(day28_leaves, day49_prop_dark_area)$estimate,
l.ci = cor.test(day28_leaves, day49_prop_dark_area)$conf.int[1],
u.ci = cor.test(day28_leaves, day49_prop_dark_area)$conf.int[2],
p = cor.test(day28_leaves, day49_prop_dark_area)$p.value)| group | est | l.ci | u.ci | p |
|---|---|---|---|---|
| ALAN | -0.416 | -0.583 | -0.215 | 0.000137 |
| Control | -0.267 | -0.46 | -0.0503 | 0.0166 |
Note: code not shown in the knitted document - no need to run this code as the processed image data is already saved and loaded for the analyses above.
Calculate pigmentation levels - Use sympmatic_area() function from pliman package to measure pigmentation in 159 samples (batch processing).