DZ_reanalysis_MM_data

# Mia Momberg - reanalysis of JEcol dataset

# importing libraries
library (readxl)
library (tidyverse)
library (modEvA)
library (car)

# importing data from Dryad repository
url <- 'https://datadryad.org/stash/downloads/file_stream/573456'
destfile <- tempfile ()
download.file (url, destfile, mode = 'wb')
mm <- read_excel (destfile, sheet = 'Momberg_et_al_data')

# preparing data for analysis
spe <- as.data.frame (mm[,12:29])  # species composition data
env <- mm[, 1:11]                  # environmental data for each quadrat
env$Grid_fact <- as.factor (env$Grid)
env$Grid_sample <- as.factor (sample (env$Grid))

# visual description of wind stress
ggplot (env, aes (x = Grid_fact, y = Max_wind_stress)) + 
  geom_boxplot () + 
  xlab ("Grid") + ylab ("Observed wind")

# GLZ - the same as paper:
GLZ_simple <- glm (Vascular_plant_cover/100 ~ poly (PDIR, 2) + poly (Soil_depth, 2) + poly (Rock_cover, 2) + poly (T_Jan, 2) + poly (T_June, 2) + poly (M_Oct, 2) + poly (pH, 2), data = env, family = quasibinomial ())
GLZ_full <- update (GLZ_simple, .~ . + poly (Max_wind_stress, 2))
anova (GLZ_simple, GLZ_full, test = 'F')  # P = 8.397e-15 ***

# deviance explained by the full model
Dsquared(GLZ_full) # 0.669, same as 66.9 reported in Table 3
D2_increment_obs <- Dsquared(GLZ_full)-Dsquared(GLZ_simple) # 0.01385595 increase of deviance by adding the wind to the model

# preparing artificial wind variable: 1) calculate mean wind stress per grid; 2) center values for wind stress of each quadrat by subtracting the mean wind stress value of given grid; 3) randomize centered quadrat values among the whole dataset; 4) to each randomized quadrat value of given grid add the mean value of another grid (ie "permute" mean wind stress values among grids)
mean_wind <- env %>% group_by (Grid_fact) %>% summarize (mean_wind = mean (Max_wind_stress)) %>% as.data.frame ()
Grid_rand <- data.frame (sapply (1:100, FUN = function (x) sample (1:9)))

env <- env %>% mutate (Max_wind_stress_center = Max_wind_stress - mean_wind[Grid, 'mean_wind'])

env_perm <- as.data.frame (matrix (nrow = 1440, ncol = 100))
for (i in 1:100)
{
  Grid_rand <- sample (1:9)
  env_perm [,i] <- sample(env$Max_wind_stress_center) + mean_wind[Grid_rand[env$Grid], 'mean_wind']
}

# display one of them
ggplot (env_perm, aes (x = env$Grid_fact, y = V1)) + 
  geom_boxplot () + 
  xlab ("Grid") + ylab ("Permuted wind")

# GLZ with "centred" wind stress
GLZ_full_centered <- update (GLZ_simple, .~ . + poly (Max_wind_stress_center, 2))
anova (GLZ_simple, GLZ_full_centered, test = 'F')  # P = 0.0664 - marginally significant! Looks like topographical estimation of wind has some weak effect
D2_increment_centered <- Dsquared(GLZ_full_centered)-Dsquared (GLZ_simple) # 0.00121418 but D2 increment is veeery low

# ... and use them to calculate GLZ:
recoded_results <- as.data.frame (matrix (ncol = 2, nrow = 100, dimnames = list (NULL, c('P_val', 'D2_increment'))))
for (i in 1:100)
  {
  GLZ_full_recoded <- update (GLZ_simple, .~ . + poly (env_perm[,i], 2))
  recoded_results[i,1] <- anova (GLZ_simple, GLZ_full_recoded, test = 'F')$`Pr(>F)`[2]
  recoded_results[i,2] <- Dsquared(GLZ_full_recoded)-Dsquared(GLZ_simple)
  }
  
sum (recoded_results$P_val < 0.001)  # 89 significant results out of 100, indicating sharply inflated Type I error rate
hist (recoded_results$D2_increment)
abline (v = 0.01385595, col = 'red')

sum (recoded_results$D2_increment >= D2_increment_obs)  # in 45 cases the incremental deviance by permuted wind is higher than that of observed wind

# Permutation P-value:
sum (c(D2_increment_obs, recoded_results$D2_increment) >= D2_increment_obs)/100 # 0.46 - not significantly higher than permuted wind stress

# Re-do Vegetation analysis
library (vegan)
spe_all <- spe[rowSums (spe)>0,]
env_all <- env[rowSums (spe)>0,]
env_perm_all <- env_perm[rowSums (spe)>0, ]

NMDS <- metaMDS (spe_all, autotransform = TRUE, trymax = 5)
ordiplot (NMDS, display = 'si', type = 'n')
points (NMDS, col = env_all$Grid, pch = 16)

ef_obs <- envfit (NMDS, env_all$Max_wind_stress, permutations = 1999)  # highly significant, R2 = 0.0227
ef_cent <- envfit (NMDS, env_all$Max_wind_stress_center, permutations = 1999) # not significant, R2 = 0.0016
ef_perm <- envfit (NMDS, env_perm_all, permutations = 1999)
sum (ef_perm$vectors$pvals < 0.001)  # 91 out of 100 significant at 0.001
sum (ef_perm$vector$r > ef_obs$vectors$r) # in 84 cases R2 by perm wind is higher than by observed wind