Chukchi Sea fish 09092025
Franz Mueter
2025-09-10
Data imports & clean-up
Assuming data files are in the working directory
arctic_traits_090424<-as.data.frame(read.csv("/Users/lauren/downloads/arctic_traits_090424_b.csv",
row.names = 1))
View(arctic_traits_090424)
names(arctic_traits_090424)[20]<-"mature.adult.1"
arctic_cpue_090424<-as.data.frame(read.csv("/Users/lauren/downloads/biomass.transformed.matrix.csv",
row.names = 1))
View(arctic_cpue_090424)
arctic_env_090424<-as.data.frame(read.csv("/Users/lauren/downloads/arctic_env_090424_b.csv",
row.names = 1))
##choose trait taxa from cpue
arctic_traits_090424_b<-as.data.frame(arctic_traits_090424) %>%
filter(rownames(arctic_traits_090424) %in%
colnames(arctic_cpue_090424))
View(arctic_traits_090424_b)
##goes to 33 taxa
##select stations in arctic_env_090424
arctic_cpue_090424.2 <- as.data.frame(arctic_cpue_090424) %>%
filter(row.names(arctic_cpue_090424) %in%
row.names(arctic_env_090424))
View(arctic_cpue_090424.2)##33 taxa, 431 stations
dim(arctic_env_090424) ## 7 environmental variables, 431 stations## [1] 431 7##create weighted traits matrix
arctic_traits_090424_b<-arctic_traits_090424_b[,-21]##removed adult 2 bc none were adult 2
arctic.traits <- as.data.frame(as.matrix(arctic_cpue_090424.2) %*%
as.matrix(arctic_traits_090424_b))
View(arctic.traits)
##create fuzzy coded traits matrix
# 'prep.fuzzy' standardizes modalities within each trait to sum to 1:
arctic.traits.fuzzy<-as.data.frame(prep.fuzzy(arctic.traits,
c(3,6,3,5,3,3)))
View(arctic.traits.fuzzy)
# Replace NAs (resulting from traits that have a 0 for all modalities)
arctic.traits.fuzzy[is.na(arctic.traits.fuzzy)] <- 0Assess multicollinearity:
Cor.all<-cor(arctic_env_090424)
corrplot(Cor.all, method="number", order = "original",
addCoef.col = "black", diag = FALSE, is.corr = TRUE,
type="lower",tl.col = "black", shade.col = "grey50")
Results show moderate correlation between salinity and temperature but do not suggest a need to eliminate any of the variables from the analysis a priori
Cluster analysis
Conduct hierarchical cluster analysis based on functional traits and Gower distances
fuzzy.dis <- vegdist(arctic.traits.fuzzy, method="gower")
fuzzy.clus <- hclust(fuzzy.dis, method = "ward.D2")
dend <- as.dendrogram(fuzzy.clus)
plot(dend)
# Cut the tree to get 3 clusters (cut off at 1.0)
groups <- as.data.frame(cutree(dend, k = 3))
View(groups)
names(groups)[1]<-"cluster"
groups$cluster<-as.factor(groups$cluster)
summary(groups)## cluster
## 1:225
## 2:139
## 3: 67mycol<-c("#FD1C26", "#0D7EF9", "#98A100")
dend_colored_branches <- color_branches(dend, k = 3,
groupLabels = FALSE,col=mycol)
dend_colored_branches<-set(dend_colored_branches,"labels_color", 'white')
plot(dend_colored_branches)
Plot results from cluster analysis
CWM stacked barplot
Create stacked barplots based on community-weighted means for each modality
CWM-RDA
Conduct redundancy analysis based on community-weighted means of functional traits. Start with a model including the depth gradient, dynamic environmental variables (temperature, day of ice retreat and salinity) and latitude/longitude to account for any remaining (linear) spatial gradients.
RDA_model.arctic <- rda(arctic.traits.fuzzy ~
Depth + doy_retreat + BottomTemp +
BottomSalinity + Latitude + Longitude,
data = arctic_env)
vif.cca(RDA_model.arctic) ## Depth doy_retreat BottomTemp BottomSalinity Latitude
## 1.185658 1.580370 2.902462 1.998060 1.971497
## Longitude
## 1.577640RDA_model.arctic## Call: rda(formula = arctic.traits.fuzzy ~ Depth + doy_retreat + BottomTemp
## + BottomSalinity + Latitude + Longitude, data = arctic_env)
##
## -- Model Summary --
##
## Inertia Proportion Rank
## Total 0.26772 1.00000
## Constrained 0.04483 0.16746 6
## Unconstrained 0.22289 0.83254 17
##
## Inertia is variance
##
## -- Eigenvalues --
##
## Eigenvalues for constrained axes:
## RDA1 RDA2 RDA3 RDA4 RDA5 RDA6
## 0.03341 0.00877 0.00147 0.00083 0.00027 0.00009
##
## Eigenvalues for unconstrained axes:
## PC1 PC2 PC3 PC4 PC5 PC6 PC7 PC8
## 0.07220 0.04100 0.03337 0.02047 0.01305 0.01201 0.00845 0.00634
## (Showing 8 of 17 unconstrained eigenvalues)RsquareAdj(RDA_model.arctic)##0.17## $r.squared
## [1] 0.1674586
##
## $adj.r.squared
## [1] 0.1556773plot(RDA_model.arctic)
The full model shows a gradient from warm, low salinity, shallow stations with an early ice retreat to the opposite, and a spatial gradient from SW to NE (Bering Strait to NE Chukchi Sea). This spatial gradient may reflect something like “distance from Bering Strait”(?)
Assess significance
Use permutation-based analysis of variance to assess the statistical significance of the model overall, of RDA axes and of individual terms:
anova(RDA_model.arctic)##significant## Permutation test for rda under reduced model
## Permutation: free
## Number of permutations: 999
##
## Model: rda(formula = arctic.traits.fuzzy ~ Depth + doy_retreat + BottomTemp + BottomSalinity + Latitude + Longitude, data = arctic_env)
## Df Variance F Pr(>F)
## Model 6 0.044832 14.214 0.001 ***
## Residual 424 0.222885
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(RDA_model.arctic, by="axis") ##RDA1, RDA2 as significant## Permutation test for rda under reduced model
## Forward tests for axes
## Permutation: free
## Number of permutations: 999
##
## Model: rda(formula = arctic.traits.fuzzy ~ Depth + doy_retreat + BottomTemp + BottomSalinity + Latitude + Longitude, data = arctic_env)
## Df Variance F Pr(>F)
## RDA1 1 0.033408 63.5531 0.001 ***
## RDA2 1 0.008766 16.6762 0.001 ***
## RDA3 1 0.001468 2.7928 0.243
## RDA4 1 0.000828 1.5758 0.587
## RDA5 1 0.000269 0.5124 0.975
## RDA6 1 0.000091 0.1737 0.992
## Residual 424 0.222885
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(RDA_model.arctic, by="margin") ##Bottom salinity not significant## Permutation test for rda under reduced model
## Marginal effects of terms
## Permutation: free
## Number of permutations: 999
##
## Model: rda(formula = arctic.traits.fuzzy ~ Depth + doy_retreat + BottomTemp + BottomSalinity + Latitude + Longitude, data = arctic_env)
## Df Variance F Pr(>F)
## Depth 1 0.001999 3.8037 0.007 **
## doy_retreat 1 0.001124 2.1382 0.067 .
## BottomTemp 1 0.002485 4.7277 0.001 ***
## BottomSalinity 1 0.000334 0.6349 0.695
## Latitude 1 0.005665 10.7760 0.001 ***
## Longitude 1 0.005650 10.7488 0.001 ***
## Residual 424 0.222885
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1The model is highly significant overall, and the first two RDA axes are significant (intrepretable).
All marginal effects (significance of individual terms after all other terms are included) were significant except for salinity. The largest variations in trait composition were associated with Temperature and the remaining spatial gradient (latitude/longitude). However, because there are likely other year-to-year differences that are not captured by the covariates, it would be best to use permutations within years for the permutation (rather than across all samples, which is the default). This is similar to considering year as a random effect as residual year-to-year differences are essentially removed by this approach:
anova(RDA_model.arctic, by="margin",
permutations = how(nperm=4999, blocks = arctic_env$Year))## Permutation test for rda under reduced model
## Marginal effects of terms
## Blocks: arctic_env$Year
## Permutation: free
## Number of permutations: 4999
##
## Model: rda(formula = arctic.traits.fuzzy ~ Depth + doy_retreat + BottomTemp + BottomSalinity + Latitude + Longitude, data = arctic_env)
## Df Variance F Pr(>F)
## Depth 1 0.001999 3.8037 0.0280 *
## doy_retreat 1 0.001124 2.1382 0.6396
## BottomTemp 1 0.002485 4.7277 0.0182 *
## BottomSalinity 1 0.000334 0.6349 0.7614
## Latitude 1 0.005665 10.7760 0.0002 ***
## Longitude 1 0.005650 10.7488 0.0002 ***
## Residual 424 0.222885
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1The results are roughly the same, but the day of ice retreat is no longer significant, which is likely because ice retreat differs much more between years than among stations within years, hence it can ‘absorbs’ some of the year to year differences.
plot(doy_retreat ~ factor(Year), data=, arctic_env)
This suggest that a year effect may be appropriate to include in the model, although it needs to be interpreted with caution because of the spatially unbalanced sampling design. Interannual differences in ice retreat are likely confounded with other sources of interannual variability and do not appear to provide the best explanation of such interannual differences when both year and ice retreat are included.
Moreover, ice retreat is also confounded with temperature (r = -0.54), thus its effects are highly uncertain
Step-wise regression
We conduct a step-wise regression to potentially simplify the model objectively (using forward & backward selection based on adjusted \(R^2\) values).
RDA_model.arctic.2 <- rda( arctic.traits.fuzzy ~ 1,
data = arctic_env)
RDA_mod.step <- ordiR2step(RDA_model.arctic.2,
scope = formula(RDA_model.arctic),
direction="both")## Step: R2.adj= 0
## Call: arctic.traits.fuzzy ~ 1
##
## R2.adjusted
## <All variables> 0.15567733
## + Latitude 0.10095382
## + BottomTemp 0.08203750
## + Longitude 0.06724984
## + doy_retreat 0.04270619
## + BottomSalinity 0.02466407
## + Depth 0.01524609
## <none> 0.00000000
##
## Df AIC F Pr(>F)
## + Latitude 1 -611.85 49.285 0.002 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: R2.adj= 0.1009538
## Call: arctic.traits.fuzzy ~ Latitude
##
## R2.adjusted
## <All variables> 0.1556773
## + BottomTemp 0.1246677
## + Longitude 0.1243322
## + Depth 0.1147063
## + BottomSalinity 0.1122931
## + doy_retreat 0.1100879
## <none> 0.1009538
##
## Df AIC F Pr(>F)
## + BottomTemp 1 -622.38 12.622 0.002 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: R2.adj= 0.1246677
## Call: arctic.traits.fuzzy ~ Latitude + BottomTemp
##
## R2.adjusted
## <All variables> 0.1556773
## + Longitude 0.1480764
## + Depth 0.1319377
## + BottomSalinity 0.1273105
## + doy_retreat 0.1269279
## <none> 0.1246677
##
## Df AIC F Pr(>F)
## + Longitude 1 -633.07 12.76 0.002 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: R2.adj= 0.1480764
## Call: arctic.traits.fuzzy ~ Latitude + BottomTemp + Longitude
##
## R2.adjusted
## <All variables> 0.1556773
## + Depth 0.1539179
## + doy_retreat 0.1502608
## <none> 0.1480764
## + BottomSalinity 0.1480588
##
## Df AIC F Pr(>F)
## + Depth 1 -635.05 3.9481 0.006 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: R2.adj= 0.1539179
## Call: arctic.traits.fuzzy ~ Latitude + BottomTemp + Longitude + Depth
##
## R2.adjusted
## + doy_retreat 0.1564027
## <All variables> 0.1556773
## <none> 0.1539179
## + BottomSalinity 0.1534161RDA_mod.step## Call: rda(formula = arctic.traits.fuzzy ~ Latitude + BottomTemp + Longitude
## + Depth, data = arctic_env)
##
## -- Model Summary --
##
## Inertia Proportion Rank
## Total 0.26772 1.00000
## Constrained 0.04331 0.16179 4
## Unconstrained 0.22440 0.83821 17
##
## Inertia is variance
##
## -- Eigenvalues --
##
## Eigenvalues for constrained axes:
## RDA1 RDA2 RDA3 RDA4
## 0.03330 0.00825 0.00127 0.00050
##
## Eigenvalues for unconstrained axes:
## PC1 PC2 PC3 PC4 PC5 PC6 PC7 PC8
## 0.07241 0.04154 0.03340 0.02067 0.01318 0.01235 0.00846 0.00636
## (Showing 8 of 17 unconstrained eigenvalues)plot(RDA_mod.step)
The model eliminates the ice retreat variable and had a similar \(R^2\) value (16.2%)
Environment only model
To assess how much of the variability can (at most) be explained by depth and dynamic ocean variables, we fit a model without the spatial pattern:
RDA_model.arctic.env <- rda(arctic.traits.fuzzy ~
Depth + doy_retreat + BottomTemp +
BottomSalinity,
data = arctic_env)
RDA_model.arctic.env## Call: rda(formula = arctic.traits.fuzzy ~ Depth + doy_retreat + BottomTemp
## + BottomSalinity, data = arctic_env)
##
## -- Model Summary --
##
## Inertia Proportion Rank
## Total 0.26772 1.00000
## Constrained 0.02835 0.10590 4
## Unconstrained 0.23937 0.89410 17
##
## Inertia is variance
##
## -- Eigenvalues --
##
## Eigenvalues for constrained axes:
## RDA1 RDA2 RDA3 RDA4
## 0.023871 0.003214 0.000904 0.000362
##
## Eigenvalues for unconstrained axes:
## PC1 PC2 PC3 PC4 PC5 PC6 PC7 PC8
## 0.08286 0.04372 0.03338 0.02078 0.01399 0.01312 0.00846 0.00677
## (Showing 8 of 17 unconstrained eigenvalues)(R2.1 <- RsquareAdj(RDA_model.arctic.env)$r.squ)## [1] 0.1058986About 10.6% of the variability in functional traits can (at most) be accounted for by bathymetry and dynamic ocean variables
Spatial model
To assess how much of the variability can (at most) be explained by the (linear) spatial gradient, we fit a model with ONLY the spatial pattern:
RDA_model.arctic.spatial <- rda(arctic.traits.fuzzy ~ Latitude + Longitude,
data = arctic_env)
RDA_model.arctic.spatial## Call: rda(formula = arctic.traits.fuzzy ~ Latitude + Longitude, data =
## arctic_env)
##
## -- Model Summary --
##
## Inertia Proportion Rank
## Total 0.26772 1.00000
## Constrained 0.03438 0.12841 2
## Unconstrained 0.23334 0.87159 17
##
## Inertia is variance
##
## -- Eigenvalues --
##
## Eigenvalues for constrained axes:
## RDA1 RDA2
## 0.028889 0.005487
##
## Eigenvalues for unconstrained axes:
## PC1 PC2 PC3 PC4 PC5 PC6 PC7 PC8
## 0.07614 0.04388 0.03475 0.02103 0.01342 0.01239 0.00858 0.00664
## (Showing 8 of 17 unconstrained eigenvalues)(R2.2 <- RsquareAdj(RDA_model.arctic.spatial)$r.sq)## [1] 0.128405About 12.8% of variability (at most) can be attributed to a linear spatial gradient.
The sum of the \(R^2\) values of the two models above (23.4%) is substantially larger than the variability explained by the full model (15.2%) that includes all variables. The difference is due to multicollinearity; that is, when only the spatial trend is included, it can “absorb” some of the variability otherwise attributed to environmental variables (and vice versa).
Based on these comparisons, we can conclude that environmental covariates may explain up to 10.6% of the variability in trait composition and that environmental covariates, combined with a linear spatial gradient, may account for 16.7% of the variability in trait composition.
Interannual variability / year effects
We can use conditioning to assess if there is additional interannual variability that is not explained by spatial pattern and oceanographic covariates by conditoning on all variables in the full model and adding year as a factor:
table(arctic_env$Year)##
## 2007 2008 2009 2010 2012 2015 2017 2019
## 30 15 84 34 37 65 125 41RDA_model.arctic.year <- rda(arctic.traits.fuzzy ~ factor(Year) +
Condition(Depth + doy_retreat + BottomTemp +
BottomSalinity + Latitude + Longitude),
data = arctic_env)
RDA_model.arctic.year## Call: rda(formula = arctic.traits.fuzzy ~ factor(Year) + Condition(Depth +
## doy_retreat + BottomTemp + BottomSalinity + Latitude + Longitude), data =
## arctic_env)
##
## -- Model Summary --
##
## Inertia Proportion Rank
## Total 0.26772 1.00000
## Conditional 0.04483 0.16746 6
## Constrained 0.02361 0.08819 7
## Unconstrained 0.19928 0.74435 17
##
## Inertia is variance
##
## -- Eigenvalues --
##
## Eigenvalues for constrained axes:
## RDA1 RDA2 RDA3 RDA4 RDA5 RDA6 RDA7
## 0.010084 0.007932 0.002290 0.001454 0.001354 0.000440 0.000057
##
## Eigenvalues for unconstrained axes:
## PC1 PC2 PC3 PC4 PC5 PC6 PC7 PC8
## 0.06661 0.03644 0.02875 0.01636 0.01215 0.01107 0.00784 0.00553
## (Showing 8 of 17 unconstrained eigenvalues)(R2.3 <- RsquareAdj(RDA_model.arctic.year)$r.sq)## [1] 0.08819173anova(RDA_model.arctic.year)## Permutation test for rda under reduced model
## Permutation: free
## Number of permutations: 999
##
## Model: rda(formula = arctic.traits.fuzzy ~ factor(Year) + Condition(Depth + doy_retreat + BottomTemp + BottomSalinity + Latitude + Longitude), data = arctic_env)
## Df Variance F Pr(>F)
## Model 7 0.02361 7.0581 0.001 ***
## Residual 417 0.19928
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1There are additional significant differences among years that are not accounted for by the full model (including the ice retreat variable) and that account for about 8.8% of the variability in trait composition
If a year effect can reasonably be estimated, we can fit a model with all covariates and use a step-wise approach for model selection to identify the reduced model that maximizes the adjusted \(R^2\) value:
RDA_model.arctic.all <- rda(arctic.traits.fuzzy ~ factor(Year) +
doy_retreat + BottomTemp + BottomSalinity +
Depth + Latitude + Longitude,
data = arctic_env)
vif.cca(RDA_model.arctic.all)## factor(Year)2008 factor(Year)2009 factor(Year)2010 factor(Year)2012
## 1.591261 3.854139 2.558639 2.301676
## factor(Year)2015 factor(Year)2017 factor(Year)2019 doy_retreat
## 3.117906 4.737821 2.616883 3.729038
## BottomTemp BottomSalinity Depth Latitude
## 3.343987 2.115824 1.247262 2.577273
## Longitude
## 1.729275RDA_model.arctic.all## Call: rda(formula = arctic.traits.fuzzy ~ factor(Year) + doy_retreat +
## BottomTemp + BottomSalinity + Depth + Latitude + Longitude, data =
## arctic_env)
##
## -- Model Summary --
##
## Inertia Proportion Rank
## Total 0.26772 1.00000
## Constrained 0.06844 0.25565 13
## Unconstrained 0.19928 0.74435 17
##
## Inertia is variance
##
## -- Eigenvalues --
##
## Eigenvalues for constrained axes:
## RDA1 RDA2 RDA3 RDA4 RDA5 RDA6 RDA7 RDA8 RDA9 RDA10
## 0.03961 0.01182 0.00864 0.00252 0.00242 0.00132 0.00110 0.00056 0.00024 0.00014
## RDA11 RDA12 RDA13
## 0.00005 0.00001 0.00000
##
## Eigenvalues for unconstrained axes:
## PC1 PC2 PC3 PC4 PC5 PC6 PC7 PC8
## 0.06661 0.03644 0.02875 0.01636 0.01215 0.01107 0.00784 0.00553
## (Showing 8 of 17 unconstrained eigenvalues)RsquareAdj(RDA_model.arctic.all)## $r.squared
## [1] 0.2556503
##
## $adj.r.squared
## [1] 0.2324452anova(RDA_model.arctic.all, by="margin",
permutations = how(nperm=4999))## Permutation test for rda under reduced model
## Marginal effects of terms
## Permutation: free
## Number of permutations: 4999
##
## Model: rda(formula = arctic.traits.fuzzy ~ factor(Year) + doy_retreat + BottomTemp + BottomSalinity + Depth + Latitude + Longitude, data = arctic_env)
## Df Variance F Pr(>F)
## factor(Year) 7 0.023610 7.0581 0.0002 ***
## doy_retreat 1 0.001707 3.5710 0.0060 **
## BottomTemp 1 0.002261 4.7313 0.0012 **
## BottomSalinity 1 0.000198 0.4143 0.8614
## Depth 1 0.001544 3.2317 0.0124 *
## Latitude 1 0.005170 10.8180 0.0002 ***
## Longitude 1 0.005716 11.9606 0.0002 ***
## Residual 417 0.199275
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1plot(RDA_model.arctic.all)
Including year in the model results in some of the variables being only marginally significant (ice retreat, depth), but the results seem otherwise consistent with the patterns when a year effect is not included.
Step-wise regression
step.best <- ordiR2step(RDA_model.arctic.2, scope = formula(RDA_model.arctic.all),
direction="both")## Step: R2.adj= 0
## Call: arctic.traits.fuzzy ~ 1
##
## R2.adjusted
## <All variables> 0.23244516
## + factor(Year) 0.11001376
## + Latitude 0.10095382
## + BottomTemp 0.08203750
## + Longitude 0.06724984
## + doy_retreat 0.04270619
## + BottomSalinity 0.02466407
## + Depth 0.01524609
## <none> 0.00000000
##
## Df AIC F Pr(>F)
## + factor(Year) 7 -610.29 8.5934 0.002 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: R2.adj= 0.1100138
## Call: arctic.traits.fuzzy ~ factor(Year)
##
## R2.adjusted
## <All variables> 0.2324452
## + Latitude 0.1851526
## + BottomTemp 0.1724627
## + Longitude 0.1608485
## + doy_retreat 0.1445540
## + BottomSalinity 0.1331800
## + Depth 0.1267104
## <none> 0.1100138
##
## Df AIC F Pr(>F)
## + Latitude 1 -647.33 40.006 0.002 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: R2.adj= 0.1851526
## Call: arctic.traits.fuzzy ~ factor(Year) + Latitude
##
## R2.adjusted
## <All variables> 0.2324452
## + Longitude 0.2059012
## + BottomTemp 0.2033283
## + Depth 0.1971948
## + BottomSalinity 0.1945589
## + doy_retreat 0.1885190
## <none> 0.1851526
##
## Df AIC F Pr(>F)
## + Longitude 1 -657.47 12.026 0.002 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: R2.adj= 0.2059012
## Call: arctic.traits.fuzzy ~ factor(Year) + Latitude + Longitude
##
## R2.adjusted
## <All variables> 0.2324452
## + BottomTemp 0.2244299
## + Depth 0.2172672
## + BottomSalinity 0.2145568
## + doy_retreat 0.2113240
## <none> 0.2059012
##
## Df AIC F Pr(>F)
## + BottomTemp 1 -666.67 11.058 0.002 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: R2.adj= 0.2244299
## Call: arctic.traits.fuzzy ~ factor(Year) + Latitude + Longitude + BottomTemp
##
## R2.adjusted
## <All variables> 0.2324452
## + doy_retreat 0.2291436
## + Depth 0.2288895
## <none> 0.2244299
## + BottomSalinity 0.2235796
##
## Df AIC F Pr(>F)
## + doy_retreat 1 -668.32 3.5682 0.012 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: R2.adj= 0.2291436
## Call: arctic.traits.fuzzy ~ factor(Year) + Latitude + Longitude + BottomTemp + doy_retreat
##
## R2.adjusted
## + Depth 0.2335206
## <All variables> 0.2324452
## <none> 0.2291436
## + BottomSalinity 0.2283473step.best## Call: rda(formula = arctic.traits.fuzzy ~ factor(Year) + Latitude +
## Longitude + BottomTemp + doy_retreat, data = arctic_env)
##
## -- Model Summary --
##
## Inertia Proportion Rank
## Total 0.26772 1.00000
## Constrained 0.06662 0.24886 11
## Unconstrained 0.20109 0.75114 17
##
## Inertia is variance
##
## -- Eigenvalues --
##
## Eigenvalues for constrained axes:
## RDA1 RDA2 RDA3 RDA4 RDA5 RDA6 RDA7 RDA8 RDA9 RDA10
## 0.03928 0.01116 0.00845 0.00251 0.00230 0.00126 0.00097 0.00047 0.00017 0.00005
## RDA11
## 0.00001
##
## Eigenvalues for unconstrained axes:
## PC1 PC2 PC3 PC4 PC5 PC6 PC7 PC8
## 0.06708 0.03706 0.02876 0.01637 0.01215 0.01125 0.00797 0.00553
## (Showing 8 of 17 unconstrained eigenvalues)anova(step.best, by="margin")## Permutation test for rda under reduced model
## Marginal effects of terms
## Permutation: free
## Number of permutations: 999
##
## Model: rda(formula = arctic.traits.fuzzy ~ factor(Year) + Latitude + Longitude + BottomTemp + doy_retreat, data = arctic_env)
## Df Variance F Pr(>F)
## factor(Year) 7 0.024281 7.2276 0.001 ***
## Latitude 1 0.005191 10.8163 0.001 ***
## Longitude 1 0.006486 13.5147 0.001 ***
## BottomTemp 1 0.005140 10.7089 0.001 ***
## doy_retreat 1 0.001713 3.5682 0.013 *
## Residual 419 0.201092
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(step.best, by="axis")## Permutation test for rda under reduced model
## Forward tests for axes
## Permutation: free
## Number of permutations: 999
##
## Model: rda(formula = arctic.traits.fuzzy ~ factor(Year) + Latitude + Longitude + BottomTemp + doy_retreat, data = arctic_env)
## Df Variance F Pr(>F)
## RDA1 1 0.039280 81.8457 0.001 ***
## RDA2 1 0.011161 23.2559 0.001 ***
## RDA3 1 0.008450 17.6070 0.001 ***
## RDA4 1 0.002507 5.2235 0.106
## RDA5 1 0.002300 4.7914 0.106
## RDA6 1 0.001255 2.6152 0.581
## RDA7 1 0.000968 2.0171 0.728
## RDA8 1 0.000474 0.9885 0.977
## RDA9 1 0.000168 0.3508 1.000
## RDA10 1 0.000051 0.1062
## RDA11 1 0.000010 0.0198
## Residual 419 0.201092
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1plot(step.best)
(R2.fun.best <- RsquareAdj(step.best)$r.sq)## [1] 0.2488631Best Model for functional traits
The resulting ‘best’ model eliminates salinity and depth and explains approximately 24.9% of overall variability. It clearly suggests that bottom temperature is a major driver of trait composition. In addition, there is a strong remaining spatial gradient, as well as considerable year-to-year variability.
Plot results (ggplot)
Plot ordination of functional traits with individual modalities and environmental variables overlaid.
## RDA1 RDA2
## 14.79 4.41
options(repos = c(CRAN = "https://cran.rstudio.com"))
Inter_south<-(cbind(summary(fuzz.simper, ordered=TRUE)$Intermediate_Southern))
View(Inter_south)
Inter_south<-as.data.frame(cbind(Inter_south,
as.data.frame(rep("Intermediate_Southern",23))))
View(Inter_south)
names(Inter_south)[8]<-"Contrast"
Inter_south$Order <- 1:23
Inter_south_new <- Inter_south %>%
mutate(Contribution = cumsum - lag(cumsum))
View(Inter_south_new)
Inter_south_new[1, 10] <- Inter_south_new[1, 6]
Inter_south_new$Modality<-row.names(Inter_south_new)
Inter_north<-as.data.frame(cbind(summary(fuzz.simper, ordered=TRUE)$Intermediate_Northern))
View(Inter_north)
Inter_north<-(as.data.frame(cbind(Inter_north,
as.data.frame(rep("Intermediate_Northern",23)))))
names(Inter_north)[8]<-"Contrast"
Inter_north$Order<-1:23
View(Inter_north)
Inter_north_new <- Inter_north %>%
mutate(Contribution = cumsum - lag(cumsum))
View(Inter_north_new)
Inter_north_new[1, 10] <- Inter_north_new[1, 6]
Inter_north_new$Modality<-row.names(Inter_north_new)
South_north<-as.data.frame(cbind(summary(fuzz.simper, ordered=TRUE)$Southern_Northern))
South_north<-(as.data.frame(cbind(South_north,
as.data.frame(rep("Southern_Northern",23)))))
names(South_north)[8]<-"Contrast"
South_north$Order<-1:23
View(South_north)
South_north_new <- South_north %>%
mutate(Contribution = cumsum - lag(cumsum))
View(South_north_new)
South_north_new[1, 10] <- South_north_new[1, 6]
South_north_new$Modality<-row.names(South_north_new)
all.contrast<-as.data.frame(rbind(Inter_south_new,
Inter_north_new,
South_north_new))
View(all.contrast)
p1 <- ggplot()+
geom_line(data =all.contrast,
aes(x = Order, y = cumsum,color=Contrast)) +
theme_bw()
p1
p2 <- ggplot()+
geom_line(data =all.contrast,
aes(x = Order, y = average,color=Contrast)) +
theme_bw()
p2
View(all.contrast)
df_summary <- all.contrast %>%
group_by(Modality) %>%
summarise(total_value = sum(Contribution)) %>%
arrange(desc(total_value)) # Arrange in descending order of total_value
View(df_summary)
all.contrast$Modality <- factor(all.contrast$Modality, levels = df_summary$Modality)
mycol## [1] "#FD1C26" "#0D7EF9" "#98A100"sig_label <- as.data.frame(case_when(
all.contrast$p < 0.001 ~ "***",
all.contrast$p < 0.01 ~ "**",
all.contrast$p < 0.05 ~ "*",
TRUE ~ "ns"
))
names(sig_label)[1]<-"Significant"
all.contrast$Significant<-sig_label$Significant
View(all.contrast)
y_coords <- all.contrast %>%
group_by(Modality) %>%
arrange(Modality) %>%
mutate(cumulative_rel_abund = cumsum(Contribution),
center_y = cumulative_rel_abund - (Contribution / 2))
View(y_coords)
# Adjust annotation data to use the centered y-coordinates
all.contrast$Contrast <- factor(all.contrast$Contrast, levels = c("Southern_Northern","Intermediate_Northern", "Intermediate_Southern"))
p3<-ggplot() +
geom_bar(data=all.contrast, aes(x = Modality, y = Contribution,
fill = Contrast),
stat = "identity") +
scale_fill_manual(values=c("#4D1FF9","#FD6C29", "#98E111"))+
geom_text(data=y_coords[y_coords$Significant=="*",],
aes(x = Modality, label = Significant, y = center_y),
vjust = 0.5, size = 6)+
theme_bw() +
theme(axis.text.x = element_text(angle = 45, hjust = 1, color="black", size = 12),
axis.text.y = element_text(color="black"),
axis.title = element_text(color="black",size=12),
legend.position = c(.8,.8),
legend.text = element_text(size=12),
legend.title = element_text(size=12),
plot.margin = margin(t = 1, r = 1, b = 1, l = 1, unit = "cm"))## Warning: A numeric `legend.position` argument in `theme()` was deprecated in ggplot2
## 3.5.0.
## ℹ Please use the `legend.position.inside` argument of `theme()` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.p3
ggarrange(p1, p2, p3,common.legend = TRUE)
View(Inter_south)
View(Inter_north)
View(South_north)
install.packages("writexl")##
## The downloaded binary packages are in
## /var/folders/62/l2kf23qd3qj5pxmqk6jprx840000gq/T//Rtmp6vCG4Z/downloaded_packageslibrary(writexl)## Warning: package 'writexl' was built under R version 4.3.3getwd()## [1] "/Users/lauren"write_xlsx(all.contrast, "/Users/lauren/Desktop/arctic fish/all.contrast.traits.xlsx")Taxonomic clustering
- Data prep & dendrogram
arctic_cpue <- as.data.frame(arctic_cpue_090424.2) %>%
filter(row.names(arctic_cpue_090424.2) %in%
row.names(arctic_env_090424))
# Bray-Curtis distances:
tax.dis <- vegdist(arctic_cpue, method="bray")
tax.clus <- hclust(tax.dis, method = "ward.D2")
plot(tax.clus, labels=F)
dend_tax <- as.dendrogram(tax.clus)
plot(dend_tax)
# Cut the tree to get,
groups_tax <- as.data.frame(cutree(dend_tax, k = 3))
View(groups_tax)
names(groups_tax)[1]<-"cluster"
groups_tax$cluster<-as.factor(groups_tax$cluster)
summary(groups_tax)## cluster
## 1:108
## 2:204
## 3:119mycol_tax<-c("#0D7EF9","#FD1C26", "#98A100")
dend_colored_branches <- color_branches(dend, k = 3,
groupLabels = FALSE,col=mycol)
dend_colored_branches_tax <- color_branches(dend_tax, k = 3,
groupLabels = FALSE,col=mycol_tax)
plot(dend_colored_branches_tax)
dend_colored_branches_tax<-set(dend_colored_branches_tax,"labels_color", 'white')
plot(dend_colored_branches_tax)
There are three (or 5) clusters. We cut this off at three
Map clusters
## [1] TRUE## Coordinate system already present. Adding new coordinate system, which will
## replace the existing one.
## Coordinate system already present. Adding new coordinate system, which will
## replace the existing one.
Stacked barplot of species compositions:
###simper
## simper species
tax.simper<-with(arctic_env_tax,
simper(arctic_cpue,
tax.cluster, permutations = 99))
tax.simper## cumulative contributions of most influential species:
##
## $Nearshore_Offshore
## Myoxocephalus.scorpius Anisarchus.medius Boreogadus.saida
## 0.06785586 0.13400120 0.20007698
## Gymnocanthus.tricuspis Hippoglossoides.robustus Artediellus.scaber
## 0.26478727 0.32893804 0.38964731
## Lumpenus.fabricii.sp. Lycodes.spp. Aspidophoroides.olrikii
## 0.44613262 0.50247270 0.55236165
## Liparis.spp. Stichaeus.punctatus Eleginus.gracilis
## 0.59748240 0.63860869 0.67477736
## Liparis.tunicatus
## 0.70750400
##
## $Nearshore_Northern
## Gymnocanthus.tricuspis Myoxocephalus.scorpius Lumpenus.fabricii.sp.
## 0.08464982 0.16748266 0.24286849
## Anisarchus.medius Lycodes.spp. Boreogadus.saida
## 0.31780052 0.38623794 0.44461352
## Artediellus.scaber Aspidophoroides.olrikii Liparis.spp.
## 0.50044838 0.54648277 0.58830494
## Stichaeus.punctatus Hippoglossoides.robustus Gymnelus.hemifasciatus
## 0.62516177 0.66007526 0.69341652
## Liparis.tunicatus
## 0.72505842
##
## $Offshore_Northern
## Lumpenus.fabricii.sp. Gymnocanthus.tricuspis Hippoglossoides.robustus
## 0.1051699 0.1875915 0.2638488
## Anisarchus.medius Myoxocephalus.scorpius Lycodes.spp.
## 0.3338862 0.4024152 0.4667623
## Aspidophoroides.olrikii Boreogadus.saida Artediellus.scaber
## 0.5266033 0.5859450 0.6367522
## Liparis.spp. Gymnelus.hemifasciatus
## 0.6854832 0.7288322Near_Off<-as.data.frame(cbind(summary(tax.simper, ordered=TRUE)$Nearshore_Offshore))
View(Near_Off)
Near_Off<-as.data.frame(cbind(Near_Off,
as.data.frame(rep("Nearshore_Offshore",33))))
View(Near_Off)
names(Near_Off)[8]<-"Contrast"
Near_Off$Order <- 1:33
Near_Off_new <- Near_Off %>%
mutate(Contribution = cumsum - lag(cumsum))
View(Near_Off_new)
Near_Off_new[1, 10] <- Near_Off_new[1, 6]
Near_Off_new$Taxa<-row.names(Near_Off_new)
Near_North<-as.data.frame(cbind(summary(tax.simper, ordered=TRUE)$Nearshore_Northern))
View(Near_North)
Near_North<-(as.data.frame(cbind(Near_North,
as.data.frame(rep("Nearshore_Northern",33)))))
names(Near_North)[8]<-"Contrast"
Near_North$Order<-1:33
View(Near_North)
Near_North_new <- Near_North %>%
mutate(Contribution = cumsum - lag(cumsum))
View(Near_North_new)
Near_North_new[1, 10] <- Near_North_new[1, 6]
Near_North_new$Taxa<-row.names(Near_North_new)
Off_North<-as.data.frame(cbind(summary(tax.simper, ordered=TRUE)$Offshore_Northern))
Off_North<-(as.data.frame(cbind(Off_North,
as.data.frame(rep("Offshore_Northern",33)))))
names(Off_North)[8]<-"Contrast"
Off_North$Order<-1:33
View(Off_North)
Off_North_new <- Off_North %>%
mutate(Contribution = cumsum - lag(cumsum))
View(Off_North_new)
Off_North_new[1, 10] <- Off_North_new[1, 6]
Off_North_new$Taxa<-row.names(Off_North_new)
all.taxa.contrast<-as.data.frame(rbind(Near_Off_new,
Near_North_new,
Off_North_new))
View(all.taxa.contrast)
p1.tax <- ggplot()+
geom_line(data =all.taxa.contrast,
aes(x = Order, y = cumsum,color=Contrast)) +
theme_bw()
p1.tax
p2.tax <- ggplot()+
geom_line(data =all.taxa.contrast,
aes(x = Order, y = average,color=Contrast)) +
theme_bw()
p2.tax
View(all.contrast)
df_summary_tax <- all.taxa.contrast %>%
group_by(Taxa) %>%
summarise(total_value = sum(Contribution)) %>%
arrange(desc(total_value)) # Arrange in descending order of total_value
View(df_summary_tax)
all.taxa.contrast$Taxa <- factor(all.taxa.contrast$Taxa,
levels = df_summary_tax$Taxa)
mycol## [1] "#FD1C26" "#0D7EF9" "#98A100"View(all.taxa.contrast)
sig_label_tax <- as.data.frame(case_when(
all.taxa.contrast$p < 0.001 ~ "***",
all.taxa.contrast$p < 0.01 ~ "**",
all.taxa.contrast$p < 0.05 ~ "*",
TRUE ~ "ns"
))
names(sig_label_tax)[1]<-"Significant"
all.taxa.contrast$Significant<-sig_label_tax$Significant
View(all.taxa.contrast)
y_coords_taxa <- all.taxa.contrast %>%
group_by(Taxa) %>%
arrange(Taxa) %>%
mutate(cumulative_rel_abund = cumsum(Contribution),
center_y = cumulative_rel_abund - (Contribution / 2))
View(y_coords_taxa)
# Adjust annotation data to use the centered y-coordinates
View(all.taxa.contrast)
all.taxa.contrast$Contrast <-
factor(all.taxa.contrast$Contrast, levels = c("Offshore_Northern","Nearshore_Northern", "Nearshore_Offshore"))
p3.taxa<-ggplot() +
geom_bar(data=all.taxa.contrast, aes(x = Taxa, y = Contribution,
fill = Contrast),
stat = "identity") +
scale_fill_manual(values=c("#4D1FF9", "#98E111","#FD6C29"))+
geom_text(data=y_coords_taxa[y_coords_taxa$Significant=="*",],
aes(x = Taxa, label = Significant, y = center_y),
vjust = 0.7, size = 6)+
theme_bw() +
theme(axis.text.x = element_text(angle = 45, hjust = 1, color="black", size = 12),
axis.text.y = element_text(color="black"),
axis.title = element_text(color="black",size=12),
legend.position = c(.8,.8),
legend.text = element_text(size=12),
legend.title = element_text(size=12),
plot.margin = margin(t = 1.5, r = 1.5, b = 1.5, l = 1.5, unit = "cm"))
p3.taxa
ggarrange(p1, p2, p3,common.legend = TRUE)
install.packages("writexl")##
## The downloaded binary packages are in
## /var/folders/62/l2kf23qd3qj5pxmqk6jprx840000gq/T//Rtmp6vCG4Z/downloaded_packageslibrary(writexl)
getwd()## [1] "/Users/lauren"write_xlsx(all.taxa.contrast, "/Users/lauren/Desktop/arctic fish/all.contrast.taxa.xlsx")Ordination based on taxa
CCA_model.arctic.tax <- cca(arctic_cpue ~
factor(Year) + Depth + doy_retreat + BottomTemp +
BottomSalinity + Latitude + Longitude,
data = arctic_env)
vif.cca(CCA_model.arctic.tax)## factor(Year)2008 factor(Year)2009 factor(Year)2010 factor(Year)2012
## 1.604783 3.314902 1.849274 2.135118
## factor(Year)2015 factor(Year)2017 factor(Year)2019 Depth
## 2.484255 3.941739 2.137603 1.373861
## doy_retreat BottomTemp BottomSalinity Latitude
## 3.527987 3.586086 2.558954 2.452021
## Longitude
## 1.807384CCA_model.arctic.tax## Call: cca(formula = arctic_cpue ~ factor(Year) + Depth + doy_retreat +
## BottomTemp + BottomSalinity + Latitude + Longitude, data = arctic_env)
##
## -- Model Summary --
##
## Inertia Proportion Rank
## Total 2.6272 1.0000
## Constrained 0.4912 0.1870 13
## Unconstrained 2.1360 0.8130 32
##
## Inertia is scaled Chi-square
##
## -- Eigenvalues --
##
## Eigenvalues for constrained axes:
## CCA1 CCA2 CCA3 CCA4 CCA5 CCA6 CCA7 CCA8 CCA9 CCA10
## 0.17185 0.08012 0.07115 0.04144 0.03229 0.02698 0.01985 0.01472 0.01289 0.00748
## CCA11 CCA12 CCA13
## 0.00668 0.00316 0.00257
##
## Eigenvalues for unconstrained axes:
## CA1 CA2 CA3 CA4 CA5 CA6 CA7 CA8
## 0.16433 0.12777 0.11988 0.10552 0.10031 0.09188 0.08794 0.08743
## (Showing 8 of 32 unconstrained eigenvalues)(R2.4 <- RsquareAdj(CCA_model.arctic.tax)$r.sq)##0.187## [1] 0.1869642plot(CCA_model.arctic.tax)
Similar to the functional traits CWM-RDA analysis, the results show a gradient from warm, low salinity, shallow stations with an early ice retreat to the opposite, and a spatial gradient from SW to NE (Bering Strait to NE Chukchi Sea).
Assess significance of model overall, CCA axes and individual terms (marginal effects)
anova(CCA_model.arctic.tax)## Permutation test for cca under reduced model
## Permutation: free
## Number of permutations: 999
##
## Model: cca(formula = arctic_cpue ~ factor(Year) + Depth + doy_retreat + BottomTemp + BottomSalinity + Latitude + Longitude, data = arctic_env)
## Df ChiSquare F Pr(>F)
## Model 13 0.49118 7.3764 0.001 ***
## Residual 417 2.13597
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(CCA_model.arctic.tax, by="axis") ##CCA1-7 significant## Permutation test for cca under reduced model
## Forward tests for axes
## Permutation: free
## Number of permutations: 999
##
## Model: cca(formula = arctic_cpue ~ factor(Year) + Depth + doy_retreat + BottomTemp + BottomSalinity + Latitude + Longitude, data = arctic_env)
## Df ChiSquare F Pr(>F)
## CCA1 1 0.17185 33.5494 0.001 ***
## CCA2 1 0.08012 15.6411 0.001 ***
## CCA3 1 0.07115 13.8904 0.001 ***
## CCA4 1 0.04144 8.0912 0.001 ***
## CCA5 1 0.03229 6.3043 0.001 ***
## CCA6 1 0.02698 5.2672 0.001 ***
## CCA7 1 0.01985 3.8759 0.008 **
## CCA8 1 0.01472 2.8747 0.310
## CCA9 1 0.01289 2.5166 0.684
## CCA10 1 0.00748 1.4600 1.000
## CCA11 1 0.00668 1.3033
## CCA12 1 0.00316 0.6175
## CCA13 1 0.00257 0.5010
## Residual 417 2.13597
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(CCA_model.arctic.tax, by="margin")## Permutation test for cca under reduced model
## Marginal effects of terms
## Permutation: free
## Number of permutations: 999
##
## Model: cca(formula = arctic_cpue ~ factor(Year) + Depth + doy_retreat + BottomTemp + BottomSalinity + Latitude + Longitude, data = arctic_env)
## Df ChiSquare F Pr(>F)
## factor(Year) 7 0.20697 5.7724 0.001 ***
## Depth 1 0.03094 6.0400 0.001 ***
## doy_retreat 1 0.01194 2.3310 0.001 ***
## BottomTemp 1 0.02124 4.1463 0.001 ***
## BottomSalinity 1 0.00897 1.7513 0.030 *
## Latitude 1 0.03760 7.3404 0.001 ***
## Longitude 1 0.02461 4.8036 0.001 ***
## Residual 417 2.13597
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1The overall model is significant and accounts for approximately 18.7% of the variability in species composition.
The first 5 CCA axes capture significant proportions of overall variability.
Step-wise regression:
CCA_model.arctic.2.tax <- cca( arctic_cpue ~ 1 ,
data = arctic_env)
step.res.tax<-ordiR2step(CCA_model.arctic.2.tax,
scope = formula(CCA_model.arctic.tax),
direction="both")## Step: R2.adj= 0
## Call: arctic_cpue ~ 1
##
## R2.adjusted
## <All variables> 0.16055523
## + factor(Year) 0.07682357
## + BottomTemp 0.04457977
## + Latitude 0.03964139
## + BottomSalinity 0.02991904
## + Depth 0.02570483
## + Longitude 0.01876224
## + doy_retreat 0.01747507
## <none> 0.00000000
##
## Df AIC F Pr(>F)
## + factor(Year) 7 1995.7 6.1574 0.002 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: R2.adj= 0.0768796
## Call: arctic_cpue ~ factor(Year)
##
## R2.adjusted
## <All variables> 0.16055523
## + BottomTemp 0.11982109
## + Latitude 0.11149302
## + BottomSalinity 0.10663844
## + Depth 0.10354539
## + doy_retreat 0.09176937
## + Longitude 0.09127404
## <none> 0.07687960
##
## Df AIC F Pr(>F)
## + BottomTemp 1 1976.1 21.699 0.002 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: R2.adj= 0.1198385
## Call: arctic_cpue ~ factor(Year) + BottomTemp
##
## R2.adjusted
## <All variables> 0.1605552
## + Latitude 0.1335469
## + Longitude 0.1319704
## + Depth 0.1315185
## + BottomSalinity 0.1260044
## + doy_retreat 0.1248217
## <none> 0.1198385
##
## Df AIC F Pr(>F)
## + Latitude 1 1970.3 7.711 0.002 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: R2.adj= 0.133541
## Call: arctic_cpue ~ factor(Year) + BottomTemp + Latitude
##
## R2.adjusted
## <All variables> 0.1605552
## + Depth 0.1468750
## + Longitude 0.1446483
## + BottomSalinity 0.1398926
## + doy_retreat 0.1375133
## <none> 0.1335410
##
## Df AIC F Pr(>F)
## + Depth 1 1964.5 7.6592 0.002 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: R2.adj= 0.1469882
## Call: arctic_cpue ~ factor(Year) + BottomTemp + Latitude + Depth
##
## R2.adjusted
## <All variables> 0.1605552
## + Longitude 0.1564859
## + BottomSalinity 0.1504629
## + doy_retreat 0.1500216
## <none> 0.1469882
##
## Df AIC F Pr(>F)
## + Longitude 1 1960.7 5.7389 0.002 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: R2.adj= 0.1564748
## Call: arctic_cpue ~ factor(Year) + BottomTemp + Latitude + Depth + Longitude
##
## R2.adjusted
## <All variables> 0.1605552
## + doy_retreat 0.1592009
## + BottomSalinity 0.1579555
## <none> 0.1564748
##
## Df AIC F Pr(>F)
## + doy_retreat 1 1960.2 2.4277 0.002 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: R2.adj= 0.1591182
## Call: arctic_cpue ~ factor(Year) + BottomTemp + Latitude + Depth + Longitude + doy_retreat
##
## R2.adjusted
## + BottomSalinity 0.1607089
## <All variables> 0.1605552
## <none> 0.1591182step.res.tax## Call: cca(formula = arctic_cpue ~ factor(Year) + BottomTemp + Latitude +
## Depth + Longitude + doy_retreat, data = arctic_env)
##
## -- Model Summary --
##
## Inertia Proportion Rank
## Total 2.6272 1.0000
## Constrained 0.4822 0.1835 12
## Unconstrained 2.1449 0.8165 32
##
## Inertia is scaled Chi-square
##
## -- Eigenvalues --
##
## Eigenvalues for constrained axes:
## CCA1 CCA2 CCA3 CCA4 CCA5 CCA6 CCA7 CCA8 CCA9 CCA10
## 0.17173 0.07929 0.07057 0.04055 0.03226 0.02691 0.01982 0.01441 0.01072 0.00720
## CCA11 CCA12
## 0.00585 0.00290
##
## Eigenvalues for unconstrained axes:
## CA1 CA2 CA3 CA4 CA5 CA6 CA7 CA8
## 0.16490 0.12912 0.12029 0.10552 0.10071 0.09236 0.08915 0.08776
## (Showing 8 of 32 unconstrained eigenvalues)anova(step.res.tax, by="margin")## Permutation test for cca under reduced model
## Marginal effects of terms
## Permutation: free
## Number of permutations: 999
##
## Model: cca(formula = arctic_cpue ~ factor(Year) + BottomTemp + Latitude + Depth + Longitude + doy_retreat, data = arctic_env)
## Df ChiSquare F Pr(>F)
## factor(Year) 7 0.21014 5.8502 0.001 ***
## BottomTemp 1 0.03439 6.7021 0.001 ***
## Latitude 1 0.03803 7.4118 0.001 ***
## Depth 1 0.03435 6.6935 0.001 ***
## Longitude 1 0.02870 5.5930 0.001 ***
## doy_retreat 1 0.01246 2.4277 0.001 ***
## Residual 418 2.14494
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1(R2.best <- RsquareAdj(step.res.tax)$r.sq)## [1] 0.1835497plot(step.res.tax)
Best Model for species composition
The ‘best’ model with the highest adjusted \(R^2\) eliminates salinity only but retains all other variables and all of the remaining variables are highly significant. It explains slightly less of the overall variability in species composition than the full model (18.4%)
Environment only model
We fit a model without the spatial pattern to assess how much of the variability can (at most) be explained by depth and dynamic ocean variables:
CCA_model.arctic.env.tax <- cca(arctic_cpue ~
Depth + doy_retreat + BottomTemp +
BottomSalinity,
data = arctic_env)
CCA_model.arctic.env.tax## Call: cca(formula = arctic_cpue ~ Depth + doy_retreat + BottomTemp +
## BottomSalinity, data = arctic_env)
##
## -- Model Summary --
##
## Inertia Proportion Rank
## Total 2.62715 1.00000
## Constrained 0.20527 0.07813 4
## Unconstrained 2.42188 0.92187 32
##
## Inertia is scaled Chi-square
##
## -- Eigenvalues --
##
## Eigenvalues for constrained axes:
## CCA1 CCA2 CCA3 CCA4
## 0.13664 0.04110 0.01552 0.01201
##
## Eigenvalues for unconstrained axes:
## CA1 CA2 CA3 CA4 CA5 CA6 CA7 CA8
## 0.18807 0.14467 0.14086 0.12835 0.11414 0.11089 0.10317 0.10123
## (Showing 8 of 32 unconstrained eigenvalues)(R2.5 <- RsquareAdj(CCA_model.arctic.env.tax)$r.sq) ##0.078## [1] 0.07813255About 7.8% of variability (at most) can be accounted for by bathymetry and dynamic ocean variables
Spatial model
Next, we fit a model with ONLY the spatial pattern to assess how much of the variability can (at most) be explained by the (linear) spatial gradient:
CCA_model.arctic.spatial.tax <- cca(arctic_cpue ~ Latitude + Longitude,
data = arctic_env)
CCA_model.arctic.spatial.tax## Call: cca(formula = arctic_cpue ~ Latitude + Longitude, data = arctic_env)
##
## -- Model Summary --
##
## Inertia Proportion Rank
## Total 2.6271 1.0000
## Constrained 0.1497 0.0570 2
## Unconstrained 2.4774 0.9430 32
##
## Inertia is scaled Chi-square
##
## -- Eigenvalues --
##
## Eigenvalues for constrained axes:
## CCA1 CCA2
## 0.11098 0.03877
##
## Eigenvalues for unconstrained axes:
## CA1 CA2 CA3 CA4 CA5 CA6 CA7 CA8
## 0.18696 0.17954 0.14929 0.12910 0.12519 0.11219 0.10270 0.10017
## (Showing 8 of 32 unconstrained eigenvalues)(R2.6 <- RsquareAdj(CCA_model.arctic.spatial.tax)$r.sq) ##0.057## [1] 0.05700268Only about 5.7% of variability (at most) can be attributed to a linear spatial gradient.
Given that the overall model (including year) explains about 18.7% of the variability, there appears to be considerable year-to-year variability not explained by any of the environmental or spatial gradients.
The additional interannual variability that is not explained by spatial pattern and oceanographic covariates can be quantified by conditoning on all other variables in the best model:
View(arctic_env)
step.res.tax2 <- cca(arctic_cpue ~ factor(Year) + Condition(BottomTemp +
Latitude + Depth + doy_retreat + Longitude + BottomSalinity),
data = arctic_env)
step.res.tax2## Call: cca(formula = arctic_cpue ~ factor(Year) + Condition(BottomTemp +
## Latitude + Depth + doy_retreat + Longitude + BottomSalinity), data =
## arctic_env)
##
## -- Model Summary --
##
## Inertia Proportion Rank
## Total 2.62715 1.00000
## Conditional 0.28421 0.10818 6
## Constrained 0.20697 0.07878 7
## Unconstrained 2.13597 0.81304 32
##
## Inertia is scaled Chi-square
##
## -- Eigenvalues --
##
## Eigenvalues for constrained axes:
## CCA1 CCA2 CCA3 CCA4 CCA5 CCA6 CCA7
## 0.07052 0.04727 0.03465 0.02537 0.01714 0.00625 0.00578
##
## Eigenvalues for unconstrained axes:
## CA1 CA2 CA3 CA4 CA5 CA6 CA7 CA8
## 0.16433 0.12777 0.11988 0.10552 0.10031 0.09188 0.08794 0.08743
## (Showing 8 of 32 unconstrained eigenvalues)anova(step.res.tax2)## Permutation test for cca under reduced model
## Permutation: free
## Number of permutations: 999
##
## Model: cca(formula = arctic_cpue ~ factor(Year) + Condition(BottomTemp + Latitude + Depth + doy_retreat + Longitude + BottomSalinity), data = arctic_env)
## Df ChiSquare F Pr(>F)
## Model 7 0.20697 5.7724 0.001 ***
## Residual 417 2.13597
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1(R2.7 <- RsquareAdj(step.res.tax2)$r.sq) ##0.078## [1] 0.07878249There are additional significant differences among years that are not accounted for by other variables in the best model that account for about 7.9% of the variability in taxonomic composition
Plot results (ggplot)
Plot CCA ordination with species and environmental variables overlaid. (We usde different scaling for environmental vectors (x4) from default CCA plot for better visualization of taxa names)
## [1] "Depth" "doy_retreat" "BottomTemp" "BottomSalinity"
## [5] "Latitude" "Longitude"## CCA1 CCA2
## 6.54 3.05##
## The downloaded binary packages are in
## /var/folders/62/l2kf23qd3qj5pxmqk6jprx840000gq/T//Rtmp6vCG4Z/downloaded_packages## Registered S3 methods overwritten by 'ggpp':
## method from
## heightDetails.titleGrob ggplot2
## widthDetails.titleGrob ggplot2##
## Attaching package: 'ggpp'## The following objects are masked from 'package:ggpubr':
##
## as_npc, as_npcx, as_npcy## The following object is masked from 'package:ggplot2':
##
## annotate
Look at diversity metrics (RaoQ and Inverse Simpson’s diversity)
compared to clusters
fd.chukchi<-dbFD(arctic_traits_090424_b,
arctic_cpue_090424.2,
w.abun = TRUE,corr = "none",
calc.FRic = TRUE)## FEVe: Could not be calculated for communities with <3 functionally singular species.
## FRic: To respect s > t, FRic could not be calculated for communities with <3 functionally singular species.
## FRic: Dimensionality reduction was required. The last 20 PCoA axes (out of 22 in total) were removed.
## FRic: Quality of the reduced-space representation = 0.3694546
## FDiv: Could not be calculated for communities with <3 functionally singular species.View(as.data.frame(fd.chukchi$FDis))
arctic_env_fd<-arctic_env
summary(arctic_env_fd)## Year Latitude Longitude doy_retreat
## Min. :2007 Min. :65.65 Min. :-169.6 Min. :119.0
## 1st Qu.:2009 1st Qu.:69.50 1st Qu.:-166.9 1st Qu.:143.0
## Median :2015 Median :70.90 Median :-165.0 Median :161.0
## Mean :2013 Mean :70.35 Mean :-164.8 Mean :163.6
## 3rd Qu.:2017 3rd Qu.:71.39 3rd Qu.:-163.0 3rd Qu.:183.0
## Max. :2019 Max. :73.00 Max. :-154.4 Max. :233.0
## Depth BottomTemp BottomSalinity cluster
## Min. :15.00 Min. :-1.760 Min. :27.15 Southern :139
## 1st Qu.:39.90 1st Qu.: 0.328 1st Qu.:32.06 Intermediate:225
## Median :43.00 Median : 2.140 Median :32.28 Northern : 67
## Mean :42.87 Mean : 2.600 Mean :32.18
## 3rd Qu.:47.00 3rd Qu.: 4.439 3rd Qu.:32.50
## Max. :90.00 Max. :10.660 Max. :34.56arctic_env_fd$raoQ<-as.numeric(fd.chukchi$RaoQ)
simpson<-as.data.frame(diversity(arctic_cpue, "simpson"))
names(simpson)<-"simpson"
arctic_env_tax$simpson<-simpson$simpson
?diversity## Help on topic 'diversity' was found in the following packages:
##
## Package Library
## igraph /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library
## vegan /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library
##
##
## Using the first match ...View(arctic_env_tax)
library(glmm)## Warning: package 'glmm' was built under R version 4.3.3## Loading required package: trust## Loading required package: mvtnorm## Loading required package: Matrix##
## Attaching package: 'Matrix'## The following objects are masked from 'package:tidyr':
##
## expand, pack, unpack## Loading required package: parallel## Loading required package: doParallel## Loading required package: foreach##
## Attaching package: 'foreach'## The following objects are masked from 'package:purrr':
##
## accumulate, when## Loading required package: iteratorslibrary(mgcv)## Loading required package: nlme##
## Attaching package: 'nlme'## The following object is masked from 'package:dplyr':
##
## collapse## This is mgcv 1.8-42. For overview type 'help("mgcv-package")'.lm.fd<-lm(raoQ~ cluster,
data=arctic_env_fd)
lm.fd.1<-lm(raoQ~1,
data=arctic_env_fd)
anova(lm.fd, lm.fd.1)## Analysis of Variance Table
##
## Model 1: raoQ ~ cluster
## Model 2: raoQ ~ 1
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 428 2470.6
## 2 430 3181.0 -2 -710.43 61.536 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1fd.aov<-aov(lm.fd)
summary(lm.fd)##
## Call:
## lm(formula = raoQ ~ cluster, data = arctic_env_fd)
##
## Residuals:
## Min 1Q Median 3Q Max
## -10.768 -1.457 0.317 1.635 6.244
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 20.8016 0.2038 102.076 <2e-16 ***
## clusterIntermediate -0.5812 0.2592 -2.242 0.0255 *
## clusterNorthern -3.8295 0.3573 -10.717 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.403 on 428 degrees of freedom
## Multiple R-squared: 0.2233, Adjusted R-squared: 0.2197
## F-statistic: 61.54 on 2 and 428 DF, p-value: < 2.2e-16TukeyHSD(fd.aov)## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = lm.fd)
##
## $cluster
## diff lwr upr p adj
## Intermediate-Southern -0.5811887 -1.190796 0.02841858 0.065407
## Northern-Southern -3.8294549 -4.669857 -2.98905320 0.000000
## Northern-Intermediate -3.2482662 -4.034698 -2.46183473 0.000000library(rcompanion)## Warning: package 'rcompanion' was built under R version 4.3.3Sum = groupwiseMean(raoQ ~ cluster,
data = arctic_env_fd,
conf = 0.95,
digits = 3)
clus.cols.fuzzy <- c("#FD1C26","#98A100", "#0D7EF9","turquoise")
fuzzy.cluster.CI<-
ggplot(data = Sum,
aes(x=cluster,y= Mean, color=cluster)) +
geom_errorbar(aes(ymin=Trad.lower,
ymax=Trad.upper,
width=0.15))+
geom_point(shape = 15,
size = 4) +
scale_color_manual(values=clus.cols.fuzzy)+
theme_bw()+
labs(x="Cluster", y="Rao's Quadratic Entropy")
fuzzy.cluster.CI
##taxonomy
View(arctic_env_fd)
lm.tax<-lm(simpson~ tax.cluster,
data=arctic_env_tax)
summary(lm.tax)##
## Call:
## lm(formula = simpson ~ tax.cluster, data = arctic_env_tax)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.36048 -0.02739 0.01854 0.04493 0.09430
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.8595921 0.0045832 187.554 < 2e-16 ***
## tax.clusterNearshore 0.0008893 0.0077899 0.114 0.909
## tax.clusterNorthern -0.0385474 0.0075508 -5.105 4.99e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.06546 on 428 degrees of freedom
## Multiple R-squared: 0.06624, Adjusted R-squared: 0.06188
## F-statistic: 15.18 on 2 and 428 DF, p-value: 4.267e-07tax.aov<-aov(lm.tax)
summary(tax.aov)## Df Sum Sq Mean Sq F value Pr(>F)
## tax.cluster 2 0.1301 0.06505 15.18 4.27e-07 ***
## Residuals 428 1.8340 0.00429
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1TukeyHSD(tax.aov)## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = lm.tax)
##
## $tax.cluster
## diff lwr upr p adj
## Nearshore-Offshore 0.000889308 -0.01743170 0.01921032 0.9928406
## Northern-Offshore -0.038547367 -0.05630609 -0.02078864 0.0000015
## Northern-Nearshore -0.039436675 -0.05989766 -0.01897569 0.0000225Sum.simpson = groupwiseMean(simpson ~ tax.cluster,
data = arctic_env_tax,
conf = 0.95,
digits = 3)
clus.cols.tax <- c("#FD1C26","#98A100", "#0D7EF9")
Sum.simpson$tax.cluster <- factor(Sum.simpson$tax.cluster ,
levels = c("Nearshore",
"Offshore",
"Northern"))
tax.cluster.CI<-
ggplot(data = Sum.simpson,
aes(x=tax.cluster,y= Mean, color=tax.cluster)) +
geom_errorbar(aes(ymin=Trad.lower,
ymax=Trad.upper,
width=0.15))+
geom_point(shape = 15,
size = 4) +
scale_color_manual(values=clus.cols.tax)+
theme_bw()+
labs(x="Cluster", y="Simpson Diversity Index")
tax.cluster.CI
Diversity metrics compared to environmental gradients for functional
diversity
View(arctic_env)
arctic.cwm.rda.sites$raoQ<-as.numeric(arctic_env_fd$raoQ)
###using linear models
install.packages("AICcmodavg")##
## The downloaded binary packages are in
## /var/folders/62/l2kf23qd3qj5pxmqk6jprx840000gq/T//Rtmp6vCG4Z/downloaded_packageslibrary(AICcmodavg)## Warning: package 'AICcmodavg' was built under R version 4.3.3lm.rda.0<-lm(raoQ~ 1,
data=arctic.cwm.rda.sites)
lm.rda.1<-lm(raoQ~ RDA1 + I(RDA1^2) + RDA2 + I(RDA2^2) ,
data=arctic.cwm.rda.sites)
anova(lm.rda.0, lm.rda.1)##r2 = 0.277 p<0.001## Analysis of Variance Table
##
## Model 1: raoQ ~ 1
## Model 2: raoQ ~ RDA1 + I(RDA1^2) + RDA2 + I(RDA2^2)
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 430 3181.0
## 2 426 2278.4 4 902.68 42.195 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(lm.rda.1)##
## Call:
## lm(formula = raoQ ~ RDA1 + I(RDA1^2) + RDA2 + I(RDA2^2), data = arctic.cwm.rda.sites)
##
## Residuals:
## Min 1Q Median 3Q Max
## -9.6797 -1.4948 0.2389 1.4690 8.2247
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 20.979 0.147 142.681 < 2e-16 ***
## RDA1 -6.880 1.509 -4.559 6.72e-06 ***
## I(RDA1^2) -112.924 13.209 -8.549 2.25e-16 ***
## RDA2 -7.944 1.218 -6.523 1.95e-10 ***
## I(RDA2^2) -49.728 7.881 -6.310 7.00e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.313 on 426 degrees of freedom
## Multiple R-squared: 0.2838, Adjusted R-squared: 0.277
## F-statistic: 42.2 on 4 and 426 DF, p-value: < 2.2e-16AIC(lm.rda.0, lm.rda.1)## df AIC
## lm.rda.0 2 2088.632
## lm.rda.1 6 1952.784lm.rda.2<-lm(raoQ~ RDA1 + I(RDA1^2), data=arctic.cwm.rda.sites)
anova(lm.rda.1, lm.rda.2)## Analysis of Variance Table
##
## Model 1: raoQ ~ RDA1 + I(RDA1^2) + RDA2 + I(RDA2^2)
## Model 2: raoQ ~ RDA1 + I(RDA1^2)
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 426 2278.4
## 2 428 2676.3 -2 -397.99 37.207 1.279e-15 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(lm.rda.2) ##r2 = 0.1547, p<0.001##
## Call:
## lm(formula = raoQ ~ RDA1 + I(RDA1^2), data = arctic.cwm.rda.sites)
##
## Residuals:
## Min 1Q Median 3Q Max
## -11.2192 -1.5152 0.3064 1.7799 5.9653
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 20.5158 0.1441 142.381 < 2e-16 ***
## RDA1 -5.7515 1.6117 -3.569 4e-04 ***
## I(RDA1^2) -108.1229 13.9503 -7.751 6.72e-14 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.501 on 428 degrees of freedom
## Multiple R-squared: 0.1587, Adjusted R-squared: 0.1547
## F-statistic: 40.36 on 2 and 428 DF, p-value: < 2.2e-16AIC(lm.rda.1, lm.rda.2) ##lm.rda.1## df AIC
## lm.rda.1 6 1952.784
## lm.rda.2 4 2018.174lm.rda.3<-lm(raoQ~ RDA1 + RDA2 ,
data=arctic.cwm.rda.sites)
anova(lm.rda.1, lm.rda.3) ##lm.rda.1## Analysis of Variance Table
##
## Model 1: raoQ ~ RDA1 + I(RDA1^2) + RDA2 + I(RDA2^2)
## Model 2: raoQ ~ RDA1 + RDA2
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 426 2278.4
## 2 428 2955.8 -2 -677.43 63.332 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1AIC(lm.rda.1, lm.rda.3)## df AIC
## lm.rda.1 6 1952.784
## lm.rda.3 4 2060.978summary(lm.rda.3) #r2 = 0.07, p<0.001##
## Call:
## lm(formula = raoQ ~ RDA1 + RDA2, data = arctic.cwm.rda.sites)
##
## Residuals:
## Min 1Q Median 3Q Max
## -12.6580 -1.5741 0.3088 1.6705 6.4160
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 19.9029 0.1266 157.232 < 2e-16 ***
## RDA1 -7.1859 1.6814 -4.274 2.37e-05 ***
## RDA2 -5.0475 1.3524 -3.732 0.000215 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.628 on 428 degrees of freedom
## Multiple R-squared: 0.07081, Adjusted R-squared: 0.06647
## F-statistic: 16.31 on 2 and 428 DF, p-value: 1.494e-07plot(lm.rda.1)
lm.rda.4<-lm(raoQ~ RDA1 + I(RDA1^2) + RDA2 + I(RDA2^2) + (RDA1*RDA2),
data=arctic.cwm.rda.sites)
anova(lm.rda.4, lm.rda.1) ##no difference## Analysis of Variance Table
##
## Model 1: raoQ ~ RDA1 + I(RDA1^2) + RDA2 + I(RDA2^2) + (RDA1 * RDA2)
## Model 2: raoQ ~ RDA1 + I(RDA1^2) + RDA2 + I(RDA2^2)
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 425 2278.3
## 2 426 2278.4 -1 -0.018545 0.0035 0.9531AIC(lm.rda.4, lm.rda.1) ##lm.rda.1 but no real difference## df AIC
## lm.rda.4 7 1954.781
## lm.rda.1 6 1952.784summary(lm.rda.4)##
## Call:
## lm(formula = raoQ ~ RDA1 + I(RDA1^2) + RDA2 + I(RDA2^2) + (RDA1 *
## RDA2), data = arctic.cwm.rda.sites)
##
## Residuals:
## Min 1Q Median 3Q Max
## -9.6823 -1.4974 0.2442 1.4698 8.2388
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 20.9797 0.1482 141.596 < 2e-16 ***
## RDA1 -6.8614 1.5424 -4.448 1.11e-05 ***
## I(RDA1^2) -112.8587 13.2707 -8.504 3.14e-16 ***
## RDA2 -7.9307 1.2403 -6.394 4.26e-10 ***
## I(RDA2^2) -49.8927 8.3717 -5.960 5.31e-09 ***
## RDA1:RDA2 0.9168 15.5873 0.059 0.953
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.315 on 425 degrees of freedom
## Multiple R-squared: 0.2838, Adjusted R-squared: 0.2753
## F-statistic: 33.68 on 5 and 425 DF, p-value: < 2.2e-16lm.rda.5<-lm(raoQ~ RDA1 + I(RDA1^2) + RDA2 + I(RDA2^2) + (I(RDA1^2)*I(RDA2^2)),
data=arctic.cwm.rda.sites)
anova(lm.rda.5, lm.rda.1) ##no difference## Analysis of Variance Table
##
## Model 1: raoQ ~ RDA1 + I(RDA1^2) + RDA2 + I(RDA2^2) + (I(RDA1^2) * I(RDA2^2))
## Model 2: raoQ ~ RDA1 + I(RDA1^2) + RDA2 + I(RDA2^2)
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 425 2278.3
## 2 426 2278.4 -1 -0.098315 0.0183 0.8923AIC(lm.rda.5, lm.rda.1) ##lm.rda.1 but no real difference## df AIC
## lm.rda.5 7 1954.766
## lm.rda.1 6 1952.784summary(lm.rda.5)##
## Call:
## lm(formula = raoQ ~ RDA1 + I(RDA1^2) + RDA2 + I(RDA2^2) + (I(RDA1^2) *
## I(RDA2^2)), data = arctic.cwm.rda.sites)
##
## Residuals:
## Min 1Q Median 3Q Max
## -9.6776 -1.5019 0.2409 1.4733 8.2384
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 20.9887 0.1645 127.582 < 2e-16 ***
## RDA1 -6.8667 1.5137 -4.536 7.46e-06 ***
## I(RDA1^2) -114.1726 16.1237 -7.081 5.96e-12 ***
## RDA2 -7.9131 1.2405 -6.379 4.66e-10 ***
## I(RDA2^2) -50.6387 10.3668 -4.885 1.47e-06 ***
## I(RDA1^2):I(RDA2^2) 79.7839 589.1317 0.135 0.892
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.315 on 425 degrees of freedom
## Multiple R-squared: 0.2838, Adjusted R-squared: 0.2754
## F-statistic: 33.68 on 5 and 425 DF, p-value: < 2.2e-16##kept with lm.rda.1
lm.step<-stepAIC(lm.rda.5)## Start: AIC=729.64
## raoQ ~ RDA1 + I(RDA1^2) + RDA2 + I(RDA2^2) + (I(RDA1^2) * I(RDA2^2))
##
## Df Sum of Sq RSS AIC
## - I(RDA1^2):I(RDA2^2) 1 0.098 2278.4 727.66
## <none> 2278.3 729.64
## - RDA1 1 110.313 2388.6 748.02
## - RDA2 1 218.141 2496.4 767.05
##
## Step: AIC=727.66
## raoQ ~ RDA1 + I(RDA1^2) + RDA2 + I(RDA2^2)
##
## Df Sum of Sq RSS AIC
## <none> 2278.4 727.66
## - RDA1 1 111.17 2389.5 746.19
## - I(RDA2^2) 1 212.95 2491.3 764.17
## - RDA2 1 227.59 2505.9 766.70
## - I(RDA1^2) 1 390.85 2669.2 793.90##Step: AIC=727.66; raoQ ~ RDA1 + I(RDA1^2) + RDA2 + I(RDA2^2)
lm.step##
## Call:
## lm(formula = raoQ ~ RDA1 + I(RDA1^2) + RDA2 + I(RDA2^2), data = arctic.cwm.rda.sites)
##
## Coefficients:
## (Intercept) RDA1 I(RDA1^2) RDA2 I(RDA2^2)
## 20.979 -6.880 -112.923 -7.944 -49.728raoQ.rda.1<-
ggplot()+
geom_smooth(data=arctic.cwm.rda.sites,
aes(x=RDA1, y=raoQ),
formula = y~poly(x,2),
method = "lm")+
geom_point(data=arctic.cwm.rda.sites,
aes(x=RDA1, y=raoQ, color=cluster))+
scale_color_manual(values=clus.cols.fuzzy)+
theme_bw()+
labs(x="RDA1", y="Rao's Quadratic Entropy")
raoQ.rda.1
raoQ.rda.2<-
ggplot()+
geom_smooth(data=arctic.cwm.rda.sites,
aes(x=RDA2, y=raoQ), formula = y~poly(x,2),
method = "lm")+
geom_point(data=arctic.cwm.rda.sites,
aes(x=RDA2, y=raoQ, color=cluster))+
scale_color_manual(values=clus.cols.fuzzy)+
theme_bw()+
labs(x="RDA2", y="Rao's Quadratic Entropy")
raoQ.rda.2
fig3.rao<-ggarrange(fuzzy.cluster.CI, raoQ.rda.1,raoQ.rda.2,
nrow=3, common.legend = TRUE, legend = "bottom",
align = "v", labels = c("b","c","d"))Diversity metrics compared to environmental gradients for taxonomic diversity
##taxonomy
arctic.cwm.cca.sites.tax$simpson<-as.numeric(arctic_env_tax$simpson)
lm.cca.0<-lm(simpson~ 1,
data=arctic.cwm.cca.sites.tax)
lm.cca.1<-lm(simpson~ CCA1 + I(CCA1^2) + CCA2 + I(CCA2^2),
data=arctic.cwm.cca.sites.tax)
anova(lm.cca.0, lm.cca.1)## Analysis of Variance Table
##
## Model 1: simpson ~ 1
## Model 2: simpson ~ CCA1 + I(CCA1^2) + CCA2 + I(CCA2^2)
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 430 1.9641
## 2 426 1.5530 4 0.41117 28.197 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(lm.cca.1)##
## Call:
## lm(formula = simpson ~ CCA1 + I(CCA1^2) + CCA2 + I(CCA2^2), data = arctic.cwm.cca.sites.tax)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.39112 -0.01927 0.01368 0.03797 0.11019
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.8759449 0.0040733 215.045 < 2e-16 ***
## CCA1 -0.0206280 0.0024860 -8.298 1.41e-15 ***
## I(CCA1^2) -0.0088179 0.0012681 -6.953 1.35e-11 ***
## CCA2 -0.0111958 0.0023405 -4.784 2.38e-06 ***
## I(CCA2^2) -0.0042692 0.0007753 -5.507 6.34e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.06038 on 426 degrees of freedom
## Multiple R-squared: 0.2093, Adjusted R-squared: 0.2019
## F-statistic: 28.2 on 4 and 426 DF, p-value: < 2.2e-16lm.cca.step<-stepAIC(lm.cca.1)## Start: AIC=-2414.78
## simpson ~ CCA1 + I(CCA1^2) + CCA2 + I(CCA2^2)
##
## Df Sum of Sq RSS AIC
## <none> 1.5530 -2414.8
## - CCA2 1 0.083416 1.6364 -2394.2
## - I(CCA2^2) 1 0.110539 1.6635 -2387.1
## - I(CCA1^2) 1 0.176258 1.7292 -2370.4
## - CCA1 1 0.251004 1.8040 -2352.2summary(lm.cca.step)##
## Call:
## lm(formula = simpson ~ CCA1 + I(CCA1^2) + CCA2 + I(CCA2^2), data = arctic.cwm.cca.sites.tax)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.39112 -0.01927 0.01368 0.03797 0.11019
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.8759449 0.0040733 215.045 < 2e-16 ***
## CCA1 -0.0206280 0.0024860 -8.298 1.41e-15 ***
## I(CCA1^2) -0.0088179 0.0012681 -6.953 1.35e-11 ***
## CCA2 -0.0111958 0.0023405 -4.784 2.38e-06 ***
## I(CCA2^2) -0.0042692 0.0007753 -5.507 6.34e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.06038 on 426 degrees of freedom
## Multiple R-squared: 0.2093, Adjusted R-squared: 0.2019
## F-statistic: 28.2 on 4 and 426 DF, p-value: < 2.2e-16lm.cca.5<-lm(simpson~ CCA1 + I(CCA1^2) + CCA2 + I(CCA2^2),
data=arctic.cwm.cca.sites.tax)
summary(lm.rda.5)##
## Call:
## lm(formula = raoQ ~ RDA1 + I(RDA1^2) + RDA2 + I(RDA2^2) + (I(RDA1^2) *
## I(RDA2^2)), data = arctic.cwm.rda.sites)
##
## Residuals:
## Min 1Q Median 3Q Max
## -9.6776 -1.5019 0.2409 1.4733 8.2384
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 20.9887 0.1645 127.582 < 2e-16 ***
## RDA1 -6.8667 1.5137 -4.536 7.46e-06 ***
## I(RDA1^2) -114.1726 16.1237 -7.081 5.96e-12 ***
## RDA2 -7.9131 1.2405 -6.379 4.66e-10 ***
## I(RDA2^2) -50.6387 10.3668 -4.885 1.47e-06 ***
## I(RDA1^2):I(RDA2^2) 79.7839 589.1317 0.135 0.892
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.315 on 425 degrees of freedom
## Multiple R-squared: 0.2838, Adjusted R-squared: 0.2754
## F-statistic: 33.68 on 5 and 425 DF, p-value: < 2.2e-16##kept with lm.rda.1
lm.cca.step<-stepAIC(lm.cca.5)## Start: AIC=-2414.78
## simpson ~ CCA1 + I(CCA1^2) + CCA2 + I(CCA2^2)
##
## Df Sum of Sq RSS AIC
## <none> 1.5530 -2414.8
## - CCA2 1 0.083416 1.6364 -2394.2
## - I(CCA2^2) 1 0.110539 1.6635 -2387.1
## - I(CCA1^2) 1 0.176258 1.7292 -2370.4
## - CCA1 1 0.251004 1.8040 -2352.2arctic.cwm.cca.sites.tax$CCA1<-as.numeric(arctic.cwm.cca.sites.tax$CCA1)
arctic.cwm.cca.sites.tax$CCA2<-as.numeric(arctic.cwm.cca.sites.tax$CCA2)
simpson.cca.1<-
ggplot()+
geom_smooth(data=arctic.cwm.cca.sites.tax,
aes(x=CCA1, y=simpson), formula = y~poly(x,2),
method = "lm" )+
geom_point(data=arctic.cwm.cca.sites.tax,
aes(x=CCA1, y=simpson, color=cluster))+
scale_color_manual(values=clus.cols.tax)+
theme_bw()+
labs(x="CCA1", y="Simpson Diversity Index")
simpson.cca.1
simpson.cca.2<-
ggplot()+
geom_smooth(data=arctic.cwm.cca.sites.tax,
aes(x=CCA2, y=simpson), formula = y~poly(x,2),
method = "lm")+
geom_point(data=arctic.cwm.cca.sites.tax,
aes(x=CCA2, y=simpson, color=cluster))+
scale_color_manual(values=clus.cols.tax)+
theme_bw()+
labs(x="CCA2", y="Simpson Diversity Index")
simpson.cca.2
fig5.simpson<-ggarrange(tax.cluster.CI, simpson.cca.1,simpson.cca.2,
nrow=3, common.legend = TRUE, legend = "bottom",
align = "v", labels = c("b","c","d"))
fig5.simpson
Supplementary Figures
arctic_env_fd_2<-as.data.frame(cbind(arctic_env_fd,arctic_env_tax$simpson))
View(arctic_env_fd_2)
names(arctic_env_fd_2)[10]<-"Simpson (1-D)"
arctic_env_fd_2$BottomTemp<-as.numeric(arctic_env_fd_2$BottomTemp)
arctic_env_fd_2$Year<-as.factor(arctic_env_fd_2$Year)
arctic_env_fd_melt<-melt(arctic_env_fd_2,
value.name = "value")## Using Year, cluster as id variablesView(arctic_env_fd_melt)
supp_fig_1<-ggplot(data=arctic_env_fd_melt,
aes(x=Year, y=value))+
geom_boxplot(aes(group=Year))+
facet_wrap(~variable, scales = "free", ncol = 2)+
theme_bw()
supp_fig_1