CorrectionModel.R
Tristan Blechinger
2023-02-07
library(ggplot2)## Warning: package 'ggplot2' was built under R version 4.2.3library(StanHeaders)## Warning: package 'StanHeaders' was built under R version 4.2.3library(rstan)## Warning: package 'rstan' was built under R version 4.2.3##
## rstan version 2.26.21 (Stan version 2.26.1)## For execution on a local, multicore CPU with excess RAM we recommend calling
## options(mc.cores = parallel::detectCores()).
## To avoid recompilation of unchanged Stan programs, we recommend calling
## rstan_options(auto_write = TRUE)
## For within-chain threading using `reduce_sum()` or `map_rect()` Stan functions,
## change `threads_per_chain` option:
## rstan_options(threads_per_chain = 1)## Do not specify '-march=native' in 'LOCAL_CPPFLAGS' or a Makevars filelibrary(brms)## Loading required package: Rcpp## Loading 'brms' package (version 2.19.0). Useful instructions
## can be found by typing help('brms'). A more detailed introduction
## to the package is available through vignette('brms_overview').##
## Attaching package: 'brms'## The following object is masked from 'package:rstan':
##
## loo## The following object is masked from 'package:stats':
##
## arlibrary(tidybayes)##
## Attaching package: 'tidybayes'## The following objects are masked from 'package:brms':
##
## dstudent_t, pstudent_t, qstudent_t, rstudent_tlibrary(tidyr)##
## Attaching package: 'tidyr'## The following object is masked from 'package:rstan':
##
## extractlibrary(bayestestR)##
## Attaching package: 'bayestestR'## The following object is masked from 'package:tidybayes':
##
## hdilibrary(bayesplot)## This is bayesplot version 1.10.0## - Online documentation and vignettes at mc-stan.org/bayesplot## - bayesplot theme set to bayesplot::theme_default()## * Does _not_ affect other ggplot2 plots## * See ?bayesplot_theme_set for details on theme setting##
## Attaching package: 'bayesplot'## The following object is masked from 'package:brms':
##
## rhatlibrary(cowplot)
library(dplyr)## Warning: package 'dplyr' was built under R version 4.2.3##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(ggpubr)##
## Attaching package: 'ggpubr'## The following object is masked from 'package:cowplot':
##
## get_legendlibrary(kableExtra)##
## Attaching package: 'kableExtra'## The following object is masked from 'package:dplyr':
##
## group_rowslibrary(loo)## Warning: package 'loo' was built under R version 4.2.3## This is loo version 2.5.1## - Online documentation and vignettes at mc-stan.org/loo## - As of v2.0.0 loo defaults to 1 core but we recommend using as many as possible. Use the 'cores' argument or set options(mc.cores = NUM_CORES) for an entire session.## - Windows 10 users: loo may be very slow if 'mc.cores' is set in your .Rprofile file (see https://github.com/stan-dev/loo/issues/94).##
## Attaching package: 'loo'## The following object is masked from 'package:rstan':
##
## loo#read in stable isotope data and scale/center
isodata <- read.csv("isotopedata.csv")
summary(isodata)## Lake DOW Site Fd.13.C
## Length:23 Min. : 1006200 Length:23 Min. :-37.23
## Class :character 1st Qu.:21014500 Class :character 1st Qu.:-31.95
## Mode :character Median :49007900 Mode :character Median :-29.73
## Mean :40531987 Mean :-30.82
## 3rd Qu.:49012700 3rd Qu.:-28.27
## Max. :86009000 Max. :-23.96
## Ed.13.C Fd.15.N Ed.15.N
## Min. :-35.12 Min. : 3.240 Min. : 3.070
## 1st Qu.:-29.95 1st Qu.: 3.940 1st Qu.: 4.105
## Median :-27.63 Median : 6.470 Median : 6.110
## Mean :-29.13 Mean : 7.089 Mean : 6.916
## 3rd Qu.:-27.36 3rd Qu.: 8.450 3rd Qu.: 8.255
## Max. :-24.28 Max. :14.440 Max. :12.750#adjust data by centering
#center data
isodata$Ed.13.Ccentered <- scale(isodata$Ed.13.C, center = T, scale = F)
isodata$Fd.13.Ccentered <- scale(isodata$Fd.13.C, center = T, scale = F)
isodata$Ed.15.Ncentered <- scale(isodata$Ed.15.N, center = T, scale = F)
isodata$Fd.15.Ncentered <- scale(isodata$Fd.15.N, center = T, scale = F)
#save values used to center fresh data
attrFC <- attributes(isodata$Fd.13.Ccentered)
attrFN <- attributes(isodata$Fd.15.Ncentered)
# plot(isodata$Ed.13.Ccentered, isodata$Fd.13.Ccentered)
# plot(isodata$Ed.13.C, isodata$Fd.13.C)#run models in BRMS using scaled data #model output available as RDS file
# set.seed(1)
# brmCcentered <- brm(Fd.13.Ccentered ~ 0 + Ed.13.Ccentered + (1 + Ed.13.Ccentered|Lake), isodata, family = gaussian, iter = 100000)
# brmNcentered <- brm(Fd.15.Ncentered ~ 0 + Ed.15.Ncentered + (1 + Ed.15.Ncentered|Lake), isodata, family = gaussian, iter = 100000)
# read in model output
brmCcentered <- readRDS("carbon_model_centered.rds")
brmNcentered <- readRDS("nitrogen_model_centered.rds")
#summary(brmCcentered, waic = TRUE)
#summary(brmNcentered, waic = TRUE)
#model output
#summary of all carbon model parameters
brmCcentered %>%
summarise_draws() #90% slope CI: 0.94 - 1.38 (1.158)## # A tibble: 17 × 10
## variable mean median sd mad q5 q95 rhat ess_bulk
## <chr> <num> <num> <num> <num> <num> <num> <num> <num>
## 1 b_Ed.13.Ccen… 1.16e+0 1.16e+0 0.137 0.104 0.940 1.38 1.00 15319.
## 2 sd_Lake__Int… 5.02e-1 4.00e-1 0.432 0.319 0.0443 1.32 1.00 8208.
## 3 sd_Lake__Ed.… 2.18e-1 1.55e-1 0.218 0.150 0.0133 0.652 1.00 3219.
## 4 cor_Lake__In… 3.23e-2 4.45e-2 0.583 0.754 -0.892 0.913 1.00 22841.
## 5 sigma 5.67e-1 5.54e-1 0.104 0.0965 0.422 0.757 1.00 5221.
## 6 r_Lake[Alexa… -8.18e-2 -2.47e-2 0.416 0.252 -0.808 0.505 1.00 11589.
## 7 r_Lake[Big.S… -6.85e-2 -1.54e-2 0.500 0.277 -0.915 0.624 1.00 6226.
## 8 r_Lake[Buffa… -3.34e-1 -2.86e-1 0.383 0.357 -0.999 0.151 1.00 3261.
## 9 r_Lake[Chipp… -8.43e-3 -3.59e-3 0.198 0.165 -0.341 0.318 1.00 66696.
## 10 r_Lake[Shami… 2.69e-1 2.19e-1 0.396 0.324 -0.234 0.926 1.00 8348.
## 11 r_Lake[Alexa… 2.87e-2 1.04e-2 0.149 0.0882 -0.189 0.294 1.00 14253.
## 12 r_Lake[Big.S… -2.84e-2 -1.30e-2 0.143 0.0885 -0.271 0.188 1.00 9025.
## 13 r_Lake[Buffa… -8.87e-2 -4.68e-2 0.236 0.135 -0.509 0.220 1.00 3355.
## 14 r_Lake[Chipp… 6.08e-3 1.30e-3 0.162 0.0946 -0.252 0.272 1.00 19413.
## 15 r_Lake[Shami… 9.00e-2 4.70e-2 0.213 0.122 -0.175 0.475 1.00 7859.
## 16 lprior -4.45e+0 -4.43e+0 0.0938 0.0313 -4.60 -4.40 1.00 5707.
## 17 lp__ -4.08e+1 -4.05e+1 3.72 3.61 -47.4 -35.3 1.00 4098.
## # ℹ 1 more variable: ess_tail <num>#summary of all nitrogen model parameters
brmNcentered %>%
summarise_draws() #90% slope CI: 0.89 - 1.28 (1.077)## # A tibble: 17 × 10
## variable mean median sd mad q5 q95 rhat ess_bulk
## <chr> <num> <num> <num> <num> <num> <num> <num> <num>
## 1 b_Ed.15.Ncen… 1.08e+0 1.08e+0 0.117 0.0825 8.95e-1 1.28 1.00 1011.
## 2 sd_Lake__Int… 3.85e-1 2.66e-1 0.435 0.244 2.43e-2 1.10 1.00 1391.
## 3 sd_Lake__Ed.… 1.78e-1 1.05e-1 0.245 0.103 9.06e-3 0.592 1.00 724.
## 4 cor_Lake__In… 5.29e-2 6.68e-2 0.582 0.754 -8.89e-1 0.922 1.00 20212.
## 5 sigma 5.55e-1 5.42e-1 0.0970 0.0889 4.22e-1 0.732 1.00 33324.
## 6 r_Lake[Alexa… -1.35e-1 -7.17e-2 0.411 0.210 -8.03e-1 0.330 1.00 1057.
## 7 r_Lake[Big.S… 1.09e-1 6.24e-2 0.224 0.163 -1.93e-1 0.528 1.00 4499.
## 8 r_Lake[Buffa… 6.02e-3 3.93e-3 0.458 0.187 -5.79e-1 0.638 1.01 504.
## 9 r_Lake[Chipp… -1.28e-1 -9.13e-2 0.199 0.165 -4.93e-1 0.143 1.00 16821.
## 10 r_Lake[Shami… 3.09e-2 1.88e-2 0.367 0.186 -5.30e-1 0.592 1.00 1006.
## 11 r_Lake[Alexa… 4.48e-2 2.07e-2 0.171 0.0814 -1.74e-1 0.332 1.00 916.
## 12 r_Lake[Big.S… -1.34e-2 -3.91e-3 0.139 0.0735 -2.46e-1 0.197 1.00 1556.
## 13 r_Lake[Buffa… 1.87e-2 1.12e-2 0.130 0.0636 -1.47e-1 0.221 1.00 930.
## 14 r_Lake[Chipp… -5.14e-2 -1.55e-2 0.158 0.0777 -3.60e-1 0.140 1.00 1039.
## 15 r_Lake[Shami… -3.64e-2 -1.21e-2 0.156 0.0691 -2.60e-1 0.152 1.00 1004.
## 16 lprior -5.50e+0 -5.48e+0 0.0618 0.0106 -5.56e+0 -5.47 1.00 951.
## 17 lp__ -4.21e+1 -4.19e+1 3.66 3.59 -4.85e+1 -36.6 1.00 1244.
## # ℹ 1 more variable: ess_tail <num># plot model output
#plot(brmCcentered)
#plot(brmNcentered)#graph fixed effects
Cfe_only <- tibble(Ed.13.Ccentered = seq(min(isodata$Ed.13.Ccentered), max(isodata$Ed.13.Ccentered), length.out=500)) %>%
add_fitted_draws(brmCcentered, re_formula = NA, scale = "response", n = 500)## Warning: `fitted_draws` and `add_fitted_draws` are deprecated as their names were confusing.
## - Use [add_]epred_draws() to get the expectation of the posterior predictive.
## - Use [add_]linpred_draws() to get the distribution of the linear predictor.
## - For example, you used [add_]fitted_draws(..., scale = "response"), which
## means you most likely want [add_]epred_draws(...).
## NOTE: When updating to the new functions, note that the `model` parameter is now
## named `object` and the `n` parameter is now named `ndraws`.Cfe_only_mean <- Cfe_only %>%
group_by(Ed.13.Ccentered) %>%
summarise(.value = mean(.value))
Cfixed <- ggplot(Cfe_only, aes(x = Ed.13.Ccentered, y = .value)) +
stat_interval(alpha = 0.05) +
guides(color = FALSE) +
geom_line(data = Cfe_only_mean, color = "red", lwd = 1) +
geom_point(data = isodata, mapping = aes(x = Ed.13.Ccentered, y = Fd.13.Ccentered, shape = Lake), color = "black", size = 3.5, stroke = .7) +
scale_shape_manual(labels = c("Alexander", "Big Sandy", "Buffalo", "Chippewa", "Shamineau"),
values = c(2,8,7,21,3)) +
geom_abline(slope = 1, intercept = 0, lty = 2, lwd = 1) +
labs(x = expression(paste("EtOH ", delta^{13}, " C (\u2030)")), y = expression(paste("Fresh ", delta^{13}, " C (\u2030)"))) +
theme_linedraw() +
theme(axis.title=element_text(size=16,face="bold"))## Warning: The `<scale>` argument of `guides()` cannot be `FALSE`. Use "none" instead as
## of ggplot2 3.3.4.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.Nfe_only <- tibble(Ed.15.Ncentered = seq(min(isodata$Ed.15.Ncentered), max(isodata$Ed.15.Ncentered), length.out=500)) %>%
add_fitted_draws(brmNcentered, re_formula = NA, scale = "response", n = 500)## Warning: `fitted_draws` and `add_fitted_draws` are deprecated as their names were confusing.
## - Use [add_]epred_draws() to get the expectation of the posterior predictive.
## - Use [add_]linpred_draws() to get the distribution of the linear predictor.
## - For example, you used [add_]fitted_draws(..., scale = "response"), which
## means you most likely want [add_]epred_draws(...).
## NOTE: When updating to the new functions, note that the `model` parameter is now
## named `object` and the `n` parameter is now named `ndraws`.Nfe_only_mean <- Nfe_only %>%
group_by(Ed.15.Ncentered) %>%
summarise(.value = mean(.value))
Nfixed <- ggplot(Nfe_only, aes(x = Ed.15.Ncentered, y = .value)) +
stat_interval(alpha = 0.05) +
guides(color = FALSE) +
geom_line(data = Nfe_only_mean, color = "red", lwd = 1) +
geom_point(data = isodata, mapping = aes(x = Ed.15.Ncentered, y = Fd.15.Ncentered, shape = Lake), color = "black", size = 3.5, stroke = .7) +
scale_shape_manual(labels = c("Alexander", "Big Sandy", "Buffalo", "Chippewa", "Shamineau"),
values = c(2,8,7,21,3)) +
geom_abline(slope = 1, intercept = 0, lty = 2, lwd = 1) +
labs(x = expression(paste("EtOH ", delta^{15}, " N (\u2030)")), y = expression(paste("Fresh ", delta^{15}, " N (\u2030)"))) +
theme_linedraw() +
theme(axis.title=element_text(size=16,face="bold"))
CNfixed <- ggarrange(Cfixed, Nfixed,
nrow = 1, ncol = 2,
common.legend = TRUE, legend = "bottom")
CNfixed
#check group-level slope/intercept effects
# find region of practical equivalence
# rope(brmCcentered) #0.38
# rope(brmNcentered) #0.36
Cslope <- brmCcentered %>%
spread_draws(r_Lake[Lake,Ed.13.Ccentered]) %>%
ggplot(aes(y = Lake, x = r_Lake, fill = stat(abs(x) < .38))) +
stat_halfeye() +
geom_vline(xintercept = c(-.38, .38), linetype = "dashed") +
scale_fill_manual(values = c("gray80", "skyblue")) +
xlab("Carbon Random Slope Offset") +
theme_linedraw()+
theme(legend.position = "none") +
scale_x_continuous(limits = c(-1.5,1.5))
Nslope <- brmNcentered %>%
spread_draws(r_Lake[Lake,Ed.15.Ncentered]) %>%
ggplot(aes(y = Lake, x = r_Lake, fill = stat(abs(x) < .36))) +
stat_halfeye() +
geom_vline(xintercept = c(-.36, .36), linetype = "dashed") +
scale_fill_manual(values = c("gray80", "skyblue")) +
xlab("Nitrogen Random Slope Offset") +
theme_linedraw()+
theme(legend.position = "none") +
scale_x_continuous(limits = c(-1.5,1.5))
Cintercept <- brmCcentered %>%
spread_draws(r_Lake[Lake,Intercept]) %>%
ggplot(aes(y = Lake, x = r_Lake, fill = stat(abs(x) < .38))) +
stat_halfeye() +
geom_vline(xintercept = c(-.38, .38), linetype = "dashed") +
scale_fill_manual(values = c("gray80", "skyblue")) +
xlab("Carbon Random Intercept Offset") +
theme_linedraw()+
theme(legend.position = "none") +
scale_x_continuous(limits = c(-1.5,1.5))
Nintercept <- brmNcentered %>%
spread_draws(r_Lake[Lake,Intercept]) %>%
ggplot(aes(y = Lake, x = r_Lake, fill = stat(abs(x) < .36))) +
stat_halfeye() +
geom_vline(xintercept = c(-.36, .36), linetype = "dashed") +
scale_fill_manual(values = c("gray80", "skyblue")) +
xlab("Nitrogen Random Intercept Offset") +
theme_linedraw()+
theme(legend.position = "none") +
scale_x_continuous(limits = c(-1.5,1.5))
CNoffsets <- plot_grid(Cslope,
Nslope + theme(axis.text.y = element_blank(),
axis.ticks.y = element_blank(),
axis.title.y = element_blank()),
Cintercept,
Nintercept + theme(axis.text.y = element_blank(),
axis.ticks.y = element_blank(),
axis.title.y = element_blank()),
nrow = 2,
labels = NA)## Warning: `stat(abs(x) < 0.38)` was deprecated in ggplot2 3.4.0.
## ℹ Please use `after_stat(abs(x) < 0.38)` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.## Warning: Removed 10264 rows containing missing values (`stat_slabinterval()`).## Warning: Removed 9591 rows containing missing values (`stat_slabinterval()`).## Warning: Removed 10264 rows containing missing values (`stat_slabinterval()`).## Warning: Removed 9591 rows containing missing values (`stat_slabinterval()`).CNoffsets## Warning: Removed 1 rows containing missing values (`geom_text()`).
#hypothesis testing
#R-squared values
bayes_R2(brmCcentered) #0.980## Estimate Est.Error Q2.5 Q97.5
## R2 0.9795584 0.004056344 0.9693175 0.9843029bayes_R2(brmNcentered) #0.979## Estimate Est.Error Q2.5 Q97.5
## R2 0.9785942 0.003518459 0.9696729 0.9827415#set up hypothesis test based on fixed effect output
#fixef(brmCcentered)
#fixef(brmNcentered)
hypothesis(brmCcentered, "Ed.13.Ccentered > 1") #91% chance slope is greater than 1 ## Hypothesis Tests for class b:
## Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
## 1 (Ed.13.Ccentered)... > 0 0.16 0.14 -0.06 0.38 9.77
## Post.Prob Star
## 1 0.91
## ---
## 'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
## '*': For one-sided hypotheses, the posterior probability exceeds 95%;
## for two-sided hypotheses, the value tested against lies outside the 95%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.hypothesis(brmNcentered, "Ed.15.Ncentered > 1") #81% chance slope is greater than 1 ## Hypothesis Tests for class b:
## Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
## 1 (Ed.15.Ncentered)... > 0 0.08 0.12 -0.11 0.28 4.14
## Post.Prob Star
## 1 0.81
## ---
## 'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
## '*': For one-sided hypotheses, the posterior probability exceeds 95%;
## for two-sided hypotheses, the value tested against lies outside the 95%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.##test model accuracy #K-fold cross validation test for carbon
#kfoldC <- kfold(brmCcentered, K = 5, folds = 'grouped', group = "Lake", save_fits = TRUE)
#kfoldN <- kfold(brmNcentered, K = 5, folds = 'grouped', group = "Lake", save_fits = TRUE)
#read in k-fold cross validation results
kfoldC <- readRDS("kfoldcarboncentered.rds")
kfoldN <- readRDS("kfoldnitrogencentered.rds")
#define loss functions
rmse <- function(y, yrep) {
yrep_mean <- colMeans(yrep)
sqrt(mean((y - yrep_mean)^2))
}
rmspe <- function(y, yrep) {
yrep_mean <- colMeans(yrep)
sqrt(mean((y - yrep_mean)/y)^2)
}
#predict responses and evaluate loss
#carbon
kfpC <- kfold_predict(kfoldC)
rmse(y = kfpC$y, yrep = kfpC$yrep) #0.701## [1] 0.7027048rmspe(y = kfpC$y, yrep = kfpC$yrep) #2.94%## [1] 2.940255#nitrogen
kfpN <- kfold_predict(kfoldN)
rmse(y = kfpN$y, yrep = kfpN$yrep) #0.59## [1] 0.5900601rmspe(y = kfpN$y, yrep = kfpN$yrep) #6.01%## [1] 6.024614#apply correction to ethanol values and rescale
#save model slopes
Ccentered <- as.data.frame(fixef(brmCcentered))
Ncentered <- as.data.frame(fixef(brmNcentered))
#create new centered fresh values by multiplying by slopes and add to initial dataframe
isodata$newfreshcenteredC <- (isodata$Ed.13.Ccentered * Ccentered[1,1])
isodata$newfreshcenteredN <- (isodata$Ed.15.Ncentered * Ncentered[1,1])
#rescale new fresh values by adding mean fresh values
isodata$returnfreshC <- isodata$newfreshcenteredC + attrFC$'scaled:center'
isodata$returnfreshN <- isodata$newfreshcenteredN + attrFN$'scaled:center'
#values used to uncenter
mean(isodata$Fd.13.C) #-30.819## [1] -30.81913mean(isodata$Fd.15.N) #7.089## [1] 7.088696#testing that model fit is more accurate than no change (1:1 line)
#create fitted values on 1:1 line
checkC <- lm(returnfreshC ~ 0 + offset(1 * Ed.13.C), data = isodata)
checkN <- lm(returnfreshN ~ 0 + offset(1 * Ed.15.N), data = isodata)
#plots of original values(black), corrected values based on brms model(blue), and data fit on 1:1 line(red)
# carbon
plot(isodata$Ed.13.C, isodata$Fd.13.C, pch = 19)
abline(coef = c(0,1))
points(isodata$Ed.13.C, checkC$fitted.values, col = "red", pch = 19)
points(isodata$Ed.13.C, isodata$returnfreshC, col = "blue", pch = 19)
legend(x = -35.5, y = -24,legend = c("observed data", "1:1 fit", "model corrected data"), fill = c("black", "red", "blue"))
# nitrogen
plot(isodata$Ed.15.N, isodata$Fd.15.N, pch = 19)
abline(coef = c(0,1))
points(isodata$Ed.15.N, checkN$fitted.values, col = "red", pch = 19)
points(isodata$Ed.15.N, isodata$returnfreshN, col = "blue", pch = 19)
legend(x = 3, y = 14,legend = c("observed data", "1:1 fit", "model corrected data"), fill = c("black", "red", "blue"))
#RMSE for observed fresh values vs. 1:1 line
sqrt(sum((isodata$Fd.13.C - checkC$fitted.values)^2)/nrow(isodata)) #1.86 carbon## [1] 1.856626sqrt(sum((isodata$Fd.15.N - checkN$fitted.values)^2)/nrow(isodata)) #0.61 nitrogen## [1] 0.614205#RMSE for observed fresh values vs. model fit: model fit is more accurate
sqrt(sum((isodata$Fd.13.C - isodata$returnfreshC)^2)/nrow(isodata)) #0.61 carbon## [1] 0.6060523sqrt(sum((isodata$Fd.15.N - isodata$returnfreshN)^2)/nrow(isodata)) #0.53 nitrogen## [1] 0.5267476#marginal predictive check using LOO probability integral transform
#Carbon
#run leave-one-out cross validation
looC <- loo(brmCcentered, save_psis = TRUE)## Warning: Found 2 observations with a pareto_k > 0.7 in model 'brmCcentered'. It
## is recommended to set 'moment_match = TRUE' in order to perform moment matching
## for problematic observations.plot(looC) # diagnostic plot
yrepC <- posterior_predict(brmCcentered, na.rm = TRUE) # model output
psisC <- looC$psis_object
lwC <- weights(psisC, na.rm = TRUE) # normalized log weights
yC <- isodata$Fd.13.Ccentered[,1] # observed fresh zooplankton values
#plot results
loo_C <- suppressWarnings(ppc_loo_intervals(yC, yrepC, psis_object = psisC))
loo_C
# Most confident between ~-26 to -39 for C13, transformed to uncentered values
4 + attrFC$'scaled:center' # -26.819## [1] -26.81913-8 + attrFC$'scaled:center' # -38.819## [1] -38.81913#RMSE value for carbon (k-fold CV RMSE more representative of data, grouped by lake)
rmse(y = yC, yrep = yrepC) #.470## [1] 0.4700717rmspe(y = yC, yrep = yrepC) #2.314%## [1] 2.306374#Nitrogen
#run leave-one-out cross validation
looN <- loo(brmNcentered, save_psis = TRUE)
plot(looN)
yrepN <- posterior_predict(brmNcentered, na.rm = TRUE)
psisN <- looN$psis_object
lwN <- weights(psisN, na.rm = TRUE) # normalized log weights
yN <- isodata$Fd.15.Ncentered[,1] # this would be fresh zooplankton values for you
#plot results
loo_N <- suppressWarnings(ppc_loo_intervals(yN, yrepN, psis_object = psisN))
loo_N
# Most confident between ~2 to 15 for N15
7.5 + attrFN$'scaled:center' # 14.587## [1] 14.5887-5 + attrFN$'scaled:center' # 2.089## [1] 2.088696#RMSE value for nitrogen (k-fold CV RMSE more representative of data, grouped by lake)
rmse(y = yN, yrep = yrepN) #.469## [1] 0.4687541rmspe(y = yN, yrep = yrepN) # 4.994%## [1] 4.988865#plot combined results
loo_check <- plot_grid(loo_C + theme(legend.position="none"), loo_N + theme(legend.position="none"), align = 'h', labels = c("C13", "N15"), label_x = .05, ncol = 2, rel_widths = c(1,1))
legend <- get_legend(loo_C)
loo_check <- plot_grid(loo_check, legend, ncol = 2, rel_widths = c(1,.1))
loo_check
#table of model output and covergence
posterior <- as.data.frame(matrix(NA, nrow = 2, ncol = 5))
posterior[1,] <- diagnostic_posterior(brmCcentered)
posterior[2,] <- diagnostic_posterior(brmNcentered)
posterior[,5] <- c(0.980, 0.979)
posterior[,1] <- c("Carbon", "Nitrogen")
colnames(posterior) <- c("Slope Parameter", "R-hat", "ESS", "MCSE", "R^2")
output <- as.data.frame(matrix(NA, nrow = 2, ncol = 4))
output[1,] <- fixef(brmCcentered)
output[2,] <- fixef(brmNcentered)
output$slope <- c("Carbon", "Nitrogen")
colnames(output) <- c("Estimate", "Est.Error", "Q2.5", "Q97.5", "Slope Parameter")
col_order <- c("Slope Parameter", "Estimate", "Est.Error", "Q2.5", "Q97.5")
output <- output[,col_order]
combined <- bind_cols(output, posterior[,-1])
write.table(combined, file = "combinedtable.txt", sep = ',')
combined %>%
kbl(caption = "Model Output", digits = 3) %>%
kable_classic(full_width = F, html_font = "Cambria")| Slope Parameter | Estimate | Est.Error | Q2.5 | Q97.5 | R-hat | ESS | MCSE | R^2 |
|---|---|---|---|---|---|---|---|---|
| Carbon | 1.158 | 0.137 | 0.866 | 1.441 | 1.000 | 10381.910 | 0.001 | 0.980 |
| Nitrogen | 1.077 | 0.117 | 0.833 | 1.359 | 1.005 | 378.018 | 0.006 | 0.979 |
combine## function (...)
## {
## lifecycle::deprecate_warn("1.0.0", "combine()", "vctrs::vec_c()",
## always = TRUE)
## args <- list2(...)
## if (length(args) == 1 && is.list(args[[1]])) {
## args <- args[[1]]
## }
## args <- keep(args, function(.x) !is.null(.x))
## names(args) <- NULL
## if (length(args) == 0) {
## logical()
## }
## else {
## vec_c(!!!args)
## }
## }
## <bytecode: 0x000001db667a5918>
## <environment: namespace:dplyr>#hypothesis test table
hyptest <- as.data.frame(matrix(NA, nrow = 2, ncol = 7))
hyptest[1,] <- c("C13 Slope > 1", 0.16, 0.14, -0.06, 0.38, 9.77, 0.91)
hyptest[2,] <- c("N15 Slope > 1", 0.08, 0.12, -0.11, 0.28, 4.14, 0.81)
colnames(hyptest) <- c("Test", "Estimate", "Est.Error", "LCI", "UCI", "Evid. Ratio", "Probability")
write.table(hyptest, file = "hyptesttable.txt", sep = ',')
hyptest <- hyptest %>%
kbl(caption = "Hypothesis Test Results", digits = 3) %>%
kable_classic(full_width = F, html_font = "Cambria")
hyptest| Test | Estimate | Est.Error | LCI | UCI | Evid. Ratio | Probability |
|---|---|---|---|---|---|---|
| C13 Slope > 1 | 0.16 | 0.14 | -0.06 | 0.38 | 9.77 | 0.91 |
| N15 Slope > 1 | 0.08 | 0.12 | -0.11 | 0.28 | 4.14 | 0.81 |
sessionInfo()## R version 4.2.2 (2022-10-31 ucrt)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 19045)
##
## Matrix products: default
##
## locale:
## [1] LC_COLLATE=English_United States.utf8
## [2] LC_CTYPE=English_United States.utf8
## [3] LC_MONETARY=English_United States.utf8
## [4] LC_NUMERIC=C
## [5] LC_TIME=English_United States.utf8
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] loo_2.5.1 kableExtra_1.3.4 ggpubr_0.6.0
## [4] dplyr_1.1.1 cowplot_1.1.1 bayesplot_1.10.0
## [7] bayestestR_0.13.0 tidyr_1.3.0 tidybayes_3.0.3
## [10] brms_2.19.0 Rcpp_1.0.10 rstan_2.26.21
## [13] StanHeaders_2.26.21 ggplot2_3.4.2
##
## loaded via a namespace (and not attached):
## [1] colorspace_2.1-0 ggsignif_0.6.4 ellipsis_0.3.2
## [4] markdown_1.5 base64enc_0.1-3 rstudioapi_0.14
## [7] farver_2.1.1 svUnit_1.0.6 DT_0.27
## [10] lubridate_1.9.2 fansi_1.0.4 mvtnorm_1.1-3
## [13] xml2_1.3.3 bridgesampling_1.1-2 codetools_0.2-18
## [16] cachem_1.0.7 knitr_1.42 shinythemes_1.2.0
## [19] jsonlite_1.8.4 broom_1.0.3 ggdist_3.2.1
## [22] shiny_1.7.4 httr_1.4.5 compiler_4.2.2
## [25] backports_1.4.1 Matrix_1.5-1 fastmap_1.1.1
## [28] cli_3.6.1 later_1.3.0 htmltools_0.5.5
## [31] prettyunits_1.1.1 tools_4.2.2 igraph_1.4.1
## [34] coda_0.19-4 gtable_0.3.3 glue_1.6.2
## [37] reshape2_1.4.4 posterior_1.4.1 V8_4.2.2
## [40] carData_3.0-5 jquerylib_0.1.4 vctrs_0.6.1
## [43] svglite_2.1.1 nlme_3.1-160 crosstalk_1.2.0
## [46] insight_0.19.0 tensorA_0.36.2 xfun_0.38
## [49] stringr_1.5.0 ps_1.7.3 rvest_1.0.3
## [52] timechange_0.2.0 mime_0.12 miniUI_0.1.1.1
## [55] lifecycle_1.0.3 gtools_3.9.4 rstatix_0.7.2
## [58] zoo_1.8-11 scales_1.2.1 colourpicker_1.2.0
## [61] promises_1.2.0.1 Brobdingnag_1.2-9 parallel_4.2.2
## [64] inline_0.3.19 shinystan_2.6.0 yaml_2.3.7
## [67] curl_5.0.0 gridExtra_2.3 sass_0.4.5
## [70] stringi_1.7.12 highr_0.10 dygraphs_1.1.1.6
## [73] checkmate_2.1.0 pkgbuild_1.4.0 systemfonts_1.0.4
## [76] rlang_1.1.0 pkgconfig_2.0.3 matrixStats_0.63.0
## [79] distributional_0.3.2 evaluate_0.20 lattice_0.20-45
## [82] purrr_1.0.1 labeling_0.4.2 rstantools_2.3.0
## [85] htmlwidgets_1.6.2 processx_3.8.0 tidyselect_1.2.0
## [88] plyr_1.8.8 magrittr_2.0.3 R6_2.5.1
## [91] generics_0.1.3 pillar_1.9.0 withr_2.5.0
## [94] xts_0.13.0 datawizard_0.6.5 abind_1.4-5
## [97] tibble_3.2.1 crayon_1.5.2 car_3.1-1
## [100] arrayhelpers_1.1-0 utf8_1.2.3 rmarkdown_2.21
## [103] grid_4.2.2 callr_3.7.3 threejs_0.3.3
## [106] webshot_0.5.4 digest_0.6.31 xtable_1.8-4
## [109] httpuv_1.6.9 RcppParallel_5.1.7 stats4_4.2.2
## [112] munsell_0.5.0 viridisLite_0.4.1 bslib_0.4.2
## [115] shinyjs_2.1.0