library("ggplot2")
library("GGally")
library("Hmisc")
library("readr")
library("dplyr")
library("networkD3")
library("tidyr")
library('tibble')
library('pals')
library("rstatix")
library("cowplot")
library("ggrepel")
sessionInfo()
## R version 4.5.0 (2025-04-11 ucrt)
## Platform: x86_64-w64-mingw32/x64
## Running under: Windows 11 x64 (build 26200)
##
## Matrix products: default
## LAPACK version 3.12.1
##
## locale:
## [1] LC_COLLATE=English_United Kingdom.utf8
## [2] LC_CTYPE=English_United Kingdom.utf8
## [3] LC_MONETARY=English_United Kingdom.utf8
## [4] LC_NUMERIC=C
## [5] LC_TIME=English_United Kingdom.utf8
##
## time zone: Europe/Berlin
## tzcode source: internal
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] ggrepel_0.9.6 cowplot_1.2.0 rstatix_1.1.0 pals_1.10
## [5] tibble_3.3.1 tidyr_1.3.2 networkD3_0.4.1 dplyr_1.2.1
## [9] readr_2.1.6 Hmisc_5.2-6 GGally_2.4.0 ggplot2_4.0.2
##
## loaded via a namespace (and not attached):
## [1] gtable_0.3.6 xfun_0.61 bslib_0.10.0
## [4] htmlwidgets_1.6.4 tzdb_0.5.0 vctrs_0.7.3
## [7] tools_4.5.0 generics_0.1.4 cluster_2.1.8.2
## [10] pkgconfig_2.0.3 data.table_1.18.2.1 checkmate_2.3.4
## [13] RColorBrewer_1.1-3 S7_0.2.1 lifecycle_1.0.5
## [16] compiler_4.5.0 farver_2.1.2 stringr_1.6.0
## [19] data.tree_1.2.0 mapproj_1.2.12 carData_3.0-6
## [22] htmltools_0.5.8.1 maps_3.4.3 sass_0.4.10
## [25] yaml_2.3.10 htmlTable_2.4.3 Formula_1.2-5
## [28] pillar_1.11.1 car_3.1-5 jquerylib_0.1.4
## [31] cachem_1.1.0 abind_1.4-8 rpart_4.1.24
## [34] ggstats_0.12.0 tidyselect_1.2.1 digest_0.6.37
## [37] stringi_1.8.7 purrr_1.2.1 fastmap_1.2.0
## [40] grid_4.5.0 colorspace_2.1-2 cli_3.6.5
## [43] magrittr_2.0.3 base64enc_0.1-6 dichromat_2.0-0.1
## [46] broom_1.0.12 foreign_0.8-91 withr_3.0.2
## [49] scales_1.4.0 backports_1.5.0 rmarkdown_2.30
## [52] igraph_2.2.2 otel_0.2.0 nnet_7.3-20
## [55] gridExtra_2.3 hms_1.1.4 evaluate_1.0.5
## [58] knitr_1.51 rlang_1.2.0 Rcpp_1.1.0
## [61] glue_1.8.0 rstudioapi_0.18.0 jsonlite_2.0.0
## [64] R6_2.6.1
Correlation matrix (Tree growth)
sample_info_no_R <-read.table("Alaska_info_no_R.txt", row.names=1, header=T, check.names=F)
df_correlation <- sample_info_no_R[, c("BAI_10y", "detBAI_10y", "rcsBAI_10y", "Age", "DBH", "Height", "pH", "Site")]
df_correlation <- na.omit(df_correlation)
# Compute correlations + p-values
cor_results <- rcorr(as.matrix(df_correlation[, 1:7])) # column 9 (Site) is used for point colours
cor_matrix <- cor_results$r
p_matrix <- cor_results$P
site_colors <- c(Alaska_Range = "#EE7733", Brooks_Range = "#009988", Interior = "#5E4FA2")
lower_fun_sig <- function(data, mapping, ...) {
ggplot(data = data, mapping = mapping) +
geom_point(aes(color = Site), alpha = 0.7, size = 1.5) +
scale_color_manual(values = site_colors,
labels = c("Alaska Range", "Brooks Range", "Interior Alaska")) +
geom_smooth(method = "lm", se = FALSE, color = "darkred", size = 0.7)}
ggpairs(df_correlation, columns = 1:7,
lower = list(continuous = lower_fun_sig),
upper = list(continuous = wrap("cor", size = 4)),
diag = list(continuous = wrap("densityDiag", alpha = 0.5, fill = "lightblue")),
legend = c(2, 1),
progress = FALSE) +
theme_bw() + theme(panel.grid = element_blank(),
strip.background = element_rect(fill = "lightgray"),
axis.text = element_text(size = 8),
legend.position = "bottom")
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
Differences between Plots
palette <- c("#EE7733","#CC6677", "#009988","#337538","#33BBEE", "#5E4FA2")
pH <- ggplot(sample_info_no_R, aes(x = Plot, y = pH, fill = Plot)) + geom_violin(trim = FALSE) +geom_point(size = 2, color = "gray10", alpha = 0.3,
position = position_jitter(width = 0.1)) + scale_fill_manual(values = palette) + ylab("Soil pH") + theme_minimal() + theme(legend.position = "none")
Age_plot <- ggplot(sample_info_no_R, aes(x = Plot, y = Age, fill = Plot)) + geom_violin(trim = FALSE) +geom_point(size = 2, color = "gray10", alpha = 0.3,
position = position_jitter(width = 0.1)) + scale_fill_manual(values = palette) + ylab("Tree age, years") + theme_minimal() + theme(legend.position = "none")
DBH_plot <- ggplot(sample_info_no_R, aes(x = Plot, y = DBH, fill = Plot)) + geom_violin(trim = FALSE) +geom_point(size = 2, color = "gray10", alpha = 0.3,
position = position_jitter(width = 0.1)) + scale_fill_manual(values = palette) + ylab("DBH, cm") + theme_minimal() + theme(legend.position = "none")
Height_plot <- ggplot(sample_info_no_R, aes(x = Plot, y = Height, fill = Plot)) + geom_violin(trim = FALSE) +geom_point(size = 2, color = "gray10", alpha = 0.3,
position = position_jitter(width = 0.1)) + scale_fill_manual(values = palette) + ylab("Tree height, m") + theme_minimal() + theme(legend.position = "none")
plot_grid(pH, Age_plot, DBH_plot, Height_plot, labels = c("Soil pH", "Tree age", "Tree DBH", "Tree height" ), ncol = 2)
## Warning: Removed 3 rows containing non-finite outside the scale range
## (`stat_ydensity()`).
## Warning: Removed 3 rows containing missing values or values outside the scale range
## (`geom_point()`).
## Warning: Removed 2 rows containing non-finite outside the scale range
## (`stat_ydensity()`).
## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_point()`).
## Warning: Removed 2 rows containing non-finite outside the scale range
## (`stat_ydensity()`).
## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_point()`).
Soil pH in BF (map)
BF_pH <- subset(sample_info_no_R, Plot == "BF" & !is.na(pH))
BF_pH$Sample <- rownames(BF_pH)
ggplot(BF_pH, aes(x = Longitude, y = Latitude)) +
geom_point(aes(fill = pH), shape = 21, size = 4, color = "gray20") +
geom_text_repel(aes(label = paste0(Sample, "\npH ", pH)), size = 2.8, point.size = 4, max.overlaps = Inf, lineheight = 0.9) +
scale_fill_viridis_c(name = "Soil pH") +
coord_quickmap() +
theme_bw()
Statistical tests
shapiro.test(sample_info_no_R$pH) # data are not normally distributed
##
## Shapiro-Wilk normality test
##
## data: sample_info_no_R$pH
## W = 0.90744, p-value = 1.111e-06
bartlett.test(pH ~ Plot, data = sample_info_no_R) # variance is not equal
##
## Bartlett test of homogeneity of variances
##
## data: pH by Plot
## Bartlett's K-squared = 70.137, df = 5, p-value = 9.597e-14
dunn_test(sample_info_no_R , pH ~ Habitat) #not significant
## # A tibble: 1 × 9
## .y. group1 group2 n1 n2 statistic p p.adj p.adj.signif
## * <chr> <chr> <chr> <int> <int> <dbl> <dbl> <dbl> <chr>
## 1 pH forest treeline 59 52 0.287 0.774 0.774 ns
dunn_test(sample_info_no_R , pH ~ Plot, p.adjust.method = "bonferroni")
## # A tibble: 15 × 9
## .y. group1 group2 n1 n2 statistic p p.adj p.adj.signif
## * <chr> <chr> <chr> <int> <int> <dbl> <dbl> <dbl> <chr>
## 1 pH AF AT 19 21 -1.47 1.41e- 1 1 e+ 0 ns
## 2 pH AF BF 19 21 3.56 3.71e- 4 5.56e- 3 **
## 3 pH AF BT 19 20 6.66 2.78e-11 4.17e-10 ****
## 4 pH AF IF 19 19 3.93 8.53e- 5 1.28e- 3 **
## 5 pH AF IT 19 11 2.92 3.48e- 3 5.22e- 2 ns
## 6 pH AT BF 21 21 5.16 2.43e- 7 3.64e- 6 ****
## 7 pH AT BT 21 20 8.32 8.89e-17 1.33e-15 ****
## 8 pH AT IF 21 19 5.50 3.83e- 8 5.75e- 7 ****
## 9 pH AT IT 21 11 4.23 2.37e- 5 3.55e- 4 ***
## 10 pH BF BT 21 20 3.22 1.29e- 3 1.93e- 2 *
## 11 pH BF IF 21 19 0.466 6.41e- 1 1 e+ 0 ns
## 12 pH BF IT 21 11 -0.0544 9.57e- 1 1 e+ 0 ns
## 13 pH BT IF 20 19 -2.68 7.39e- 3 1.11e- 1 ns
## 14 pH BT IT 20 11 -2.73 6.28e- 3 9.42e- 2 ns
## 15 pH IF IT 19 11 -0.443 6.58e- 1 1 e+ 0 ns
dunn_test(sample_info_no_R , pH ~ Site, p.adjust.method = "bonferroni")
## # A tibble: 3 × 9
## .y. group1 group2 n1 n2 statistic p p.adj p.adj.signif
## * <chr> <chr> <chr> <int> <int> <dbl> <dbl> <dbl> <chr>
## 1 pH Alaska_Range Brook… 40 41 8.38 5.25e-17 1.58e-16 ****
## 2 pH Alaska_Range Inter… 40 30 6.04 1.57e- 9 4.72e- 9 ****
## 3 pH Brooks_Range Inter… 41 30 -1.68 9.22e- 2 2.77e- 1 ns
shapiro.test(sample_info_no_R$Age) # data are not normally distributed
##
## Shapiro-Wilk normality test
##
## data: sample_info_no_R$Age
## W = 0.92442, p-value = 1.209e-05
bartlett.test(Age ~ Habitat, data = sample_info_no_R) # variance is equal
##
## Bartlett test of homogeneity of variances
##
## data: Age by Habitat
## Bartlett's K-squared = 1.0282, df = 1, p-value = 0.3106
kruskal_test(sample_info_no_R , Age ~ Habitat) #significant
## # A tibble: 1 × 6
## .y. n statistic df p method
## * <chr> <int> <dbl> <int> <dbl> <chr>
## 1 Age 108 8.43 1 0.00369 Kruskal-Wallis
shapiro.test(sample_info_no_R$Age) # data are not normally distributed
##
## Shapiro-Wilk normality test
##
## data: sample_info_no_R$Age
## W = 0.92442, p-value = 1.209e-05
bartlett.test(Age ~ Plot, data = sample_info_no_R) # variance is not equal
##
## Bartlett test of homogeneity of variances
##
## data: Age by Plot
## Bartlett's K-squared = 54.315, df = 5, p-value = 1.805e-10
res.kruskal <-kruskal_test(sample_info_no_R , Age ~ Plot)
res.kruskal #significant
## # A tibble: 1 × 6
## .y. n statistic df p method
## * <chr> <int> <dbl> <int> <dbl> <chr>
## 1 Age 108 22.4 5 0.000440 Kruskal-Wallis
dunn_test(sample_info_no_R , Age ~ Plot, p.adjust.method = "bonferroni")
## # A tibble: 15 × 9
## .y. group1 group2 n1 n2 statistic p p.adj p.adj.signif
## * <chr> <chr> <chr> <int> <int> <dbl> <dbl> <dbl> <chr>
## 1 Age AF AT 19 18 -3.50 0.000465 0.00697 **
## 2 Age AF BF 19 21 -1.43 0.154 1 ns
## 3 Age AF BT 19 20 -2.82 0.00479 0.0718 ns
## 4 Age AF IF 19 19 -3.52 0.000432 0.00648 **
## 5 Age AF IT 19 11 -3.50 0.000459 0.00688 **
## 6 Age AT BF 18 21 2.18 0.0294 0.441 ns
## 7 Age AT BT 18 20 0.762 0.446 1 ns
## 8 Age AT IF 18 19 0.0282 0.977 1 ns
## 9 Age AT IT 18 11 -0.460 0.645 1 ns
## 10 Age BF BT 21 20 -1.45 0.148 1 ns
## 11 Age BF IF 21 19 -2.18 0.0292 0.439 ns
## 12 Age BF IT 21 11 -2.35 0.0186 0.279 ns
## 13 Age BT IF 20 19 -0.744 0.457 1 ns
## 14 Age BT IT 20 11 -1.13 0.259 1 ns
## 15 Age IF IT 19 11 -0.490 0.624 1 ns
DBH vs Habitat/Plot
shapiro.test(sample_info_no_R$DBH) # data are not normally distributed
##
## Shapiro-Wilk normality test
##
## data: sample_info_no_R$DBH
## W = 0.92854, p-value = 1.919e-05
bartlett.test(DBH ~ Habitat, data = sample_info_no_R) # variance is not equal
##
## Bartlett test of homogeneity of variances
##
## data: DBH by Habitat
## Bartlett's K-squared = 5.1545, df = 1, p-value = 0.02319
dunn_test(sample_info_no_R , DBH ~ Habitat) # not significant
## # A tibble: 1 × 9
## .y. group1 group2 n1 n2 statistic p p.adj p.adj.signif
## * <chr> <chr> <chr> <int> <int> <dbl> <dbl> <dbl> <chr>
## 1 DBH forest treeline 59 50 -1.51 0.131 0.131 ns
shapiro.test(sample_info_no_R$DBH) # data are not normally distributed
##
## Shapiro-Wilk normality test
##
## data: sample_info_no_R$DBH
## W = 0.92854, p-value = 1.919e-05
bartlett.test(DBH ~ Plot, data = sample_info_no_R) # variance is not equal
##
## Bartlett test of homogeneity of variances
##
## data: DBH by Plot
## Bartlett's K-squared = 23.109, df = 5, p-value = 0.0003217
res.kruskal <-kruskal_test(sample_info_no_R , DBH ~ Plot)
res.kruskal #significant
## # A tibble: 1 × 6
## .y. n statistic df p method
## * <chr> <int> <dbl> <int> <dbl> <chr>
## 1 DBH 109 30.0 5 0.0000149 Kruskal-Wallis
dunn_test(sample_info_no_R , DBH ~ Plot, p.adjust.method = "bonferroni")
## # A tibble: 15 × 9
## .y. group1 group2 n1 n2 statistic p p.adj p.adj.signif
## * <chr> <chr> <chr> <int> <int> <dbl> <dbl> <dbl> <chr>
## 1 DBH AF AT 19 19 -4.37 0.0000126 0.000188 ***
## 2 DBH AF BF 19 21 -2.26 0.0238 0.358 ns
## 3 DBH AF BT 19 20 -2.72 0.00657 0.0985 ns
## 4 DBH AF IF 19 19 -2.78 0.00536 0.0805 ns
## 5 DBH AF IT 19 11 0.616 0.538 1 ns
## 6 DBH AT BF 19 21 2.22 0.0267 0.401 ns
## 7 DBH AT BT 19 20 1.71 0.0881 1 ns
## 8 DBH AT IF 19 19 1.58 0.113 1 ns
## 9 DBH AT IT 19 11 4.36 0.0000133 0.000199 ***
## 10 DBH BF BT 21 20 -0.497 0.619 1 ns
## 11 DBH BF IF 21 19 -0.593 0.553 1 ns
## 12 DBH BF IT 21 11 2.55 0.0108 0.162 ns
## 13 DBH BT IF 20 19 -0.102 0.919 1 ns
## 14 DBH BT IT 20 11 2.94 0.00327 0.0491 *
## 15 DBH IF IT 19 11 3.00 0.00270 0.0405 *
Height vs Habitat/Plot
shapiro.test(sample_info_no_R$Height) # data are not normally distributed
##
## Shapiro-Wilk normality test
##
## data: sample_info_no_R$Height
## W = 0.94182, p-value = 0.0001269
bartlett.test(Height ~ Habitat, data = sample_info_no_R) # variance is equal
##
## Bartlett test of homogeneity of variances
##
## data: Height by Habitat
## Bartlett's K-squared = 2.6759, df = 1, p-value = 0.1019
dunn_test(sample_info_no_R , Height ~ Habitat) # significant
## # A tibble: 1 × 9
## .y. group1 group2 n1 n2 statistic p p.adj p.adj.signif
## * <chr> <chr> <chr> <int> <int> <dbl> <dbl> <dbl> <chr>
## 1 Height forest treeline 59 50 -4.08 0.0000455 0.0000455 ****
shapiro.test(sample_info_no_R$Height) # data are not normally distributed
##
## Shapiro-Wilk normality test
##
## data: sample_info_no_R$Height
## W = 0.94182, p-value = 0.0001269
bartlett.test(Height ~ Plot, data = sample_info_no_R) # variance is not equal
##
## Bartlett test of homogeneity of variances
##
## data: Height by Plot
## Bartlett's K-squared = 35.151, df = 5, p-value = 1.404e-06
res.kruskal <-kruskal_test(sample_info_no_R , Height ~ Plot)
res.kruskal #significant
## # A tibble: 1 × 6
## .y. n statistic df p method
## * <chr> <int> <dbl> <int> <dbl> <chr>
## 1 Height 109 47.6 5 0.00000000424 Kruskal-Wallis
dunn_test(sample_info_no_R , Height ~ Plot, p.adjust.method = "bonferroni")
## # A tibble: 15 × 9
## .y. group1 group2 n1 n2 statistic p p.adj p.adj.signif
## * <chr> <chr> <chr> <int> <int> <dbl> <dbl> <dbl> <chr>
## 1 Height AF AT 19 19 -4.74 0.00000217 3.25e-5 ****
## 2 Height AF BF 19 21 -1.52 0.129 1 e+0 ns
## 3 Height AF BT 19 20 -3.37 0.000739 1.11e-2 *
## 4 Height AF IF 19 19 0.272 0.786 1 e+0 ns
## 5 Height AF IT 19 11 1.08 0.280 1 e+0 ns
## 6 Height AT BF 19 21 3.34 0.000845 1.27e-2 *
## 7 Height AT BT 19 20 1.42 0.155 1 e+0 ns
## 8 Height AT IF 19 19 5.01 0.000000546 8.19e-6 ****
## 9 Height AT IT 19 11 5.14 0.000000277 4.16e-6 ****
## 10 Height BF BT 21 20 -1.92 0.0545 8.17e-1 ns
## 11 Height BF IF 21 19 1.80 0.0726 1 e+0 ns
## 12 Height BF IT 21 11 2.39 0.0168 2.52e-1 ns
## 13 Height BT IF 20 19 3.65 0.000262 3.93e-3 **
## 14 Height BT IT 20 11 3.97 0.0000714 1.07e-3 **
## 15 Height IF IT 19 11 0.848 0.396 1 e+0 ns