Skip to content

Commit

Permalink
rm old commented code
Browse files Browse the repository at this point in the history
  • Loading branch information
@jsta
jsta committed Apr 23, 2020
1 parent 4ab2c9a commit 2637921
Showing 1 changed file with 0 additions and 228 deletions.
228 changes: 0 additions & 228 deletions scripts/03_model.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11,84 +11,6 @@ good_hu4s <- readRDS("data/dt.rds") %>%
dplyr::filter(n >= 3)
dt <- dplyr::filter(dt, hu4vzoneid %in% good_hu4s$hu4_zoneid)

# evaulate fixed effects
# (model_forms_fe <- list(
# "tp" = bf(tp ~ ag),
# "tp_nolulc" = bf(tp ~ maxdepth + iwslavratio +
# soilvorgvcarbon + wetlandvpotential + hu12vpptvmean +
# clayvpct + hu12vbaseflowvmean +
# nitrogenvfertilizervuse + nitrogenvlivestockvmanure +
# hu4vnitrogenvatmosphericvdeposition +
# phosphorusvfertilizervuse +
# phosphorusvlivestockvmanure),
# "tp_nuts" = bf(tp ~ nitrogenvfertilizervuse + nvinput +
# nitrogenvlivestockvmanure +
# hu4vnitrogenvatmosphericvdeposition +
# phosphorusvfertilizervuse + pvinput +
# phosphorusvlivestockvmanure),
# "tp_depth" = bf(tp ~ ag + maxdepth), # lake
# "tp_bf" = bf(tp ~ ag + maxdepth + hu12vbaseflowvmean), # transport
# "tp_nfert" = bf(tp ~ ag + maxdepth + hu12vbaseflowvmean +
# phosphorusvfertilizervuse), # sources
# "tp_buffer" = bf(tp ~ ag + maxdepth + hu12vbaseflowvmean +
# phosphorusvfertilizervuse + buffervcultivatedvcrops), # buffer
# "tn" = bf(tn ~ ag), # proxy
# "tn_nolulc" = bf(tn ~ maxdepth + iwslavratio +
# soilvorgvcarbon + wetlandvpotential + hu12vpptvmean +
# clayvpct + hu12vbaseflowvmean +
# nitrogenvfertilizervuse + nvinput + nitrogenvlivestockvmanure +
# hu4vnitrogenvatmosphericvdeposition +
# phosphorusvfertilizervuse + pvinput +
# phosphorusvlivestockvmanure),
# "tn_depth" = bf(tn ~ ag + maxdepth), # lake
# "tn_sc" = bf(tn ~ ag + maxdepth + soilvorgvcarbon), # transport
# "tn_nfert" = bf(tn ~ ag + maxdepth + soilvorgvcarbon +
# nitrogenvfertilizervuse), # sources
# "tn_buffer" = bf(tn ~ ag + maxdepth + soilvorgvcarbon +
# nitrogenvfertilizervuse + buffervcultivatedvcrops), # buffer
# "tn_ndep" = bf(tn ~ ag + maxdepth + soilvorgvcarbon +
# nitrogenvfertilizervuse + buffervcultivatedvcrops +
# hu4vnitrogenvatmosphericvdeposition)
# ))
#
# fe_brms <- list()
# i <- 3
# print(paste0("Fitting ", paste0("data/mcmc/fe/", names(model_forms_fe)[i])))
# fe_brms[[i]] <- brm(formula = model_forms_fe[[i]], data = dt,
# prior = set_prior("horseshoe(2)"), family = gaussian(),
# control = list(adapt_delta = 0.99))
# summary(fe_brms[[i]])

# saveRDS(re_brms[[i]], paste0("data/mcmc/re/", names(model_forms_re)[i]))
# }# else{
# re_brms[[i]] <- readRDS(paste0("data/mcmc/re/", names(model_forms_re)[i]))
# }
# }

# fe_brms <-
# lapply(seq_along(model_forms_fe), function(i)
# get_if_not_exists(brm_fit,
# paste0("data/mcmc/fe/", names(model_forms_fe)[i]),
# formula = model_forms_fe[[i]],
# data = dt))

# r2_fe <- dplyr::bind_rows(
# lapply(fe_brms, function(x) data.frame(brms::bayes_R2(x)))) %>%
# mutate(Model = names(model_forms_fe),
# Estimate = round(Estimate, 2)) %>%
# dplyr::select(Model, Estimate)

# r2_fe$`Proxy` <- c(rep("Ag", 5),
# rep("Ag", 5))
# r2_fe$`Lake` <- c(NA, rep("maxdepth", 4),
# NA, rep("maxdepth", 4))
# r2_fe$`Transport` <- c(NA, NA, "Baseflow", "Baseflow", "Baseflow",
# NA, NA, "Soil Org Carbon", "Soil Org Carbon", "Soil Org Carbon")
# r2_fe$`Source` <- c(NA, NA, NA, "P fertilizer", "P fertilizer",
# NA, NA, NA, "N fertilizer", "N fertilizer")
# r2_fe$`Buffer` <- c(NA, NA, NA, NA, "Row crop",
# NA, NA, NA, NA, "Row crop")

# evaulate spatial random effects
# {fixed effects} + 1/ag, 1/soybeans, 1/corn
(model_forms_re <- list(
Expand Down Expand Up @@ -233,51 +155,6 @@ re_brms <-
formula = model_forms_re[[i]],
data = dt))

# (model_forms_resuper <- list("tn_crops" = bf(tn ~
# maxdepth + iwslavratio +
# soilvorgvcarbon + wetlandvpotential + hu12vpptvmean +
# clayvpct + hu12vbaseflowvmean +
# nitrogenvfertilizervuse + nvinput + nitrogenvlivestockvmanure +
# hu4vnitrogenvatmosphericvdeposition +
# phosphorusvfertilizervuse + pvinput + phosphorusvlivestockvmanure +
# buffervcultivatedvcrops + buffervnatural +
# (1 + soybeans + corn + pasture | hu4vzoneid)),
# "tn_ag" = bf(tn ~
# maxdepth + iwslavratio +
# soilvorgvcarbon + wetlandvpotential + hu12vpptvmean +
# clayvpct + hu12vbaseflowvmean +
# nitrogenvfertilizervuse + nvinput + nitrogenvlivestockvmanure +
# hu4vnitrogenvatmosphericvdeposition +
# phosphorusvfertilizervuse + pvinput + phosphorusvlivestockvmanure +
# buffervcultivatedvcrops + buffervnatural +
# (1 + ag | hu4vzoneid)),
# "tp_crops" = bf(tp ~
# maxdepth + iwslavratio +
# soilvorgvcarbon + wetlandvpotential + hu12vpptvmean +
# clayvpct + hu12vbaseflowvmean +
# nitrogenvfertilizervuse + nvinput + nitrogenvlivestockvmanure +
# hu4vnitrogenvatmosphericvdeposition +
# phosphorusvfertilizervuse + pvinput + phosphorusvlivestockvmanure +
# buffervcultivatedvcrops + buffervnatural +
# (1 + soybeans + corn + pasture | hu4vzoneid)),
# "tp_ag" = bf(tp ~
# maxdepth + iwslavratio +
# soilvorgvcarbon + wetlandvpotential + hu12vpptvmean +
# clayvpct + hu12vbaseflowvmean +
# nitrogenvfertilizervuse + nvinput + nitrogenvlivestockvmanure +
# hu4vnitrogenvatmosphericvdeposition +
# phosphorusvfertilizervuse + pvinput + phosphorusvlivestockvmanure +
# buffervcultivatedvcrops + buffervnatural +
# (1 + ag | hu4vzoneid))
# ))
#
# resuper_brms <-
# lapply(seq_along(model_forms_resuper), function(i)
# get_if_not_exists(brm_fit,
# paste0("data/mcmc/resuper/", names(model_forms_resuper)[i]),
# formula = model_forms_resuper[[i]],
# data = dt))

# evaluate hu4 re slope significance
re_brms <- lapply(re_brms, function(x) get_re_signif(x))

Expand Down Expand Up @@ -424,108 +301,3 @@ dt %>%
#
# pairs(~maxdepth + hu12vbaseflowvmean + iwslavratio + , data = dt)
}

# ---- placeholder ----

# cat("predictors:
# lake: depth
# iws: area, ws:lk,
# pasture, wheat, corn, soybean, non-ag, other ag,
# org c, hydro conduc, runoff pot, erod,
# N fert, P fert, manure, N dep
# reponse:
# tn, tp")

# ---- exploratory_models ----
# #### Find optimal random effects structure
# library(lme4)
# library(ggeffects)

# model_form <- as.formula(tp ~ 1 + ag + (1 + ag | hu4vzoneid))
# fit0 <- lmer(model_form, data = dt, na.action = na.omit)
#
# model_form <- as.formula(tp ~ 1 + streamvcultivatedvcrops +
# (1 + streamvcultivatedvcrops | hu4vzoneid))
# fit1 <- lmer(model_form, data = dt, na.action = na.omit)
#
# sjstats::r2(fit0); sjstats::r2(fit1);
#
# model_form <- as.formula(tp ~ 1 +
# streamvcultivatedvcrops + maxdepth +
# (1 + streamvcultivatedvcrops | hu4vzoneid))
# fit2 <- lmer(model_form, data = dt, na.action = na.omit)
#
# # #### Fit unconditional model
# formula_unconditional <- as.formula(tp ~ (1|hu4_zoneid) - 1)
# fit <- lmer(formula_unconditional,
# data = dt, na.action = na.omit)
# getME(fit, "X")
# getME(fit, "Z")
#
# re <- ggpredict(fit, type = "re")
# sjstats::r2(fit)

# #####
#
# dt <- readRDS("data/dt.rds")
#
# #### explore mediation effects
# # https://en.wikipedia.org/wiki/Mediation_(statistics)
# fit1 <- lm(tn ~ row_crop_pct, data = dt)
# fit2 <- lm(tn ~ maxdepth, data = dt)
# fit3 <- lm(tn ~ row_crop_pct + maxdepth, data = dt)
#
# summary(fit1)
# summary(fit2)
# summary(fit3)
#
# ####
#
# plot(re)$hu4_zoneid +
# geom_hline(aes(yintercept = mean(
# dplyr::filter(dt, hu4_zoneid == "HU4_16")$tp))) +
# theme(axis.text.x = element_text(angle = 90))
#
# # shinyMer(fit)
#
# fit_ranef <- left_join(
# data.frame(ranef(fit), stringsAsFactors = FALSE),
# dplyr::select(dt, hu4_zoneid, geometry),
# by = c("grp" = "hu4_zoneid"))
# st_geometry(fit_ranef) <- fit_ranef$geometry
# fit_ranef <- mutate(fit_ranef,
# grp_id = stringr::str_extract(grp, "(?<=_)(\\d*)$"))
#
# # mapview::mapview(fit_ranef, zcol = "condval")
#
# dt_mean <- mean(dt$tp, na.rm = TRUE)
#
# test <- data.frame(dt) %>%
# group_by(hu4_zoneid) %>%
# summarize(tp_grp = mean(tp, na.rm = TRUE),
# tp_diff = tp_grp - dt_mean) %>%
# left_join(
# dplyr::select(data.frame(fit_ranef), grp, condval),
# by = c("hu4_zoneid" = "grp")) %>%
# distinct()
#
# mean(test$tp_grp)
#
# plot(test$tp_grp, test$condval)
# plot(test$tp_diff, test$condval)
# abline(0, 1)
# abline(h = dt_mean)
#
# #### Fit random intercept model with covariate model
# fit <- lmer(tp ~ 1 + iws_ag_2006_pcent + (1|hu6), data=dt, na.action=na.omit)
#
# #### Fit random intercept and slope with covariate
# fit <- lmer(tp ~ 1 + iws_ag_2006_pcent + (1+iws_ag_2006_pcent|hu6), data=dt, na.action=na.omit)
#
# # shinyMer(fit)
#
# fit_ranef <- left_join(data.frame(ranef(fit)),
# dplyr::select(crops, location_desc, geom),
# by = c("grp" = "location_desc"))
#
# # mapview(st_sf(fit_ranef), zcol = "condval")

0 comments on commit 2637921

Please sign in to comment.