Data imports & clean-up

Assuming data files are in the working directory

getwd()
## [1] "/Users/lauren/Downloads"
arctic_traits_090424<-as.data.frame(read.csv("arctic_traits_090424_b.csv",
                               row.names = 1))
View(arctic_traits_090424)
arctic_cpue_090424<-as.data.frame(read.csv("biomass.transformed.matrix.csv",
                             row.names = 1))
arctic_env_090424<-as.data.frame(read.csv("arctic_env_090424_b.csv",
                            row.names = 1))

# Six species with no info on age-at-maturity:
arctic_traits_090424[,1:3][rowSums(arctic_traits_090424[,1:3])==0,]
##                          low.age.at.maturity medium.age.at.maturity
## Anisarchus.medius                          0                      0
## Enophrys.diceraus                          0                      0
## Eumesogrammus.praecisus                    0                      0
## Lycodes.spp.                               0                      0
## Stichaeus.punctatus                        0                      0
## Trichocottus.brashnikovi                   0                      0
##                          high.age.at.maturity
## Anisarchus.medius                           0
## Enophrys.diceraus                           0
## Eumesogrammus.praecisus                     0
## Lycodes.spp.                                0
## Stichaeus.punctatus                         0
## Trichocottus.brashnikovi                    0
##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)

##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))
dim(arctic_env_090424)
## [1] 431   7
##create weighted traits matrix
arctic_traits_090424_b<-arctic_traits_090424_b[,-21]
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)))

# Replace NAs (resulting from traits that have a 0 for all modalities)
arctic.traits.fuzzy[is.na(arctic.traits.fuzzy)] <- 0

Fuzzy correspondence analysis

arctic.fca <- dudi.fca(arctic.traits.fuzzy, scannf = F, nf=2)

scatter(arctic.fca, csub = 3, clab.moda = 1.5) 

summary(arctic.fca) 
## Class: fca dudi
## Call: dudi.fca(df = arctic.traits.fuzzy, scannf = F, nf = 2)
## 
## Total inertia: 0.1995
## 
## Eigenvalues:
##     Ax1     Ax2     Ax3     Ax4     Ax5 
## 0.07332 0.03414 0.01921 0.01633 0.01144 
## 
## Projected inertia (%):
##     Ax1     Ax2     Ax3     Ax4     Ax5 
##  36.758  17.116   9.633   8.187   5.737 
## 
## Cumulative projected inertia (%):
##     Ax1   Ax1:2   Ax1:3   Ax1:4   Ax1:5 
##   36.76   53.87   63.51   71.69   77.43 
## 
## (Only 5 dimensions (out of 17) are shown)
arctic.fca$cr
##            RS1         RS2
## FV1 0.05085099 0.008495159
## FV2 0.11368902 0.074642624
## FV3 0.10386336 0.005476967
## FV4 0.06769199 0.064917568
## FV5 0.05066136 0.024185292
## FV6 0.05316375 0.027127278

Assess 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")
plot(fuzzy.clus, labels=F)

fviz_nbclust(arctic.traits.fuzzy, kmeans,
             method = "sil", diss = fuzzy.dis, k.max = 6)

fviz_nbclust(arctic.traits.fuzzy, kmeans, 
             method = "wss", diss = fuzzy.dis, k.max = 6)

fviz_nbclust(arctic.traits.fuzzy, kmeans, 
             method = "gap_stat", diss = fuzzy.dis, k.max = 6)

# All methods point towards 3 clusters!
set.seed(1)
fuzzy.clus.kmeans<-kmeans(arctic.traits.fuzzy, centers = 3)

fviz_cluster(fuzzy.clus.kmeans, data = arctic.traits.fuzzy, labelsize=0) +
  scale_color_manual(values = c("#0D7EF9", "#98A100", "#FD1C26","orange4")) +
  scale_fill_manual(values = c("#0D7EF9", "#98A100", "#FD1C26","orange4"))

There is good separation among the three clusters in just the first two PCs of functional trait composition (CWM matrix), although a PCA based on Euclidean distances is not really an appropriate ordination method here.

Plot results from cluster analysis

## [1] TRUE

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.577640
RDA_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)
## $r.squared
## [1] 0.1674586
## 
## $adj.r.squared
## [1] 0.1556773
plot(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 ' ' 1
anova(RDA_model.arctic, 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 ~ 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.250    
## RDA4       1 0.000828  1.5758  0.601    
## RDA5       1 0.000269  0.5124  0.981    
## RDA6       1 0.000091  0.1737  0.992    
## Residual 424 0.222885                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(RDA_model.arctic, 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 ~ Depth + doy_retreat + BottomTemp + BottomSalinity + Latitude + Longitude, data = arctic_env)
##                 Df Variance       F Pr(>F)    
## Depth            1 0.001999  3.8037  0.004 ** 
## doy_retreat      1 0.001124  2.1382  0.053 .  
## BottomTemp       1 0.002485  4.7277  0.003 ** 
## BottomSalinity   1 0.000334  0.6349  0.657    
## 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 ' ' 1

The 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. 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.0200 *  
## doy_retreat      1 0.001124  2.1382 0.6354    
## BottomTemp       1 0.002485  4.7277 0.0208 *  
## BottomSalinity   1 0.000334  0.6349 0.7632    
## 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 ' ' 1

The 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.1534161
RDA_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.1058986

About 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.128405

About 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   41
RDA_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.08819173
anova(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 ' ' 1

There 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.729275
RDA_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.2324452
anova(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.0082 ** 
## BottomTemp       1 0.002261  4.7313 0.0030 ** 
## BottomSalinity   1 0.000198  0.4143 0.8590    
## Depth            1 0.001544  3.2317 0.0148 *  
## 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 ' ' 1
plot(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.004 **
## ---
## 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.2283473
step.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.007 ** 
## Residual     419 0.201092                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(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.096 .  
## RDA5       1 0.002300  4.7914  0.099 .  
## RDA6       1 0.001255  2.6152  0.564    
## RDA7       1 0.000968  2.0171  0.718    
## RDA8       1 0.000474  0.9885  0.978    
## RDA9       1 0.000168  0.3508  0.999    
## RDA10      1 0.000051  0.1062  1.000    
## RDA11      1 0.000010  0.0198           
## Residual 419 0.201092                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(step.best)

(R2.fun.best <- RsquareAdj(step.best)$r.sq)
## [1] 0.2488631

Best Model for functional traits

The resulting ‘best’ model eliminates day of ice retreat 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

## cumulative contributions of most influential species:
## 
## $Intermediate_Southern
##     fusiform...normal   very.high.longevity medium.high.fecundity 
##             0.0700453             0.1368361             0.1954519 
##              juvenile      medium.body.size   very.high.fecundity 
##             0.2529954             0.3096603             0.3643692 
##       small.body.size                  flat      medium.longevity 
##             0.4189860             0.4712934             0.5230296 
##  high.age.at.maturity      medium.fecundity             sub.adult 
##             0.5737627             0.6238831             0.6705561 
##              eel.like 
##             0.7155096 
## 
## $Intermediate_Northern
##   very.high.longevity     fusiform...normal              juvenile 
##            0.08500226            0.15082343            0.21580073 
##              eel.like medium.high.fecundity       small.body.size 
##            0.27892946            0.33938269            0.39969286 
##      medium.longevity      medium.fecundity      medium.body.size 
##            0.45861152            0.51526972            0.56902686 
##        high.longevity             sub.adult                adult1 
##            0.61623060            0.66079559            0.70041729 
## 
## $Southern_Northern
##  very.high.longevity             juvenile      small.body.size 
##            0.1052709            0.1825539            0.2533913 
##     medium.longevity     medium.body.size high.age.at.maturity 
##            0.3233804            0.3740326            0.4235775 
##            sub.adult       high.longevity  very.high.fecundity 
##            0.4703233            0.5156984            0.5590295 
##             eel.like     medium.fecundity                 flat 
##            0.6023426            0.6423036            0.6816705 
##      large.body.size 
##            0.7206724
## 
## Contrast: Intermediate_Southern 
## 
##                         average       sd    ratio      ava      avb cumsum    p
## fusiform...normal      0.014822 0.011310 1.310500 0.469200 0.315400  0.070 0.01
## very.high.longevity    0.014133 0.011613 1.217000 0.357800 0.507100  0.137 1.00
## medium.high.fecundity  0.012403 0.010149 1.222200 0.382500 0.275000  0.196 0.02
## juvenile               0.012177 0.008748 1.392000 0.617100 0.727000  0.253 1.00
## medium.body.size       0.011991 0.009730 1.232300 0.565100 0.640200  0.310 0.54
## very.high.fecundity    0.011577 0.008540 1.355600 0.046300 0.177700  0.364 0.01
## small.body.size        0.011557 0.008369 1.381000 0.261200 0.164500  0.419 1.00
## flat                   0.011069 0.007287 1.519000 0.037600 0.162600  0.471 0.01
## medium.longevity       0.010948 0.007887 1.388100 0.322400 0.222400  0.523 1.00
## high.age.at.maturity   0.010735 0.007192 1.492600 0.108000 0.222700  0.574 0.01
## medium.fecundity       0.010606 0.008038 1.319400 0.281200 0.336000  0.624 0.90
## sub.adult              0.009876 0.007167 1.378000 0.269000 0.176500  0.671 0.62
## eel.like               0.009512 0.007772 1.223900 0.376200 0.416200  0.716 1.00
## medium.age.at.maturity 0.009450 0.007265 1.300700 0.699200 0.630500  0.760 0.05
## large.body.size        0.008828 0.006742 1.309500 0.173600 0.195300  0.802 0.19
## high.longevity         0.007561 0.006524 1.159000 0.319800 0.270500  0.838 1.00
## low.medium.fecundity   0.007211 0.005204 1.385700 0.116300 0.054000  0.872 0.97
## low.age.at.maturity    0.006942 0.005767 1.203800 0.192800 0.146800  0.904 1.00
## high.fecundity         0.006550 0.005080 1.289300 0.164800 0.152300  0.935 0.23
## adult1                 0.006065 0.004457 1.360700 0.113900 0.096600  0.964 1.00
## deep.and.or.short      0.005005 0.004131 1.211400 0.105000 0.093700  0.988 0.98
## compressed             0.001531 0.001805 0.848000 0.012100 0.012100  0.995 0.01
## low.fecundity          0.001056 0.001947 0.542700 0.009000 0.005000  1.000 1.00
##                          
## fusiform...normal      **
## very.high.longevity      
## medium.high.fecundity  * 
## juvenile                 
## medium.body.size         
## very.high.fecundity    **
## small.body.size          
## flat                   **
## medium.longevity         
## high.age.at.maturity   **
## medium.fecundity         
## sub.adult                
## eel.like                 
## medium.age.at.maturity * 
## large.body.size          
## high.longevity           
## low.medium.fecundity     
## low.age.at.maturity      
## high.fecundity           
## adult1                   
## deep.and.or.short        
## compressed             **
## low.fecundity            
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Contrast: Intermediate_Northern 
## 
##                         average       sd    ratio      ava      avb cumsum    p
## very.high.longevity    0.018713 0.011361 1.647100 0.357800 0.156200  0.085 0.01
## fusiform...normal      0.014490 0.010393 1.394300 0.469200 0.316500  0.151 0.01
## juvenile               0.014305 0.010877 1.315100 0.617100 0.472800  0.216 0.04
## eel.like               0.013898 0.009773 1.422000 0.376200 0.527800  0.279 0.01
## medium.high.fecundity  0.013309 0.009902 1.344100 0.382500 0.260200  0.339 0.01
## small.body.size        0.013277 0.010715 1.239100 0.261200 0.396200  0.400 0.13
## medium.longevity       0.012971 0.010257 1.264600 0.322400 0.446400  0.459 0.05
## medium.fecundity       0.012473 0.009135 1.365500 0.281200 0.395800  0.515 0.01
## medium.body.size       0.011835 0.009921 1.192900 0.565100 0.505800  0.569 0.67
## high.longevity         0.010392 0.008591 1.209600 0.319800 0.397400  0.616 0.01
## sub.adult              0.009811 0.007666 1.279800 0.269000 0.318700  0.661 0.71
## adult1                 0.008723 0.006506 1.340800 0.113900 0.208500  0.700 0.01
## medium.age.at.maturity 0.008454 0.006792 1.244700 0.699200 0.706800  0.739 0.98
## low.age.at.maturity    0.008355 0.006731 1.241300 0.192800 0.221900  0.777 0.02
## large.body.size        0.008139 0.006352 1.281400 0.173600 0.098000  0.814 0.95
## low.medium.fecundity   0.007641 0.006106 1.251400 0.116300 0.166500  0.848 0.41
## high.age.at.maturity   0.007276 0.005070 1.435000 0.108000 0.071300  0.881 1.00
## high.fecundity         0.006654 0.005175 1.285900 0.164800 0.101600  0.912 0.17
## deep.and.or.short      0.005724 0.004947 1.157200 0.105000 0.107900  0.938 0.10
## very.high.fecundity    0.004821 0.004508 1.069300 0.046300 0.041400  0.960 1.00
## flat                   0.004719 0.004820 0.979100 0.037600 0.043700  0.981 1.00
## low.fecundity          0.002966 0.003363 0.881900 0.009000 0.034500  0.995 0.01
## compressed             0.001204 0.001879 0.641000 0.012100 0.004100  1.000 0.90
##                          
## very.high.longevity    **
## fusiform...normal      **
## juvenile               * 
## eel.like               **
## medium.high.fecundity  **
## small.body.size          
## medium.longevity       * 
## medium.fecundity       **
## medium.body.size         
## high.longevity         **
## sub.adult                
## adult1                 **
## medium.age.at.maturity   
## low.age.at.maturity    * 
## large.body.size          
## low.medium.fecundity     
## high.age.at.maturity     
## high.fecundity           
## deep.and.or.short      . 
## very.high.fecundity      
## flat                     
## low.fecundity          **
## compressed               
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Contrast: Southern_Northern 
## 
##                         average       sd    ratio      ava      avb cumsum    p
## very.high.longevity    0.029450 0.013988 2.105500 0.507100 0.156200  0.105 0.01
## juvenile               0.021621 0.012562 1.721100 0.727000 0.472800  0.183 0.01
## small.body.size        0.019817 0.012428 1.594600 0.164500 0.396200  0.253 0.01
## medium.longevity       0.019580 0.012083 1.620400 0.222400 0.446400  0.323 0.01
## medium.body.size       0.014170 0.011617 1.219800 0.640200 0.505800  0.374 0.01
## high.age.at.maturity   0.013861 0.008229 1.684300 0.222700 0.071300  0.424 0.01
## sub.adult              0.013078 0.008926 1.465100 0.176500 0.318700  0.470 0.01
## high.longevity         0.012694 0.010043 1.264000 0.270500 0.397400  0.516 0.01
## very.high.fecundity    0.012122 0.008612 1.407500 0.177700 0.041400  0.559 0.01
## eel.like               0.012117 0.008401 1.442400 0.416200 0.527800  0.602 0.01
## medium.fecundity       0.011179 0.008436 1.325200 0.336000 0.395800  0.642 0.23
## flat                   0.011013 0.007369 1.494600 0.162600 0.043700  0.682 0.01
## large.body.size        0.010911 0.007737 1.410300 0.195300 0.098000  0.721 0.01
## medium.age.at.maturity 0.010608 0.008449 1.255500 0.630500 0.706800  0.759 0.01
## low.medium.fecundity   0.010485 0.007006 1.496600 0.054000 0.166500  0.796 0.01
## medium.high.fecundity  0.010298 0.007879 1.307100 0.275000 0.260200  0.833 1.00
## adult1                 0.010125 0.007197 1.406900 0.096600 0.208500  0.869 0.01
## fusiform...normal      0.009932 0.007716 1.287100 0.315400 0.316500  0.905 1.00
## low.age.at.maturity    0.009488 0.007607 1.247300 0.146800 0.221900  0.939 0.01
## high.fecundity         0.007183 0.005197 1.382100 0.152300 0.101600  0.964 0.01
## deep.and.or.short      0.005842 0.004870 1.199600 0.093700 0.107900  0.985 0.02
## low.fecundity          0.002961 0.003555 0.832800 0.005000 0.034500  0.996 0.01
## compressed             0.001223 0.002088 0.586000 0.012100 0.004100  1.000 0.82
##                          
## very.high.longevity    **
## juvenile               **
## small.body.size        **
## medium.longevity       **
## medium.body.size       **
## high.age.at.maturity   **
## sub.adult              **
## high.longevity         **
## very.high.fecundity    **
## eel.like               **
## medium.fecundity         
## flat                   **
## large.body.size        **
## medium.age.at.maturity **
## low.medium.fecundity   **
## medium.high.fecundity    
## adult1                 **
## fusiform...normal        
## low.age.at.maturity    **
## high.fecundity         **
## deep.and.or.short      * 
## low.fecundity          **
## compressed               
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation: free
## Number of permutations: 99

Taxonomic clustering

  1. 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)

There are two obvious clusters, but it is unclear if the 2nd cluster should be grouped into more smaller clusters and, if so, into how many. There appears to be a fairly obvious cutoff at 5 clusters (at a height of about 3.5)

Determine number of clusters

fviz_nbclust(arctic_cpue, kmeans, 
             method = "silh", diss = tax.dis, k.max = 8) # 2 (or 4)

fviz_nbclust(arctic_cpue, kmeans, 
             method = "wss", diss = tax.dis, k.max = 8) # 2 or 5?

fviz_nbclust(arctic_cpue, kmeans, 
             method = "gap", diss = tax.dis, k.max = 8) # 5

There is some ambiguity but results suggest either 2 or 5 clusters, consistent with the dendrogram. Nevertheless, for consistency with the functional clusters and for ease of comparison, we use three clusters:

set.seed(1)
tax.clus.kmeans.1 <- kmeans(arctic_cpue, centers = 4, 
                            iter.max = 100)

# Or, if using 5 clusters (messier!)
# tax.clus.kmeans<-kmeans(arctic_cpue, centers = 5)

# Visualize clusters using a PCA of the taxonomic data (which is unlikely to represent data very well)
fviz_cluster(tax.clus.kmeans.1, data = arctic_cpue, labelsize=0) +
  scale_color_manual(values = c("#FD1C26", "#98A100", "#0D7EF9","orange")) +
  scale_fill_manual(values = c("#FD1C26", "#98A100", "#0D7EF9","orange"))

Compared to the functional clusters, the PCA of the cluster results accounts for a limited proportion of the variability in species abundances and the clusters overlap substantially in the space of the first two PCs. This could be a result of the PCA not being an appropriate ordination method for the species abundance data (many zeros!).

Map clusters

## [1] TRUE

Stacked barplot of species compositions:

Ordination based on taxa

There was a comment from someone (Caitlin?): “I added distance to shore in addition to latitude and longitude as our spatial component we are specifically thinking about the Alaska Coastal Current compared to Anadyr so distance to shore in tandem with temp and sal should capture that?”

Franz: This makes sense, but Distance from shore was not in the data frame!

Also, I made a few changes below. I used the reduced data set (23 taxa) and used cca instead of rda, as cca allows for a Gaussian (dome-shaped) response of taxa to environmental gradients.

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.807384
CCA_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)
## [1] 0.1869642
plot(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 ' ' 1
anova(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.011 *  
## CCA8       1   0.01472  2.8747  0.306    
## CCA9       1   0.01289  2.5166  0.674    
## 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 ' ' 1
anova(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.048 *  
## 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 ' ' 1

The 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.

All individual terms except salinity are significant.

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.16058896
## + factor(Year)    0.07685339
## + BottomTemp      0.04456610
## + Latitude        0.03966945
## + BottomSalinity  0.02988888
## + Depth           0.02575235
## + Longitude       0.01878999
## + doy_retreat     0.01744261
## <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.07679814 
## Call: arctic_cpue ~ factor(Year) 
##  
##                  R2.adjusted
## <All variables>   0.16058896
## + BottomTemp      0.11985292
## + Latitude        0.11144451
## + BottomSalinity  0.10666028
## + Depth           0.10356078
## + doy_retreat     0.09172164
## + Longitude       0.09130296
## <none>            0.07679814
## 
##              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.1198166 
## Call: arctic_cpue ~ factor(Year) + BottomTemp 
##  
##                  R2.adjusted
## <All variables>    0.1605890
## + Latitude         0.1334319
## + Longitude        0.1319710
## + Depth            0.1315062
## + BottomSalinity   0.1260115
## + doy_retreat      0.1246509
## <none>             0.1198166
## 
##            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.133516 
## Call: arctic_cpue ~ factor(Year) + BottomTemp + Latitude 
##  
##                  R2.adjusted
## <All variables>    0.1605890
## + Depth            0.1468811
## + Longitude        0.1445690
## + BottomSalinity   0.1398196
## + doy_retreat      0.1374759
## <none>             0.1335160
## 
##         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.1469695 
## Call: arctic_cpue ~ factor(Year) + BottomTemp + Latitude + Depth 
##  
##                  R2.adjusted
## <All variables>    0.1605890
## + Longitude        0.1563189
## + BottomSalinity   0.1503971
## + doy_retreat      0.1499596
## <none>             0.1469695
## 
##             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.1563866 
## Call: arctic_cpue ~ factor(Year) + BottomTemp + Latitude + Depth +      Longitude 
##  
##                  R2.adjusted
## <All variables>    0.1605890
## + doy_retreat      0.1592537
## + BottomSalinity   0.1579646
## <none>             0.1563866
## 
##               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.1591602 
## Call: arctic_cpue ~ factor(Year) + BottomTemp + Latitude + Depth +      Longitude + doy_retreat 
##  
##                  R2.adjusted
## <All variables>    0.1605890
## + BottomSalinity   0.1605691
## <none>             0.1591602
## 
##                  Df    AIC      F Pr(>F)  
## + BottomSalinity  1 1960.3 1.7513  0.036 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Step: R2.adj= 0.1605687 
## Call: arctic_cpue ~ factor(Year) + BottomTemp + Latitude + Depth +      Longitude + doy_retreat + BottomSalinity 
## 
step.res.tax
## Call: cca(formula = arctic_cpue ~ factor(Year) + BottomTemp + Latitude +
## Depth + Longitude + doy_retreat + BottomSalinity, 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)
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 + BottomSalinity, data = arctic_env)
##                 Df ChiSquare      F Pr(>F)    
## factor(Year)     7   0.20697 5.7724  0.001 ***
## BottomTemp       1   0.02124 4.1463  0.001 ***
## Latitude         1   0.03760 7.3404  0.001 ***
## Depth            1   0.03094 6.0400  0.001 ***
## Longitude        1   0.02461 4.8036  0.001 ***
## doy_retreat      1   0.01194 2.3310  0.003 ** 
## BottomSalinity   1   0.00897 1.7513  0.044 *  
## Residual       417   2.13597                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(R2.best <- RsquareAdj(step.res.tax)$r.sq)
## [1] 0.1869642
plot(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.7%)

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)
## [1] 0.07813255

About 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)
## [1] 0.05700268

Only 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:

step.res.tax2 <- cca(arctic_cpue ~ factor(Year) + Condition(BottomTemp + 
                           Latitude + Depth + doy_retreat + Longitude), 
                         data = arctic_env)
step.res.tax2
## Call: cca(formula = arctic_cpue ~ factor(Year) + Condition(BottomTemp +
## Latitude + Depth + doy_retreat + Longitude), data = arctic_env)
## 
## -- Model Summary --
## 
##               Inertia Proportion Rank
## Total         2.62715    1.00000     
## Conditional   0.27208    0.10356    5
## Constrained   0.21014    0.07999    7
## Unconstrained 2.14494    0.81645   32
## 
## Inertia is scaled Chi-square
## 
## -- Eigenvalues --
## 
## Eigenvalues for constrained axes:
##    CCA1    CCA2    CCA3    CCA4    CCA5    CCA6    CCA7 
## 0.07257 0.04729 0.03567 0.02491 0.01740 0.00642 0.00588 
## 
## 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.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), data = arctic_env)
##           Df ChiSquare      F Pr(>F)    
## Model      7   0.21014 5.8502  0.001 ***
## Residual 418   2.14494                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(R2.7 <- RsquareAdj(step.res.tax2)$r.sq)
## [1] 0.07998691

There are additional significant differences among years that are not accounted for by other variables in the best model that account for about 8% of the variability in taxonomic composition

Plot results (ggplot)

Plot CCA ordination with species and environmental variables overlaid. (this seems to use different scaling for environmental vectors from default CCA plot?)

## [1] "BottomTemp"  "Latitude"    "Depth"       "Longitude"   "doy_retreat"
## CCA1 CCA2 
## 6.54 3.05

## 
## Contrast: Nearshore_Intermediate 
## 
##                               average      sd   ratio     ava     avb cumsum
## Myoxocephalus.scorpius        0.05948 0.03392 1.75320 6.48000 1.85100  0.090
## Gymnocanthus.tricuspis        0.05697 0.03755 1.51700 6.54300 2.38200  0.176
## Artediellus.scaber            0.04323 0.03514 1.23010 3.74000 0.70700  0.242
## Boreogadus.saida              0.04253 0.03551 1.19740 4.49300 3.32000  0.306
## Lumpenus.fabricii.sp.         0.03781 0.03302 1.14520 6.01600 4.40300  0.363
## Stichaeus.punctatus           0.03190 0.03290 0.96980 2.63800 0.76000  0.411
## Liparis.spp.                  0.03190 0.03394 0.93970 1.44600 2.19600  0.460
## Aspidophoroides.olrikii       0.03013 0.02921 1.03130 2.26200 1.62200  0.505
## Hippoglossoides.robustus      0.02870 0.02924 0.98160 1.38900 2.05600  0.549
## Eleginus.gracilis             0.02798 0.03339 0.83810 2.25100 0.63200  0.591
## Lycodes.spp.                  0.02333 0.02689 0.86750 1.26700 1.42300  0.626
## Nautichthys.pribilovius       0.02331 0.02660 0.87630 1.85100 0.44200  0.662
## Triglops.pingelii             0.02133 0.02691 0.79250 1.71900 0.38100  0.694
## Liparis.tunicatus             0.02082 0.02800 0.74330 1.11600 1.04600  0.725
## Podothecus.veternus           0.01867 0.02769 0.67430 1.49200 0.14400  0.754
## Anisarchus.medius             0.01649 0.02513 0.65610 0.59800 0.99700  0.779
## Limanda.aspera                0.01542 0.02659 0.58000 1.18300 0.28700  0.802
## Gymnelus.hemifasciatus        0.01408 0.02520 0.55880 0.91900 0.47300  0.823
## Liparis.gibbus                0.01369 0.02341 0.58460 0.86100 0.49100  0.844
## Trichocottus.brashnikovi      0.01234 0.02188 0.56400 0.96500 0.20800  0.863
## Ammodytes.hexapterus          0.01185 0.02463 0.48130 0.87000 0.19100  0.881
## Gymnelus.viridis              0.01045 0.02065 0.50590 0.66700 0.32300  0.896
## Icelus.spatula                0.01026 0.02202 0.46600 0.31300 0.63700  0.912
## Gadus.chalcogrammus           0.01025 0.02133 0.48050 0.33300 0.59100  0.927
## Pleuronectidae                0.00769 0.01992 0.38580 0.35800 0.34600  0.939
## Gadus.macrocephalus           0.00691 0.01835 0.37680 0.29100 0.34400  0.950
## Eumesogrammus.praecisus       0.00671 0.01783 0.37630 0.46600 0.15600  0.960
## Podothecus.accipenserinus     0.00550 0.01604 0.34260 0.37600 0.07700  0.968
## Aspidophoroides.monopterygius 0.00538 0.01406 0.38280 0.46700 0.07400  0.976
## Hemilepidotus.papilio         0.00488 0.01577 0.30960 0.39400 0.04200  0.984
## Limanda.proboscidea           0.00462 0.01442 0.32020 0.29500 0.11600  0.991
## Gymnelus.sp.                  0.00390 0.01322 0.29530 0.17600 0.17900  0.996
## Leptoclinus.maculatus         0.00239 0.01046 0.22850 0.14800 0.08900  1.000
##                                  p   
## Myoxocephalus.scorpius        0.01 **
## Gymnocanthus.tricuspis        0.01 **
## Artediellus.scaber            0.01 **
## Boreogadus.saida              0.01 **
## Lumpenus.fabricii.sp.         1.00   
## Stichaeus.punctatus           0.01 **
## Liparis.spp.                  0.02 * 
## Aspidophoroides.olrikii       1.00   
## Hippoglossoides.robustus      1.00   
## Eleginus.gracilis             0.01 **
## Lycodes.spp.                  1.00   
## Nautichthys.pribilovius       0.01 **
## Triglops.pingelii             0.01 **
## Liparis.tunicatus             0.34   
## Podothecus.veternus           0.01 **
## Anisarchus.medius             1.00   
## Limanda.aspera                0.01 **
## Gymnelus.hemifasciatus        1.00   
## Liparis.gibbus                0.38   
## Trichocottus.brashnikovi      0.01 **
## Ammodytes.hexapterus          0.05 * 
## Gymnelus.viridis              0.41   
## Icelus.spatula                0.98   
## Gadus.chalcogrammus           0.33   
## Pleuronectidae                0.01 **
## Gadus.macrocephalus           0.01 **
## Eumesogrammus.praecisus       0.02 * 
## Podothecus.accipenserinus     0.01 **
## Aspidophoroides.monopterygius 0.92   
## Hemilepidotus.papilio         0.05 * 
## Limanda.proboscidea           0.09 . 
## Gymnelus.sp.                  0.06 . 
## Leptoclinus.maculatus         1.00   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Contrast: Nearshore_Offshore 
## 
##                               average      sd   ratio     ava     avb cumsum
## Anisarchus.medius             0.04419 0.02419 1.82690 0.59800 4.97000  0.083
## Hippoglossoides.robustus      0.04029 0.02479 1.62560 1.38900 5.15600  0.158
## Lycodes.spp.                  0.03334 0.02236 1.49150 1.26700 4.15900  0.220
## Artediellus.scaber            0.03215 0.02636 1.21990 3.74000 1.78900  0.280
## Myoxocephalus.scorpius        0.03023 0.02487 1.21560 6.48000 4.13800  0.337
## Boreogadus.saida              0.02972 0.02491 1.19310 4.49300 5.14800  0.393
## Aspidophoroides.olrikii       0.02755 0.02180 1.26420 2.26200 3.75300  0.444
## Stichaeus.punctatus           0.02449 0.02523 0.97080 2.63800 0.43000  0.490
## Lumpenus.fabricii.sp.         0.02264 0.02098 1.07920 6.01600 6.51300  0.532
## Eleginus.gracilis             0.02263 0.02524 0.89670 2.25100 0.96300  0.575
## Liparis.spp.                  0.02127 0.02486 0.85570 1.44600 1.56800  0.614
## Gymnocanthus.tricuspis        0.01951 0.01582 1.23380 6.54300 5.74800  0.651
## Liparis.tunicatus             0.01829 0.02227 0.82130 1.11600 1.43900  0.685
## Nautichthys.pribilovius       0.01773 0.01994 0.88930 1.85100 0.43300  0.718
## Triglops.pingelii             0.01760 0.02027 0.86810 1.71900 0.78600  0.751
## Podothecus.veternus           0.01711 0.02144 0.79810 1.49200 0.84800  0.783
## Gymnelus.hemifasciatus        0.01434 0.02048 0.70010 0.91900 1.01700  0.810
## Liparis.gibbus                0.01185 0.01792 0.66110 0.86100 0.72700  0.832
## Limanda.aspera                0.01174 0.02032 0.57790 1.18300 0.23200  0.854
## Aspidophoroides.monopterygius 0.01029 0.01658 0.62030 0.46700 0.90600  0.873
## Ammodytes.hexapterus          0.00959 0.01916 0.50070 0.87000 0.24700  0.891
## Trichocottus.brashnikovi      0.00933 0.01689 0.55230 0.96500 0.14000  0.909
## Gymnelus.viridis              0.00656 0.01482 0.44240 0.66700 0.10900  0.921
## Gadus.chalcogrammus           0.00566 0.01326 0.42680 0.33300 0.36100  0.932
## Leptoclinus.maculatus         0.00528 0.01314 0.40170 0.14800 0.48200  0.942
## Hemilepidotus.papilio         0.00523 0.01367 0.38270 0.39400 0.21800  0.951
## Icelus.spatula                0.00482 0.01251 0.38510 0.31300 0.27200  0.960
## Eumesogrammus.praecisus       0.00480 0.01356 0.35390 0.46600 0.09900  0.969
## Gadus.macrocephalus           0.00387 0.01272 0.30450 0.29100 0.14900  0.977
## Podothecus.accipenserinus     0.00363 0.01159 0.31350 0.37600 0.00000  0.983
## Limanda.proboscidea           0.00355 0.01064 0.33340 0.29500 0.10700  0.990
## Pleuronectidae                0.00354 0.01237 0.28630 0.35800 0.03700  0.997
## Gymnelus.sp.                  0.00186 0.00786 0.23710 0.17600 0.03000  1.000
##                                  p   
## Anisarchus.medius             0.02 * 
## Hippoglossoides.robustus      0.01 **
## Lycodes.spp.                  1.00   
## Artediellus.scaber            0.80   
## Myoxocephalus.scorpius        1.00   
## Boreogadus.saida              1.00   
## Aspidophoroides.olrikii       1.00   
## Stichaeus.punctatus           0.01 **
## Lumpenus.fabricii.sp.         1.00   
## Eleginus.gracilis             0.01 **
## Liparis.spp.                  1.00   
## Gymnocanthus.tricuspis        1.00   
## Liparis.tunicatus             0.91   
## Nautichthys.pribilovius       0.01 **
## Triglops.pingelii             0.01 **
## Podothecus.veternus           0.01 **
## Gymnelus.hemifasciatus        1.00   
## Liparis.gibbus                0.86   
## Limanda.aspera                0.01 **
## Aspidophoroides.monopterygius 0.01 **
## Ammodytes.hexapterus          0.46   
## Trichocottus.brashnikovi      0.01 **
## Gymnelus.viridis              1.00   
## Gadus.chalcogrammus           1.00   
## Leptoclinus.maculatus         0.84   
## Hemilepidotus.papilio         0.01 **
## Icelus.spatula                1.00   
## Eumesogrammus.praecisus       0.32   
## Gadus.macrocephalus           0.50   
## Podothecus.accipenserinus     0.03 * 
## Limanda.proboscidea           0.31   
## Pleuronectidae                0.53   
## Gymnelus.sp.                  0.77   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Contrast: Nearshore_Northern 
## 
##                               average      sd   ratio     ava     avb cumsum
## Gymnocanthus.tricuspis        0.07042 0.03349 2.10290 6.54300 0.79100  0.098
## Myoxocephalus.scorpius        0.06259 0.03246 1.92820 6.48000 1.34200  0.185
## Lumpenus.fabricii.sp.         0.05546 0.03568 1.55410 6.01600 1.80100  0.262
## Anisarchus.medius             0.05354 0.03062 1.74890 0.59800 4.90000  0.336
## Lycodes.spp.                  0.04496 0.03004 1.49680 1.26700 4.49100  0.399
## Artediellus.scaber            0.03940 0.03278 1.20180 3.74000 2.01700  0.454
## Boreogadus.saida              0.03396 0.02888 1.17610 4.49300 5.39500  0.501
## Aspidophoroides.olrikii       0.03028 0.02816 1.07530 2.26200 2.17600  0.543
## Stichaeus.punctatus           0.03007 0.03265 0.92090 2.63800 0.03200  0.585
## Eleginus.gracilis             0.02535 0.03288 0.77110 2.25100 0.13400  0.620
## Hippoglossoides.robustus      0.02420 0.02729 0.88680 1.38900 1.53400  0.654
## Liparis.spp.                  0.02410 0.03126 0.77100 1.44600 1.16000  0.687
## Nautichthys.pribilovius       0.02170 0.02544 0.85280 1.85100 0.22400  0.718
## Gymnelus.hemifasciatus        0.02124 0.02717 0.78200 0.91900 1.49500  0.747
## Triglops.pingelii             0.02024 0.02567 0.78840 1.71900 0.33200  0.775
## Podothecus.veternus           0.01738 0.02625 0.66210 1.49200 0.05700  0.799
## Liparis.tunicatus             0.01715 0.02636 0.65070 1.11600 0.57900  0.823
## Icelus.spatula                0.01527 0.02355 0.64830 0.31300 1.24100  0.845
## Ammodytes.hexapterus          0.01386 0.02546 0.54420 0.87000 0.48000  0.864
## Limanda.aspera                0.01352 0.02486 0.54380 1.18300 0.10100  0.882
## Liparis.gibbus                0.01351 0.02224 0.60770 0.86100 0.56600  0.901
## Gymnelus.viridis              0.01275 0.02156 0.59140 0.66700 0.69800  0.919
## Trichocottus.brashnikovi      0.01076 0.02060 0.52240 0.96500 0.02900  0.934
## Gadus.chalcogrammus           0.00733 0.01717 0.42700 0.33300 0.37200  0.944
## Eumesogrammus.praecisus       0.00608 0.01668 0.36470 0.46600 0.14500  0.953
## Leptoclinus.maculatus         0.00552 0.01466 0.37620 0.14800 0.41500  0.960
## Aspidophoroides.monopterygius 0.00549 0.01454 0.37800 0.46700 0.08400  0.968
## Podothecus.accipenserinus     0.00473 0.01475 0.32100 0.37600 0.02500  0.975
## Hemilepidotus.papilio         0.00428 0.01458 0.29350 0.39400 0.00000  0.981
## Limanda.proboscidea           0.00411 0.01358 0.30270 0.29500 0.07300  0.986
## Pleuronectidae                0.00395 0.01466 0.26970 0.35800 0.00000  0.992
## Gadus.macrocephalus           0.00320 0.01352 0.23690 0.29100 0.00000  0.996
## Gymnelus.sp.                  0.00274 0.01112 0.24630 0.17600 0.07200  1.000
##                                  p   
## Gymnocanthus.tricuspis        0.01 **
## Myoxocephalus.scorpius        0.01 **
## Lumpenus.fabricii.sp.         0.01 **
## Anisarchus.medius             0.01 **
## Lycodes.spp.                  0.01 **
## Artediellus.scaber            0.01 **
## Boreogadus.saida              0.97   
## Aspidophoroides.olrikii       1.00   
## Stichaeus.punctatus           0.01 **
## Eleginus.gracilis             0.01 **
## Hippoglossoides.robustus      1.00   
## Liparis.spp.                  1.00   
## Nautichthys.pribilovius       0.01 **
## Gymnelus.hemifasciatus        0.32   
## Triglops.pingelii             0.01 **
## Podothecus.veternus           0.01 **
## Liparis.tunicatus             0.96   
## Icelus.spatula                0.33   
## Ammodytes.hexapterus          0.01 **
## Limanda.aspera                0.01 **
## Liparis.gibbus                0.49   
## Gymnelus.viridis              0.03 * 
## Trichocottus.brashnikovi      0.01 **
## Gadus.chalcogrammus           0.95   
## Eumesogrammus.praecisus       0.05 * 
## Leptoclinus.maculatus         0.84   
## Aspidophoroides.monopterygius 0.96   
## Podothecus.accipenserinus     0.01 **
## Hemilepidotus.papilio         0.04 * 
## Limanda.proboscidea           0.09 . 
## Pleuronectidae                0.53   
## Gadus.macrocephalus           0.77   
## Gymnelus.sp.                  0.26   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Contrast: Intermediate_Offshore 
## 
##                               average      sd   ratio     ava     avb cumsum
## Anisarchus.medius             0.05374 0.03214 1.67210 0.99700 4.97000  0.091
## Hippoglossoides.robustus      0.04645 0.03269 1.42100 2.05600 5.15600  0.169
## Gymnocanthus.tricuspis        0.04640 0.03233 1.43510 2.38200 5.74800  0.248
## Boreogadus.saida              0.04261 0.03454 1.23360 3.32000 5.14800  0.320
## Myoxocephalus.scorpius        0.04212 0.03185 1.32250 1.85100 4.13800  0.391
## Lycodes.spp.                  0.04195 0.02909 1.44210 1.42300 4.15900  0.462
## Aspidophoroides.olrikii       0.03769 0.02806 1.34330 1.62200 3.75300  0.525
## Lumpenus.fabricii.sp.         0.03505 0.03156 1.11050 4.40300 6.51300  0.585
## Liparis.spp.                  0.03119 0.03154 0.98880 2.19600 1.56800  0.637
## Artediellus.scaber            0.02440 0.03010 0.81060 0.70700 1.78900  0.678
## Liparis.tunicatus             0.02286 0.02766 0.82630 1.04600 1.43900  0.717
## Gymnelus.hemifasciatus        0.01594 0.02665 0.59830 0.47300 1.01700  0.744
## Eleginus.gracilis             0.01594 0.02452 0.65020 0.63200 0.96300  0.771
## Liparis.gibbus                0.01254 0.02158 0.58120 0.49100 0.72700  0.792
## Stichaeus.punctatus           0.01243 0.02253 0.55150 0.76000 0.43000  0.813
## Triglops.pingelii             0.01179 0.02044 0.57700 0.38100 0.78600  0.833
## Podothecus.veternus           0.01110 0.02194 0.50590 0.14400 0.84800  0.852
## Aspidophoroides.monopterygius 0.01057 0.01993 0.53040 0.07400 0.90600  0.870
## Gadus.chalcogrammus           0.01055 0.02136 0.49390 0.59100 0.36100  0.887
## Icelus.spatula                0.00995 0.02127 0.46770 0.63700 0.27200  0.904
## Nautichthys.pribilovius       0.00902 0.01855 0.48640 0.44200 0.43300  0.919
## Leptoclinus.maculatus         0.00651 0.01695 0.38400 0.08900 0.48200  0.930
## Limanda.aspera                0.00581 0.01696 0.34270 0.28700 0.23200  0.940
## Gadus.macrocephalus           0.00559 0.01587 0.35210 0.34400 0.14900  0.950
## Gymnelus.viridis              0.00502 0.01491 0.33630 0.32300 0.10900  0.958
## Ammodytes.hexapterus          0.00471 0.01464 0.32180 0.19100 0.24700  0.966
## Pleuronectidae                0.00452 0.01546 0.29210 0.34600 0.03700  0.974
## Trichocottus.brashnikovi      0.00356 0.01180 0.30190 0.20800 0.14000  0.980
## Hemilepidotus.papilio         0.00321 0.01249 0.25740 0.04200 0.21800  0.985
## Eumesogrammus.praecisus       0.00306 0.01156 0.26460 0.15600 0.09900  0.990
## Limanda.proboscidea           0.00253 0.01162 0.21800 0.11600 0.10700  0.995
## Gymnelus.sp.                  0.00226 0.01023 0.22110 0.17900 0.03000  0.999
## Podothecus.accipenserinus     0.00092 0.00631 0.14520 0.07700 0.00000  1.000
##                                  p   
## Anisarchus.medius             0.01 **
## Hippoglossoides.robustus      0.01 **
## Gymnocanthus.tricuspis        0.01 **
## Boreogadus.saida              0.01 **
## Myoxocephalus.scorpius        0.22   
## Lycodes.spp.                  0.01 **
## Aspidophoroides.olrikii       0.01 **
## Lumpenus.fabricii.sp.         1.00   
## Liparis.spp.                  0.10 . 
## Artediellus.scaber            1.00   
## Liparis.tunicatus             0.03 * 
## Gymnelus.hemifasciatus        1.00   
## Eleginus.gracilis             0.68   
## Liparis.gibbus                0.71   
## Stichaeus.punctatus           0.99   
## Triglops.pingelii             0.95   
## Podothecus.veternus           0.48   
## Aspidophoroides.monopterygius 0.01 **
## Gadus.chalcogrammus           0.26   
## Icelus.spatula                1.00   
## Nautichthys.pribilovius       0.99   
## Leptoclinus.maculatus         0.57   
## Limanda.aspera                0.98   
## Gadus.macrocephalus           0.07 . 
## Gymnelus.viridis              1.00   
## Ammodytes.hexapterus          1.00   
## Pleuronectidae                0.22   
## Trichocottus.brashnikovi      0.98   
## Hemilepidotus.papilio         0.33   
## Eumesogrammus.praecisus       0.92   
## Limanda.proboscidea           0.73   
## Gymnelus.sp.                  0.63   
## Podothecus.accipenserinus     0.97   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Contrast: Intermediate_Northern 
## 
##                               average      sd   ratio     ava     avb cumsum
## Anisarchus.medius             0.07063 0.04353 1.62270 0.99700 4.90000  0.101
## Lycodes.spp.                  0.06166 0.04320 1.42750 1.42300 4.49100  0.190
## Lumpenus.fabricii.sp.         0.06053 0.04706 1.28630 4.40300 1.80100  0.276
## Boreogadus.saida              0.05560 0.04715 1.17910 3.32000 5.39500  0.356
## Hippoglossoides.robustus      0.04021 0.04164 0.96550 2.05600 1.53400  0.414
## Liparis.spp.                  0.04015 0.04346 0.92390 2.19600 1.16000  0.472
## Aspidophoroides.olrikii       0.03985 0.03936 1.01230 1.62200 2.17600  0.529
## Gymnocanthus.tricuspis        0.03976 0.03975 1.00030 2.38200 0.79100  0.586
## Myoxocephalus.scorpius        0.03761 0.04085 0.92080 1.85100 1.34200  0.640
## Artediellus.scaber            0.03544 0.04144 0.85530 0.70700 2.01700  0.691
## Gymnelus.hemifasciatus        0.02713 0.03913 0.69330 0.47300 1.49500  0.729
## Icelus.spatula                0.02444 0.03537 0.69080 0.63700 1.24100  0.764
## Liparis.tunicatus             0.02195 0.03341 0.65700 1.04600 0.57900  0.796
## Liparis.gibbus                0.01482 0.02852 0.51980 0.49100 0.56600  0.817
## Gadus.chalcogrammus           0.01481 0.03001 0.49360 0.59100 0.37200  0.838
## Gymnelus.viridis              0.01403 0.02751 0.51000 0.32300 0.69800  0.859
## Stichaeus.punctatus           0.01241 0.02784 0.44590 0.76000 0.03200  0.876
## Eleginus.gracilis             0.01136 0.02636 0.43100 0.63200 0.13400  0.893
## Triglops.pingelii             0.01008 0.02321 0.43450 0.38100 0.33200  0.907
## Ammodytes.hexapterus          0.00977 0.02581 0.37870 0.19100 0.48000  0.921
## Nautichthys.pribilovius       0.00942 0.02310 0.40770 0.44200 0.22400  0.935
## Leptoclinus.maculatus         0.00718 0.02023 0.35470 0.08900 0.41500  0.945
## Limanda.aspera                0.00568 0.01987 0.28610 0.28700 0.10100  0.953
## Pleuronectidae                0.00544 0.01986 0.27370 0.34600 0.00000  0.961
## Gadus.macrocephalus           0.00527 0.01794 0.29380 0.34400 0.00000  0.968
## Eumesogrammus.praecisus       0.00438 0.01568 0.27930 0.15600 0.14500  0.975
## Gymnelus.sp.                  0.00352 0.01522 0.23160 0.17900 0.07200  0.980
## Podothecus.veternus           0.00329 0.01580 0.20840 0.14400 0.05700  0.985
## Trichocottus.brashnikovi      0.00314 0.01275 0.24660 0.20800 0.02900  0.989
## Limanda.proboscidea           0.00288 0.01593 0.18060 0.11600 0.07300  0.993
## Aspidophoroides.monopterygius 0.00266 0.01382 0.19270 0.07400 0.08400  0.997
## Podothecus.accipenserinus     0.00150 0.00904 0.16540 0.07700 0.02500  0.999
## Hemilepidotus.papilio         0.00065 0.00644 0.10170 0.04200 0.00000  1.000
##                                  p   
## Anisarchus.medius             0.01 **
## Lycodes.spp.                  0.01 **
## Lumpenus.fabricii.sp.         0.01 **
## Boreogadus.saida              0.01 **
## Hippoglossoides.robustus      0.01 **
## Liparis.spp.                  0.01 **
## Aspidophoroides.olrikii       0.01 **
## Gymnocanthus.tricuspis        1.00   
## Myoxocephalus.scorpius        1.00   
## Artediellus.scaber            0.01 **
## Gymnelus.hemifasciatus        0.01 **
## Icelus.spatula                0.01 **
## Liparis.tunicatus             0.10 . 
## Liparis.gibbus                0.13   
## Gadus.chalcogrammus           0.01 **
## Gymnelus.viridis              0.01 **
## Stichaeus.punctatus           0.96   
## Eleginus.gracilis             1.00   
## Triglops.pingelii             1.00   
## Ammodytes.hexapterus          0.44   
## Nautichthys.pribilovius       0.99   
## Leptoclinus.maculatus         0.23   
## Limanda.aspera                0.99   
## Pleuronectidae                0.01 **
## Gadus.macrocephalus           0.06 . 
## Eumesogrammus.praecisus       0.40   
## Gymnelus.sp.                  0.06 . 
## Podothecus.veternus           1.00   
## Trichocottus.brashnikovi      1.00   
## Limanda.proboscidea           0.69   
## Aspidophoroides.monopterygius 1.00   
## Podothecus.accipenserinus     0.88   
## Hemilepidotus.papilio         1.00   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Contrast: Offshore_Northern 
## 
##                               average      sd   ratio     ava     avb cumsum
## Gymnocanthus.tricuspis        0.06051 0.02664 2.27130 5.74800 0.79100  0.109
## Lumpenus.fabricii.sp.         0.05940 0.03456 1.71870 6.51300 1.80100  0.215
## Hippoglossoides.robustus      0.04909 0.03155 1.55580 5.15600 1.53400  0.304
## Myoxocephalus.scorpius        0.04240 0.03156 1.34330 4.13800 1.34200  0.380
## Aspidophoroides.olrikii       0.03418 0.02735 1.24960 3.75300 2.17600  0.441
## Artediellus.scaber            0.02952 0.03000 0.98400 1.78900 2.01700  0.494
## Boreogadus.saida              0.02920 0.02750 1.06150 5.14800 5.39500  0.547
## Anisarchus.medius             0.02807 0.02631 1.06670 4.97000 4.90000  0.597
## Lycodes.spp.                  0.02738 0.02636 1.03860 4.15900 4.49100  0.646
## Liparis.spp.                  0.02394 0.02834 0.84480 1.56800 1.16000  0.689
## Gymnelus.hemifasciatus        0.02225 0.02756 0.80740 1.01700 1.49500  0.729
## Liparis.tunicatus             0.01969 0.02633 0.74760 1.43900 0.57900  0.764
## Icelus.spatula                0.01503 0.02306 0.65180 0.27200 1.24100  0.791
## Liparis.gibbus                0.01242 0.02050 0.60590 0.72700 0.56600  0.814
## Eleginus.gracilis             0.01192 0.02205 0.54080 0.96300 0.13400  0.835
## Triglops.pingelii             0.01094 0.01920 0.56960 0.78600 0.33200  0.855
## Aspidophoroides.monopterygius 0.01041 0.01970 0.52850 0.90600 0.08400  0.874
## Podothecus.veternus           0.00987 0.02017 0.48960 0.84800 0.05700  0.891
## Leptoclinus.maculatus         0.00895 0.01841 0.48620 0.48200 0.41500  0.907
## Gymnelus.viridis              0.00822 0.01805 0.45560 0.10900 0.69800  0.922
## Ammodytes.hexapterus          0.00774 0.01899 0.40770 0.24700 0.48000  0.936
## Gadus.chalcogrammus           0.00765 0.01732 0.44160 0.36100 0.37200  0.950
## Nautichthys.pribilovius       0.00673 0.01573 0.42770 0.43300 0.22400  0.962
## Stichaeus.punctatus           0.00503 0.01519 0.33140 0.43000 0.03200  0.971
## Limanda.aspera                0.00362 0.01313 0.27600 0.23200 0.10100  0.978
## Hemilepidotus.papilio         0.00265 0.01124 0.23600 0.21800 0.00000  0.982
## Eumesogrammus.praecisus       0.00252 0.01004 0.25150 0.09900 0.14500  0.987
## Limanda.proboscidea           0.00209 0.01073 0.19460 0.10700 0.07300  0.991
## Gadus.macrocephalus           0.00177 0.00931 0.19010 0.14900 0.00000  0.994
## Trichocottus.brashnikovi      0.00173 0.00827 0.20930 0.14000 0.02900  0.997
## Gymnelus.sp.                  0.00109 0.00728 0.14910 0.03000 0.07200  0.999
## Pleuronectidae                0.00047 0.00493 0.09450 0.03700 0.00000  1.000
## Podothecus.accipenserinus     0.00023 0.00272 0.08470 0.00000 0.02500  1.000
##                                  p   
## Gymnocanthus.tricuspis        0.01 **
## Lumpenus.fabricii.sp.         0.01 **
## Hippoglossoides.robustus      0.01 **
## Myoxocephalus.scorpius        0.17   
## Aspidophoroides.olrikii       0.17   
## Artediellus.scaber            0.99   
## Boreogadus.saida              1.00   
## Anisarchus.medius             1.00   
## Lycodes.spp.                  1.00   
## Liparis.spp.                  1.00   
## Gymnelus.hemifasciatus        0.13   
## Liparis.tunicatus             0.56   
## Icelus.spatula                0.48   
## Liparis.gibbus                0.83   
## Eleginus.gracilis             1.00   
## Triglops.pingelii             1.00   
## Aspidophoroides.monopterygius 0.01 **
## Podothecus.veternus           0.88   
## Leptoclinus.maculatus         0.01 **
## Gymnelus.viridis              0.93   
## Ammodytes.hexapterus          0.91   
## Gadus.chalcogrammus           0.94   
## Nautichthys.pribilovius       1.00   
## Stichaeus.punctatus           1.00   
## Limanda.aspera                1.00   
## Hemilepidotus.papilio         0.69   
## Eumesogrammus.praecisus       0.98   
## Limanda.proboscidea           0.92   
## Gadus.macrocephalus           0.99   
## Trichocottus.brashnikovi      1.00   
## Gymnelus.sp.                  0.98   
## Pleuronectidae                1.00   
## Podothecus.accipenserinus     1.00   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation: free
## Number of permutations: 99

Look at diversity metrics (RaoQ and 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_tax
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    :133  
##  1st Qu.:39.90   1st Qu.: 0.328   1st Qu.:32.06   Intermediate:170  
##  Median :43.00   Median : 2.140   Median :32.28   Northern    :128  
##  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.56                     
##        tax.cluster 
##  Nearshore   : 86  
##  Intermediate: 95  
##  Offshore    :111  
##  Northern    :139  
##                    
## 
arctic_env_fd$raoQ<-as.numeric(fd.chukchi$RaoQ)

simpson<-as.data.frame(diversity(arctic_cpue, "simpson"))
names(simpson)<-"simpson"
arctic_env_fd$simpson<-simpson$simpson

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: iterators
library(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 3084.3                                
## 2    430 3181.0 -2   -96.714 6.7103 0.001351 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
fd.aov<-aov(lm.fd)
summary(lm.fd)
## 
## Call:
## lm(formula = raoQ ~ cluster, data = arctic_env_fd)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -12.9821  -1.5045   0.3644   1.8041   6.3515 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          20.3232     0.2328  87.309  < 2e-16 ***
## clusterIntermediate  -0.2099     0.3108  -0.675 0.499822    
## clusterNorthern      -1.1367     0.3324  -3.420 0.000687 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.684 on 428 degrees of freedom
## Multiple R-squared:  0.0304, Adjusted R-squared:  0.02587 
## F-statistic:  6.71 on 2 and 428 DF,  p-value: 0.001351
TukeyHSD(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.2098722 -0.9407545  0.5210101 0.7779705
## Northern-Southern     -1.1366887 -1.9184350 -0.3549424 0.0019820
## Northern-Intermediate -0.9268165 -1.6656645 -0.1879685 0.0093702
library(rcompanion)
## Warning: package 'rcompanion' was built under R version 4.3.3
Sum = groupwiseMean(raoQ ~ cluster,
                    data   = arctic_env_fd,
                    conf   = 0.95,
                    digits = 3)
clus.cols.fuzzy <- c("#FD1C26","#98A100", "#0D7EF9")
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
lm.tax<-lm(simpson~ tax.cluster, 
           data=arctic_env_fd)
summary(lm.tax)
## 
## Call:
## lm(formula = simpson ~ tax.cluster, data = arctic_env_fd)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.32198 -0.01993  0.00883  0.03188  0.11637 
## 
## Coefficients:
##                          Estimate Std. Error t value Pr(>|t|)    
## (Intercept)              0.884585   0.006214 142.344  < 2e-16 ***
## tax.clusterIntermediate -0.075286   0.008578  -8.777  < 2e-16 ***
## tax.clusterOffshore      0.005333   0.008279   0.644     0.52    
## tax.clusterNorthern     -0.062609   0.007906  -7.919  2.1e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.05763 on 427 degrees of freedom
## Multiple R-squared:  0.278,  Adjusted R-squared:  0.2729 
## F-statistic:  54.8 on 3 and 427 DF,  p-value: < 2.2e-16
tax.aov<-aov(lm.tax)
summary(tax.aov)
##              Df Sum Sq Mean Sq F value Pr(>F)    
## tax.cluster   3  0.546 0.18199    54.8 <2e-16 ***
## Residuals   427  1.418 0.00332                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
TukeyHSD(tax.aov)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = lm.tax)
## 
## $tax.cluster
##                                diff          lwr         upr     p adj
## Intermediate-Nearshore -0.075285514 -0.097409213 -0.05316182 0.0000000
## Offshore-Nearshore      0.005332905 -0.016019749  0.02668556 0.9175544
## Northern-Nearshore     -0.062609293 -0.083001494 -0.04221709 0.0000000
## Offshore-Intermediate   0.080618419  0.059843479  0.10139336 0.0000000
## Northern-Intermediate   0.012676221 -0.007110243  0.03246269 0.3504801
## Northern-Offshore      -0.067942198 -0.086862610 -0.04902179 0.0000000
Sum.simpson = groupwiseMean(simpson ~ tax.cluster,
                    data   = arctic_env_fd,
                    conf   = 0.95,
                    digits = 3)
clus.cols.tax <- c("#FD1C26","#98A100", "#0D7EF9","orange")

Sum.simpson$tax.cluster  <- factor(Sum.simpson$tax.cluster ,
                                            levels = c("Nearshore",
                                                       "Intermediate",
                                                       "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)+
  theme_bw()+
  labs(x="Cluster", y="Simpson Diversity Index")
tax.cluster.CI

Diversity metrics compared to environmental gradients for functional diversity

arctic.cwm.rda.sites$raoQ<-as.numeric(arctic_env_fd$raoQ)

lm.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)
## 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 ' ' 1
summary(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-16
lm.cca.step<-stepAIC(lm.rda.1)
## Start:  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
raoQ.rda.1<-
  ggplot()+
  geom_smooth(data=arctic.cwm.rda.sites,
              aes(x=RDA1, y=raoQ), formula = y~x,
              method = "loess")+
  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~x,
              method = "loess")+
  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

Diversity metrics compared to environmental gradients for taxonomic diversity

##taxonomy
arctic.cwm.cca.sites.tax$simpson<-as.numeric(arctic_env_fd$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 ' ' 1
summary(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-16
lm.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.2
summary(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-16
arctic.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~x,
              method = "loess" )+
  geom_point(data=arctic.cwm.cca.sites.tax,
             aes(x=CCA1, y=simpson, color=cluster))+
  scale_color_manual(values=clus.cols)+
  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~x,
              method = "loess")+
  geom_point(data=arctic.cwm.cca.sites.tax,
             aes(x=CCA2, y=simpson, color=cluster))+
  scale_color_manual(values=clus.cols)+
  theme_bw()+
  labs(x="CCA2", y="Simpson Diversity Index")
simpson.cca.2

Supplementary Figures

arctic_env_fd_2<-arctic_env_fd
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, tax.cluster as id variables
View(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)
supp_fig_1