pkgBuild/manuscript/manuscript_stats.R

In rBatt/trawlDiversity: Long-term changes in species richness along the North American coastline

library('trawlDiversity')
library('rbLib')
library('lme4')
library('car')
library('multcomp')

setwd("~/Documents/School&Work/pinskyPost/trawlDiversity")

# ========================
# = Richness Time Series =
# ========================
# ---- mean richness ----
# smallest and largest long-term averages
comm_master[,mean(reg_rich), by='reg'][reg%in%c("gmex","shelf")]

# ---- long-term variability ----
comm_master[,stats::sd(reg_rich),by='reg']#[,mean(V1)]
comm_master[,stats::sd(reg_rich),by='reg'][,mean(V1)]

# ---- long-term trend ----
load("pkgBuild/results/rich_trend_kendall.RData")
load("pkgBuild/results/rich_naive_trend_kendall.RData")
rich_trend_kendall[reg!="wcann",BH:=p.adjust(pvalue, method="BH")]
rich_trend_kendall <- rich_trend_kendall[reg!="wcann",list(reg=reg, estimate=tau, BH=BH, p.value=pvalue)]
rich_naive_trend_kendall[reg!="wcann",BH:=p.adjust(taup, method='BH')]
rich_naive_trend_kendall <- rich_naive_trend_kendall[reg!="wcann",list(reg=reg, estimate=tau, BH=BH, p.value=taup)]
print(rich_naive_trend_kendall)
print(rich_trend_kendall)


# =======================================
# = Prevalence and Proximity to Absence =
# =======================================
r2.corr.mer <- function(m) {
   lmfit <-  lm(model.response(model.frame(m)) ~ fitted(m))
   summary(lmfit)$r.squared
}
omega2 <- function(m){
	1-var(residuals(m))/(var(model.response(model.frame(m))))
}

# ---- testing on ebs and column that i think might be relevant ----
# b <- spp_master[reg=='ebs' & !is.na(stretch_type) & propStrata!=0, list(spp, year, propStrata, ext_dist_sign, stretch_type, event_year, ce_categ)]

# too complicated to think about stretch_type while making sure i have less-obvious components right
# (tb <- lmer(propStrata~stretch_type*ext_dist_sign+(ext_dist_sign|spp:stretch_type), data=b))
# summary(tb)
#
# (tb_pc <- lmer(propStrata~ext_dist_sign+(ext_dist_sign|spp), data=b[stretch_type=="post_col"]))
# summary(tb_pc)

big <- spp_master[!is.na(stretch_type) & propStrata!=0]
big <- big[,list(prevalence=propStrata, time=ext_dist, type=as.character(stretch_type), reg=reg, event=as.character(event_year), spp=spp)]

# LMER Model Components:
# Begin Model

# mm_form <- formula(
# 	# Repsonse Variable
# 	prevalence ~
#
# 	# Fixed Effects
# 	time + time:type +
# 	# ext_dist_sign*stretch_type +  # stretch type has no discernable effect on the intercept, which makes sense, because 0 years for either colonization or extinction means 0% strata occupied!
#
# 	 # Random Effects
# 	(1 + time + type:time | reg) + (1 + time + type:time | reg:spp) # + (1 | reg:event))
# )

# mm_form <- formula(
# 	# Repsonse Variable
# 	prevalence ~
#
# 	# Fixed Effects
# 	time*reg +
#
# 	 # Random Effects
# 	(1 + time | spp)
# )

mm_form <- formula(
	# Repsonse Variable
	prevalence ~

	# Fixed Effects
	time + 

	 # Random Effects
	(1 + time | reg) + (1 + time | reg:spp) # + (1 | reg:event))
)

# am <- afex::mixed(mm_form,data=big)
(prev_prox_mod <- lmer(mm_form,data=big))

# LMER Model Metrics
summary(prev_prox_mod) # summary :p
r2.corr.mer(prev_prox_mod) # R^2, v1
omega2(prev_prox_mod) # Omega_naught_squared (R^2 v2)
Anova(prev_prox_mod)
anova(prev_prox_mod)
coef(prev_prox_mod)$reg


# ---- Meaning of intercept in prevalence - proximity model ----
spp_master[ce_categ!='neither' & propStrata!=0,min(propStrata),by=c("spp","reg")][,summary(V1)]
spp_master[ce_categ!='neither' & propStrata!=0,min(propStrata),by=c("spp","reg")][,ecdf(V1)(0.07)]


# ======================
# = Spatial Statistics =
# ======================


# ==========================
# = Within-site prevalence =
# ==========================
 # ---- variation among species vs among years ----
 prevSY_dat <- spp_master[present==1, list(reg=reg, spp=spp, year=as.character(year), propTow_occ=propTow_occ)]
(prevSY_mod <- lmer(propTow_occ ~ (1|reg/spp) + (1|reg:year), data=prevSY_dat))
summary(prevSY_mod)
omega2(prevSY_mod)
(re_sd <- c(Residual_stddev=summary(prevSY_mod)$sigma, sapply(summary(prevSY_mod)$varcor, attr, which="stddev")))
 
 # ---- maybe correlate within-site prevalence with biomass, as supplemental? ----
 
 # ---- correlated within-site prevalence with detectability (for Discussion) ----
 
 # ---- correlation between within-site prevalence and species richness ----
(wsprev_rich_mod <- lmer(reg_rich ~ propTow_occ_avg + (propTow_occ_avg | reg) ,data=comm_master))
summary(wsprev_rich_mod)
Anova(wsprev_rich_mod)
omega2(wsprev_rich_mod)