exec/statistical_tests.R

In bmcnellis/RSFIA: FIA Tree Mortality Analysis

# old tests

mort_df <- FIA_mortality_with_explanatory
cols <- colnames(mort_df)[grep('ndead', colnames(mort_df))]
cols <- cols[-grep('of', cols)]
cols <- cols[-grep('Missing', cols)]

df <- matrix(nrow = nrow(mort_df), ncol = length(cols))
c <- 0
for (i in cols) {
  c <- c + 1
  df[, c] <- ifelse(mort_df[, i] == 0, 0, 1)
}
colnames(df) <- cols
df <- df[-which(rowSums(df) == 0), ]

df_veg <- vegdist(df, 'binomial', binary = T)
df_MDS <- metaMDS(df_veg, 'binomial')

#### gamlss

library(gamlss)
library(RSFIA)
library(ClelandEcoregions)

mort_df <- FIA_mortality_with_explanatory
IQR_groups <- CalcCategoricalMortByIQR(
  x = mort_df$mort_rate, group = mort_df[['Cleland_section']]
)[, 4]
mort_df$is_elevated <- IQR_groups
mort_df$sections <- as.factor(KeyClelandCode(mort_df$Cleland_section, 'section'))
mort_df <- mort_df[which(mort_df$is_elevated == F), ]
#mort_df <- mort_df[which(mort_df$dominant_AGENTCD == 'Weather'), ]
mort_df <- mort_df[, c('mort_rate', 'sections', 'MAT_10', 'MAP_10',
                       'BDTICM', 'CECSOL', 'BLDFIE', 'PHIHOX', 'CRFVOL',
                       'max_temp_by_Q_val_10', 'mean_dry_period_10', 'is_elevated',
                       'dominant_AGENTCD', 'Disease_frac_of_ndead',
                       'Vegetation_rate_wrt_ntree'
)]
mort_df <- na.omit(mort_df)

BEINF_mod <- gamlss(
  Vegetation_rate_wrt_ntree ~ sections + MAT_10 + BDTICM + CECSOL + BLDFIE + PHIHOX + CRFVOL + max_temp_by_Q_val_10 + mean_dry_period_10 + dominant_AGENTCD,
  data = mort_df, family = BEOI)
pred <- meanBEINF(BEINF_mod)
#pred <- meanBEZI(BEINF_mod)
#pred <- meanBEOI(BEINF_mod)
means_m1 <- lpred(BEINF_mod, type='response', what='mu', se.fit=T)
plot(mort_df$mort_rate, pred)
#plot(mort_df$mort_rate, means_m1[[1]])
fit <- lm(mort_df$mort_rate ~ pred)

### zoib
library(RSFIA)
library(ClelandEcoregions)
data("FIA_mortality_with_explanatory")
mort_df <- FIA_mortality_with_explanatory
mort_df$sections <- as.factor(KeyClelandCode(mort_df$Cleland_section, 'section'))
eg.fixed <- zoib(mort_rate ~ sections + MAT_10 | 1,
                 data=mort_df, joint = FALSE,  random = 0,
                 EUID = 1:nrow(d), zero.inflation = T,
                 one.inflation = T, n.iter = 1100, n.thin = 5,
                 n.burn=100)
coeff <- eg.fixed$coef
summary(coeff)
eg.fixed$Xb; eg.fixed$Xd; eg.fixed$Xb0; eg.fixed$Xb1
ypred = rbind(eg.fixed$ypred[[1]],eg.fixed$ypred[[2]])
post.mean= apply(ypred,2,mean)
plot(GasolineYield$yield, post.mean, col='blue',pch=2)
abline(0,1,col='red')