library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.2
## ✔ ggplot2   4.0.0     ✔ tibble    3.3.0
## ✔ lubridate 1.9.4     ✔ tidyr     1.3.1
## ✔ purrr     1.1.0     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(ggpubr)
library(ggplot2)

set.seed(123)

# --- Simulated dataset with LST and NDVI ---
dat <- tibble(
  group = rep(c("Urban", "NonUrban"), each = 30),
  LST   = c(rnorm(30, mean = 34.5, sd = 2),  # higher in Urban
            rnorm(30, mean = 30.0, sd = 2)),
  NDVI  = runif(60, 0.2, 0.8)  # random vegetation index
)

# --- Frequentist model 1: t-test as lm() ---
fit_group <- lm(LST ~ group, data = dat)
summary(fit_group)
## 
## Call:
## lm(formula = LST ~ group, data = dat)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.8390 -1.1272 -0.2254  1.1182  3.9812 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  30.3567     0.3327  91.257  < 2e-16 ***
## groupUrban    4.0491     0.4704   8.607 5.93e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.822 on 58 degrees of freedom
## Multiple R-squared:  0.5609, Adjusted R-squared:  0.5533 
## F-statistic: 74.08 on 1 and 58 DF,  p-value: 5.935e-12
# Plot: LST by land use group
ggboxplot(dat, x = "group", y = "LST",
          color = "group", palette = c("steelblue", "goldenrod"),
          add = "jitter", shape = "group") +
  theme_minimal() +
  labs(title = "LST by Land Use Group",
       x = "Land Use Group", y = "Land Surface Temperature (°C)")

# --- Frequentist model 2: regression LST ~ NDVI ---
fit_ndvi <- lm(LST ~ NDVI, data = dat)
summary(fit_ndvi)
## 
## Call:
## lm(formula = LST ~ NDVI, data = dat)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.3206 -2.1460  0.1752  1.9968  5.4897 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   31.180      1.187   26.27   <2e-16 ***
## NDVI           2.465      2.326    1.06    0.294    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.723 on 58 degrees of freedom
## Multiple R-squared:  0.019,  Adjusted R-squared:  0.002087 
## F-statistic: 1.123 on 1 and 58 DF,  p-value: 0.2936
# Plot: regression of LST on NDVI
ggscatter(dat, x = "NDVI", y = "LST",
          add = "reg.line",          # Adds regression line
          conf.int = TRUE,           # Adds CI ribbon
          cor.coef = TRUE,           # Show correlation R
          size = 3, alpha = 0.6,
          shape = 21,
          color = "steelblue",
          fill = "skyblue") +
  theme_minimal(base_size = 14) +
  labs(title = "Regression of LST on NDVI",
       x = "NDVI", y = "Land Surface Temperature (°C)")

# Define the helper function once
action_box <- function(title, items, border = "#e91e63", bg = "#fff5f7") {
  html <- paste0(
    '<div style="border-left:6px solid ', border,
    ';padding:12px 14px;background:', bg,
    ';border-radius:10px;margin:14px 0 6px 0">',
    '<h4 style="margin:0 0 8px 0;color:', border, ';font-weight:bold;">', title, '</h4>',
    '<ul style="margin:0;list-style-type:disc;color:black;">',
    paste0("<li>", items, "</li>", collapse = ""),
    "</ul></div>"
  )
  cat(html)
}
action_box(
  title  = "Action Items — Frequentist Models",
  items  = c(
    "What is the coefficient? The coefficient for NDVI is 2.465",
    "What is the p-value? p value for NDVI is 0.294", 
    "What does the CI mean? The CI is not shown", 
    "How would you describe the result if you were writing a results section? Do your best here. I mainly would like to see how you would present your results. We have not gone over this, so, feel free to use Google, etc for help. But, don’t blindly ask ChatGPT for it - I want to challenge you to start working on writing up results. We’ll come back to this as a class next week. A liner regression was used to look at whether NDVI predicts LST. The model estimated a positive, but not significant relationship"
  ),
  border = "#2e7d32", bg = "#eefaf0"
)

Action Items — Frequentist Models

# Load required packages
library(brms)      # for Bayesian regression
## Loading required package: Rcpp
## Loading 'brms' package (version 2.23.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:stats':
## 
##     ar
library(tidyverse) # for data handling
library(bayesplot) # for posterior checks
## This is bayesplot version 1.14.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':
## 
##     rhat
# Bayesian model: group ~ LST
bayes_group <- brm(
  LST ~ group,
  data   = dat,
  family = gaussian(),
  prior  = c(set_prior("normal(0, 5)", class = "b"),
             set_prior("normal(0, 10)", class = "Intercept"),
             set_prior("student_t(3, 0, 5)", class = "sigma")),
  iter   = 2000,
  chains = 4,
  warmup = 500,
  seed   = 123
)
## Compiling Stan program...
## WARNING: Rtools is required to build R packages, but is not currently installed.
## 
## Please download and install the appropriate version of Rtools for 4.5.1 from
## https://cran.r-project.org/bin/windows/Rtools/.
## Trying to compile a simple C file
## Running "C:/PROGRA~1/R/R-45~1.1/bin/x64/Rcmd.exe" SHLIB foo.c
## using C compiler: 'gcc.exe (GCC) 14.2.0'
## gcc  -I"C:/PROGRA~1/R/R-45~1.1/include" -DNDEBUG   -I"C:/Users/orkne/AppData/Local/R/win-library/4.5/Rcpp/include/"  -I"C:/Users/orkne/AppData/Local/R/win-library/4.5/RcppEigen/include/"  -I"C:/Users/orkne/AppData/Local/R/win-library/4.5/RcppEigen/include/unsupported"  -I"C:/Users/orkne/AppData/Local/R/win-library/4.5/BH/include" -I"C:/Users/orkne/AppData/Local/R/win-library/4.5/StanHeaders/include/src/"  -I"C:/Users/orkne/AppData/Local/R/win-library/4.5/StanHeaders/include/"  -I"C:/Users/orkne/AppData/Local/R/win-library/4.5/RcppParallel/include/" -DRCPP_PARALLEL_USE_TBB=1 -I"C:/Users/orkne/AppData/Local/R/win-library/4.5/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DUSE_STANC3 -DSTRICT_R_HEADERS  -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION  -D_HAS_AUTO_PTR_ETC=0  -include "C:/Users/orkne/AppData/Local/R/win-library/4.5/StanHeaders/include/stan/math/prim/fun/Eigen.hpp"  -std=c++1y    -I"C:/rtools45/x86_64-w64-mingw32.static.posix/include"      -O2 -Wall -std=gnu2x  -mfpmath=sse -msse2 -mstackrealign   -c foo.c -o foo.o
## cc1.exe: warning: command-line option '-std=c++14' is valid for C++/ObjC++ but not for C
## In file included from C:/Users/orkne/AppData/Local/R/win-library/4.5/RcppEigen/include/Eigen/Core:19,
##                  from C:/Users/orkne/AppData/Local/R/win-library/4.5/RcppEigen/include/Eigen/Dense:1,
##                  from C:/Users/orkne/AppData/Local/R/win-library/4.5/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22,
##                  from <command-line>:
## C:/Users/orkne/AppData/Local/R/win-library/4.5/RcppEigen/include/Eigen/src/Core/util/Macros.h:679:10: fatal error: cmath: No such file or directory
##   679 | #include <cmath>
##       |          ^~~~~~~
## compilation terminated.
## make: *** [C:/PROGRA~1/R/R-45~1.1/etc/x64/Makeconf:289: foo.o] Error 1
## WARNING: Rtools is required to build R packages, but is not currently installed.
## 
## Please download and install the appropriate version of Rtools for 4.5.1 from
## https://cran.r-project.org/bin/windows/Rtools/.
## Start sampling
## 
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 1).
## Chain 1: 
## Chain 1: Gradient evaluation took 3.2e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.32 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1: 
## Chain 1: 
## Chain 1: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 1: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 1: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 1: Iteration:  501 / 2000 [ 25%]  (Sampling)
## Chain 1: Iteration:  700 / 2000 [ 35%]  (Sampling)
## Chain 1: Iteration:  900 / 2000 [ 45%]  (Sampling)
## Chain 1: Iteration: 1100 / 2000 [ 55%]  (Sampling)
## Chain 1: Iteration: 1300 / 2000 [ 65%]  (Sampling)
## Chain 1: Iteration: 1500 / 2000 [ 75%]  (Sampling)
## Chain 1: Iteration: 1700 / 2000 [ 85%]  (Sampling)
## Chain 1: Iteration: 1900 / 2000 [ 95%]  (Sampling)
## Chain 1: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 1: 
## Chain 1:  Elapsed Time: 0.017 seconds (Warm-up)
## Chain 1:                0.02 seconds (Sampling)
## Chain 1:                0.037 seconds (Total)
## Chain 1: 
## 
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 2).
## Chain 2: 
## Chain 2: Gradient evaluation took 6e-06 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.06 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2: 
## Chain 2: 
## Chain 2: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 2: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 2: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 2: Iteration:  501 / 2000 [ 25%]  (Sampling)
## Chain 2: Iteration:  700 / 2000 [ 35%]  (Sampling)
## Chain 2: Iteration:  900 / 2000 [ 45%]  (Sampling)
## Chain 2: Iteration: 1100 / 2000 [ 55%]  (Sampling)
## Chain 2: Iteration: 1300 / 2000 [ 65%]  (Sampling)
## Chain 2: Iteration: 1500 / 2000 [ 75%]  (Sampling)
## Chain 2: Iteration: 1700 / 2000 [ 85%]  (Sampling)
## Chain 2: Iteration: 1900 / 2000 [ 95%]  (Sampling)
## Chain 2: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 2: 
## Chain 2:  Elapsed Time: 0.011 seconds (Warm-up)
## Chain 2:                0.022 seconds (Sampling)
## Chain 2:                0.033 seconds (Total)
## Chain 2: 
## 
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 3).
## Chain 3: 
## Chain 3: Gradient evaluation took 7e-06 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.07 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3: 
## Chain 3: 
## Chain 3: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 3: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 3: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 3: Iteration:  501 / 2000 [ 25%]  (Sampling)
## Chain 3: Iteration:  700 / 2000 [ 35%]  (Sampling)
## Chain 3: Iteration:  900 / 2000 [ 45%]  (Sampling)
## Chain 3: Iteration: 1100 / 2000 [ 55%]  (Sampling)
## Chain 3: Iteration: 1300 / 2000 [ 65%]  (Sampling)
## Chain 3: Iteration: 1500 / 2000 [ 75%]  (Sampling)
## Chain 3: Iteration: 1700 / 2000 [ 85%]  (Sampling)
## Chain 3: Iteration: 1900 / 2000 [ 95%]  (Sampling)
## Chain 3: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 3: 
## Chain 3:  Elapsed Time: 0.017 seconds (Warm-up)
## Chain 3:                0.022 seconds (Sampling)
## Chain 3:                0.039 seconds (Total)
## Chain 3: 
## 
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 4).
## Chain 4: 
## Chain 4: Gradient evaluation took 1.1e-05 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.11 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4: 
## Chain 4: 
## Chain 4: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 4: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 4: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 4: Iteration:  501 / 2000 [ 25%]  (Sampling)
## Chain 4: Iteration:  700 / 2000 [ 35%]  (Sampling)
## Chain 4: Iteration:  900 / 2000 [ 45%]  (Sampling)
## Chain 4: Iteration: 1100 / 2000 [ 55%]  (Sampling)
## Chain 4: Iteration: 1300 / 2000 [ 65%]  (Sampling)
## Chain 4: Iteration: 1500 / 2000 [ 75%]  (Sampling)
## Chain 4: Iteration: 1700 / 2000 [ 85%]  (Sampling)
## Chain 4: Iteration: 1900 / 2000 [ 95%]  (Sampling)
## Chain 4: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 4: 
## Chain 4:  Elapsed Time: 0.014 seconds (Warm-up)
## Chain 4:                0.02 seconds (Sampling)
## Chain 4:                0.034 seconds (Total)
## Chain 4:
summary(bayes_group)
##  Family: gaussian 
##   Links: mu = identity 
## Formula: LST ~ group 
##    Data: dat (Number of observations: 60) 
##   Draws: 4 chains, each with iter = 2000; warmup = 500; thin = 1;
##          total post-warmup draws = 6000
## 
## Regression Coefficients:
##            Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept     30.35      0.34    29.69    31.01 1.00     6217     4133
## groupUrban     4.02      0.47     3.10     4.95 1.00     5714     4047
## 
## Further Distributional Parameters:
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     1.86      0.18     1.55     2.25 1.00     5287     4325
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
bayes_R2(bayes_group)
##     Estimate  Est.Error     Q2.5     Q97.5
## R2 0.5488601 0.06017632 0.411371 0.6421825
pp_check(bayes_group, ndraws = 100)

# Bayesian model: NDVI ~ LST
bayes_reg <- brm(
  LST ~ NDVI,
  data   = dat,
  family = gaussian(),
  prior  = c(set_prior("normal(0, 5)", class = "b"),
             set_prior("normal(0, 10)", class = "Intercept"),
             set_prior("student_t(3, 0, 5)", class = "sigma")),
  iter   = 2000,
  chains = 4,
  warmup = 500,
  seed   = 123
)
## Compiling Stan program...
## WARNING: Rtools is required to build R packages, but is not currently installed.
## 
## Please download and install the appropriate version of Rtools for 4.5.1 from
## https://cran.r-project.org/bin/windows/Rtools/.
## Trying to compile a simple C file
## Running "C:/PROGRA~1/R/R-45~1.1/bin/x64/Rcmd.exe" SHLIB foo.c
## using C compiler: 'gcc.exe (GCC) 14.2.0'
## gcc  -I"C:/PROGRA~1/R/R-45~1.1/include" -DNDEBUG   -I"C:/Users/orkne/AppData/Local/R/win-library/4.5/Rcpp/include/"  -I"C:/Users/orkne/AppData/Local/R/win-library/4.5/RcppEigen/include/"  -I"C:/Users/orkne/AppData/Local/R/win-library/4.5/RcppEigen/include/unsupported"  -I"C:/Users/orkne/AppData/Local/R/win-library/4.5/BH/include" -I"C:/Users/orkne/AppData/Local/R/win-library/4.5/StanHeaders/include/src/"  -I"C:/Users/orkne/AppData/Local/R/win-library/4.5/StanHeaders/include/"  -I"C:/Users/orkne/AppData/Local/R/win-library/4.5/RcppParallel/include/" -DRCPP_PARALLEL_USE_TBB=1 -I"C:/Users/orkne/AppData/Local/R/win-library/4.5/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DUSE_STANC3 -DSTRICT_R_HEADERS  -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION  -D_HAS_AUTO_PTR_ETC=0  -include "C:/Users/orkne/AppData/Local/R/win-library/4.5/StanHeaders/include/stan/math/prim/fun/Eigen.hpp"  -std=c++1y    -I"C:/rtools45/x86_64-w64-mingw32.static.posix/include"      -O2 -Wall -std=gnu2x  -mfpmath=sse -msse2 -mstackrealign   -c foo.c -o foo.o
## cc1.exe: warning: command-line option '-std=c++14' is valid for C++/ObjC++ but not for C
## In file included from C:/Users/orkne/AppData/Local/R/win-library/4.5/RcppEigen/include/Eigen/Core:19,
##                  from C:/Users/orkne/AppData/Local/R/win-library/4.5/RcppEigen/include/Eigen/Dense:1,
##                  from C:/Users/orkne/AppData/Local/R/win-library/4.5/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22,
##                  from <command-line>:
## C:/Users/orkne/AppData/Local/R/win-library/4.5/RcppEigen/include/Eigen/src/Core/util/Macros.h:679:10: fatal error: cmath: No such file or directory
##   679 | #include <cmath>
##       |          ^~~~~~~
## compilation terminated.
## make: *** [C:/PROGRA~1/R/R-45~1.1/etc/x64/Makeconf:289: foo.o] Error 1
## WARNING: Rtools is required to build R packages, but is not currently installed.
## 
## Please download and install the appropriate version of Rtools for 4.5.1 from
## https://cran.r-project.org/bin/windows/Rtools/.
## Start sampling
## 
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 1).
## Chain 1: 
## Chain 1: Gradient evaluation took 7.5e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.75 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1: 
## Chain 1: 
## Chain 1: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 1: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 1: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 1: Iteration:  501 / 2000 [ 25%]  (Sampling)
## Chain 1: Iteration:  700 / 2000 [ 35%]  (Sampling)
## Chain 1: Iteration:  900 / 2000 [ 45%]  (Sampling)
## Chain 1: Iteration: 1100 / 2000 [ 55%]  (Sampling)
## Chain 1: Iteration: 1300 / 2000 [ 65%]  (Sampling)
## Chain 1: Iteration: 1500 / 2000 [ 75%]  (Sampling)
## Chain 1: Iteration: 1700 / 2000 [ 85%]  (Sampling)
## Chain 1: Iteration: 1900 / 2000 [ 95%]  (Sampling)
## Chain 1: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 1: 
## Chain 1:  Elapsed Time: 0.037 seconds (Warm-up)
## Chain 1:                0.058 seconds (Sampling)
## Chain 1:                0.095 seconds (Total)
## Chain 1: 
## 
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 2).
## Chain 2: 
## Chain 2: Gradient evaluation took 1.4e-05 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.14 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2: 
## Chain 2: 
## Chain 2: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 2: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 2: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 2: Iteration:  501 / 2000 [ 25%]  (Sampling)
## Chain 2: Iteration:  700 / 2000 [ 35%]  (Sampling)
## Chain 2: Iteration:  900 / 2000 [ 45%]  (Sampling)
## Chain 2: Iteration: 1100 / 2000 [ 55%]  (Sampling)
## Chain 2: Iteration: 1300 / 2000 [ 65%]  (Sampling)
## Chain 2: Iteration: 1500 / 2000 [ 75%]  (Sampling)
## Chain 2: Iteration: 1700 / 2000 [ 85%]  (Sampling)
## Chain 2: Iteration: 1900 / 2000 [ 95%]  (Sampling)
## Chain 2: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 2: 
## Chain 2:  Elapsed Time: 0.033 seconds (Warm-up)
## Chain 2:                0.05 seconds (Sampling)
## Chain 2:                0.083 seconds (Total)
## Chain 2: 
## 
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 3).
## Chain 3: 
## Chain 3: Gradient evaluation took 1.4e-05 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.14 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3: 
## Chain 3: 
## Chain 3: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 3: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 3: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 3: Iteration:  501 / 2000 [ 25%]  (Sampling)
## Chain 3: Iteration:  700 / 2000 [ 35%]  (Sampling)
## Chain 3: Iteration:  900 / 2000 [ 45%]  (Sampling)
## Chain 3: Iteration: 1100 / 2000 [ 55%]  (Sampling)
## Chain 3: Iteration: 1300 / 2000 [ 65%]  (Sampling)
## Chain 3: Iteration: 1500 / 2000 [ 75%]  (Sampling)
## Chain 3: Iteration: 1700 / 2000 [ 85%]  (Sampling)
## Chain 3: Iteration: 1900 / 2000 [ 95%]  (Sampling)
## Chain 3: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 3: 
## Chain 3:  Elapsed Time: 0.039 seconds (Warm-up)
## Chain 3:                0.055 seconds (Sampling)
## Chain 3:                0.094 seconds (Total)
## Chain 3: 
## 
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 4).
## Chain 4: 
## Chain 4: Gradient evaluation took 1.7e-05 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.17 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4: 
## Chain 4: 
## Chain 4: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 4: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 4: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 4: Iteration:  501 / 2000 [ 25%]  (Sampling)
## Chain 4: Iteration:  700 / 2000 [ 35%]  (Sampling)
## Chain 4: Iteration:  900 / 2000 [ 45%]  (Sampling)
## Chain 4: Iteration: 1100 / 2000 [ 55%]  (Sampling)
## Chain 4: Iteration: 1300 / 2000 [ 65%]  (Sampling)
## Chain 4: Iteration: 1500 / 2000 [ 75%]  (Sampling)
## Chain 4: Iteration: 1700 / 2000 [ 85%]  (Sampling)
## Chain 4: Iteration: 1900 / 2000 [ 95%]  (Sampling)
## Chain 4: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 4: 
## Chain 4:  Elapsed Time: 0.041 seconds (Warm-up)
## Chain 4:                0.067 seconds (Sampling)
## Chain 4:                0.108 seconds (Total)
## Chain 4:
summary(bayes_reg)
##  Family: gaussian 
##   Links: mu = identity 
## Formula: LST ~ NDVI 
##    Data: dat (Number of observations: 60) 
##   Draws: 4 chains, each with iter = 2000; warmup = 500; thin = 1;
##          total post-warmup draws = 6000
## 
## Regression Coefficients:
##           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept    31.35      1.11    29.17    33.56 1.00     5195     4115
## NDVI          2.03      2.14    -2.24     6.25 1.00     5365     3902
## 
## Further Distributional Parameters:
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     2.77      0.27     2.31     3.36 1.00     5265     4168
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
bayes_R2(bayes_reg)
##     Estimate  Est.Error         Q2.5     Q97.5
## R2 0.0254301 0.02912638 3.293696e-05 0.1066629
pp_check(bayes_reg, ndraws = 100)

# Define the helper function once
action_box <- function(title, items, border = "#e91e63", bg = "#fff5f7") {
  html <- paste0(
    '<div style="border-left:6px solid ', border,
    ';padding:12px 14px;background:', bg,
    ';border-radius:10px;margin:14px 0 6px 0">',
    '<h4 style="margin:0 0 8px 0;color:', border, ';font-weight:bold;">', title, '</h4>',
    '<ul style="margin:0;list-style-type:disc;color:black;">',
    paste0("<li>", items, "</li>", collapse = ""),
    "</ul></div>"
  )
  cat(html)
}
action_box(
  title  = "Action Items — Bayesian Models",
  items  = c(
    "Interpretation: What are the coefficient values? Intercept = 31.35, NDVI = 2.03",
    "What does the credible interval say? [-2.24, 6.25] crosses 0 so no clear effect",
    "How is this different from the lm() output? Positive estimate, but Bayesian focuses on uncertainty instead of p-values",
    "Are the models significant (in a Bayesian sense?). How do you tell? Not significant because CI includes 0",
    "How would you describe your results if you were writing them up? NDVI shows a weak, uncertain effect on LST"
  ),
  border = "#FF6EC7", bg = "#ffe6f5"
)
## <div style="border-left:6px solid #FF6EC7;padding:12px 14px;background:#ffe6f5;border-radius:10px;margin:14px 0 6px 0"><h4 style="margin:0 0 8px 0;color:#FF6EC7;font-weight:bold;">Action Items — Bayesian Models</h4><ul style="margin:0;list-style-type:disc;color:black;"><li>Interpretation: What are the coefficient values? Intercept = 31.35, NDVI = 2.03</li><li>What does the credible interval say? [-2.24, 6.25] crosses 0 so no clear effect</li><li>How is this different from the lm() output? Positive estimate, but Bayesian focuses on uncertainty instead of p-values</li><li>Are the models significant (in a Bayesian sense?). How do you tell? Not significant because CI includes 0</li><li>How would you describe your results if you were writing them up? NDVI shows a weak, uncertain effect on LST</li></ul></div>