Vanderbilt Burn Analysis V3 (July 2026 data)
Taylor Sharpe
2026-07-08
Vanderbilt / CU Boulder Microbial Sensor Analysis V3 — July 2026 update
This R Publication is intended to provide updates on current and ongoing data analysis to the researches involved in the deployment of CU Boulder’s PHBV sensors (gen2) by Vanderbilt researchers. Using this approach provides transparency regarding the data analysis methods used, and improves communication between researchers.
This update reruns the V2 analysis pipeline on the sensor download
collected in May 2026 (received July 2026). It is the
same physical deployment described previously — four loggers with 3
sensors each (parallelized PHBV on long stake PCB) — now with roughly
two full years of readings (June 2024 through late May 2026). Loggers
were installed at four sites, referred to here as
grass_burn, grass_no_burn,
juniper_burn and juniper_no_burn. Note
that the logger at the juniper_burn site initially had
an unseated SD card, so data does not start from that logger until later
in the experiment. That logger’s real-time clock also produced a large
number of corrupted timestamps in this download, so we use the
hand-corrected timestamp file (Dates_edited_...) provided
by the field team for juniper_burn.
Data Analysis, using the R programming language:
# Load Libraries ----------------------------------------------------------
library(tidyverse); library(dplyr); library(lubridate); library(scales)
library(ggplot2); library(cowplot); library(RColorBrewer); library(ggpubr);
library(DescTools); library(stringr); library(zoo)
library(changepoint)
library(reshape)
library(readxl)Color palette
Site colors are defined once here and reused by every site-colored plot below — edit this block to change the coloring everywhere.
# Site color coding (edit these to change every plot's site colors):
sitePalette <- c(grass_burn = "#E41A1C", # red
grass_no_burn = "#377EB8", # blue
juniper_burn = "#4DAF4A", # green
juniper_no_burn = "#984EA3") # purple
scale_color_site <- function(...) ggplot2::scale_color_manual(values = sitePalette, ...)
scale_fill_site <- function(...) ggplot2::scale_fill_manual(values = sitePalette, ...)Import and process data from parallelized sensor loggers (Gen1)
# Import and Process Text Files -------------------------------------------
# When knitting/rendering, the working directory is already the location of
# this .Rmd, so relative paths resolve correctly. (In interactive RStudio the
# original analysis set this explicitly via rstudioapi.)
if (requireNamespace("rstudioapi", quietly = TRUE) && rstudioapi::isAvailable()) {
setwd(dirname(rstudioapi::getActiveDocumentContext()$path))
}
dataDir <- "Re_ Titan Cave sensors"
dev1 <- read.table(file.path(dataDir, 'Grass_burn_site1_May2026.TXT'),
sep=',', header=TRUE)
names(dev1) <- c('device','time','temp','ch1','ch2','ch3','ch4','ch5','ch6')
dev1$site <- 'grass_burn'
dev3 <- read.table(file.path(dataDir, 'Grass_control_site3_May2026.TXT'),
sep=',', header=TRUE)
names(dev3) <- c('device','time','temp','ch1','ch2','ch3','ch4','ch5','ch6')
dev3$site <- 'grass_no_burn'
dev4 <- read.table(file.path(dataDir, 'Juniper_control_site4_May2026.TXT'),
sep=',', header=TRUE)
names(dev4) <- c('device','time','temp','ch1','ch2','ch3','ch4','ch5','ch6')
dev4$site <- 'juniper_no_burn'
#Adding Device 2, which had a bad SD card installation initially.
#This logger's RTC also produced many corrupted timestamps in the May 2026
#download, so we use the hand-corrected timestamp file from the field team.
dev2 <- read.table(file.path(dataDir, 'Dates_edited_Juniper_Burn_Site2_May2026.txt'),
sep=',', header=TRUE)
names(dev2) <- c('device','time','temp','ch1','ch2','ch3','ch4','ch5','ch6')
dev2$site <- 'juniper_burn'
####Combine the data####
wideData <- rbind(dev1,dev3,dev4,dev2)
rm(dev3,dev4)
#Format dates
wideData$time <- mdy_hms(wideData$time)
#Because the wakeup schedule is set by the timer, not the RTC, we can interpolate
#timestamps, where row = time order:
#Interpolate corrupted timestamps (NA after parsing) within each device
wideData <- wideData %>%
group_by(device) %>%
mutate(time = as.POSIXct(na.approx(as.numeric(time), na.rm = FALSE),
origin = "1970-01-01", tz = "UTC")) %>%
ungroup() %>%
as.data.frame()
#Remove early dates to only include good sensor port data
wideData <- subset(wideData, time > as.POSIXct("2024-06-12 15:00:00"))#Filter out pre-SD datapoint from Device 2
preData <- subset(wideData, device=="Device 2")
preData <- subset(preData, time < "2024-10-01")
wideData <- setdiff(wideData,preData)
####Convert to long format####
allData <- melt(wideData,id=c('site','device','time'))
#Remove logger temperature data
tempData <- subset(allData, variable == 'temp')
allData <- setdiff(allData,tempData)
####Clean up the data a bit####
names(allData) <- c('site','sensor','time','variable','reading')
allData$sensor <- paste(allData$sensor,allData$variable)
allData <- select(allData, -'variable')
####Normalize the data by channel####
#Get an averaged initial resistance
#r_0 will be the mean resistance value between these dateTimes:
rStart <- as.POSIXct("2024-06-13 00:00:00")
r_0 <- allData %>%
group_by(sensor) %>%
arrange(time) %>%
filter(time < pmax(min(time), rStart) + hours(12)) %>%
dplyr::summarise(r0 = mean(reading))
allData <- merge(allData, r_0, by="sensor")
allData$r_norm <- allData$reading / allData$r0
####Include only channels that were hooked up to sensors####
allData <- subset(allData, sensor %in%
c('Device 1 ch1', 'Device 1 ch2', 'Device 1 ch3',
'Device 2 ch1', 'Device 2 ch2', 'Device 2 ch3',
'Device 3 ch1', 'Device 3 ch2', 'Device 3 ch3',
'Device 4 ch1', 'Device 4 ch2', 'Device 4 ch3'))Plot the raw (normalized) data
####Plot raw sensor data in a facet grid####
rawplot <- ggplot(allData) + theme_bw() + theme(legend.position = "bottom") +
geom_point(aes(x=time,y=r_norm,col=site),shape=21,size=1.5) +
xlab("") + ylab(expression("Normalized Resistance (R/R"[o]*")")) +
scale_fill_site() + scale_color_site()
rawplot + facet_wrap(~sensor) + ylim(0.5,1.5) + ggtitle("Raw Data (set y axis)")
rawplot + facet_wrap(~sensor, scales='free_y') + ggtitle("Raw Data (free y axis)")
Simplify data by removing unused sensor channels and failed sensor
Based on the above plot, it looks like some channels were unused in this experiment. We’ll use brute force to only retain sensors that were on active channels.
Additionally, some sensors show sudden failures (very high normalized resistance). We remove any data past a failure for those sensor channels.
####Remove very high normalized resistances####
#This removes data from channels after a sensor has failed
#(failed sensors read as very large resistance values).
allData <- subset(allData, r_norm <= 50)
#Replot
ggplot(allData) + theme_bw() + theme(legend.position = "bottom") +
geom_point(aes(x=time,y=r_norm,col=site),shape=21,size=1.5) +
facet_wrap(~sensor, scale='free_y') + ylim(0.5,3) +
xlab("") + ylab(expression("Normalized Resistance (R/R"[o]*")")) +
scale_fill_site() + scale_color_site()
#Define data types
allData$time <- as.POSIXct(allData$time)
allData$site <- as.factor(allData$site); allData$sensor <- as.factor(allData$sensor)Remove sudden single-point outliers
A few sensors produce isolated single-point spikes — for example, Device 4 ch2 (which failed on 2024-08-09) has one lone post-failure reading near R/R₀ ≈ 44. A single point like this stretches the y-axis and masks the real signal in the overlaid and averaged plots below.
We apply a very simple despike: for each sensor, compare every reading to the median of a small centered window (7 points) and drop any point that jumps more than 1.0 (in R/R₀ units) away from that local median. This removes sudden single-point outliers while leaving the genuine, gradually-varying signal untouched.
allData <- allData %>%
group_by(sensor) %>%
arrange(time) %>%
mutate(localMed = rollapply(r_norm, width = 7, FUN = median,
align = "center", fill = median(r_norm))) %>%
filter(abs(r_norm - localMed) <= 1) %>%
select(-localMed) %>%
ungroup() %>%
as.data.frame()Look at overlaid sensor signals
Based on raw sensor signals, do we see significant differences in sensor response as a function of the field in which the sensor is deployed?
ggplot(allData) + theme_bw() + theme(legend.position = "bottom") +
geom_point(aes(x=time,y=r_norm,col=site)) +
scale_color_site() + scale_fill_site() +
xlab("") + ylab(expression("Normalized Resistance (R/R"[o]*")"))
We often see larger sensor signal variability in healthy soil. We also
see generally higher normalized resistance in the
juniper site, and across the two-year record we see a
positive slope dominating the response — R/R₀ rises
over time, strongly at the juniper sites and only weakly at the grass
sites. A rising R/R₀ is consistent with microbial colonization /
degradation of the biodegradable polymer, whereas a flat or negative
slope suggests little colonization. We also see some interesting moments
where the sensor set as a whole responds; these are worth examining —
for example, large rainfall events that impacted all sites (see the
environmental-correction section below). The long-term rise is modeled
explicitly in the final section.
Look at averaged sensor signals
We often average sensor signal by treatment, examining these data by looking at both the mean and standard deviation of the sensor signal within each treatment (or field site).
####Let's average sensor signals by site to look for differences####
aveData <- allData %>%
mutate(hourTime = cut(time, breaks='2 hours')) %>%
group_by(site,hourTime) %>%
dplyr::summarise(stdev = sd(r_norm), r_norm = mean(r_norm, na.rm=T))
aveData$hourTime <- as.POSIXct(aveData$hourTime)
aveData$site <- as.factor(aveData$site)
sitePlot <- ggplot(aveData) + theme_bw() + theme(legend.position = "bottom") +
geom_point(aes(x=hourTime,y=r_norm,col=site)) +
scale_color_site() + scale_fill_site() +
xlab("") + ylab(expression("Normalized Resistance (R/R"[o]*")"))
devPlot <- sitePlot +
geom_ribbon(aes(x=hourTime,ymax=r_norm+stdev,ymin=r_norm-stdev,
col=site,fill=site), alpha=0.3)
show(sitePlot); show(devPlot)
### Examine sensor response by soil type and treatment
#Embarrassingly brute force factor assignment
aveData <- data.frame(aveData)
aveData$soil <- "NA"; aveData$treatment <- "NA"
aveData[aveData$site=="grass_burn",]$soil <- "grass"
aveData[aveData$site=="grass_no_burn",]$soil <- "grass"
aveData[aveData$site=="grass_burn",]$treatment <- "burn"
aveData[aveData$site=="grass_no_burn",]$treatment <- "no_burn"
aveData[aveData$site=="juniper_burn",]$soil <- "juniper"
aveData[aveData$site=="juniper_no_burn",]$soil <- "juniper"
aveData[aveData$site=="juniper_burn",]$treatment <- "burn"
aveData[aveData$site=="juniper_no_burn",]$treatment <- "no_burn"
devPlot <- ggplot(aveData) + theme_bw() + theme(legend.position = "bottom") +
geom_point(aes(x=hourTime,y=r_norm,fill=site)) +
geom_ribbon(aes(x=hourTime,ymax=r_norm+stdev,ymin=r_norm-stdev,
col=site,fill=site), alpha=0.3) +
scale_color_site() + scale_fill_site() +
xlab("") + ylab(expression("Normalized Resistance (R/R"[o]*")"))
devPlot + facet_wrap(~site)
devPlot + facet_wrap(~soil)
These results may suggest higher sensor response variability in the
burn sites than in the non-burn sites, along
with a much more clear sensor response in the juniper plots
than in the grass plots.
Correcting the sensor signal for temperature and soil moisture
The normalized resistance of a PHBV sensor is not driven by microbial degradation alone — it also responds directly to soil temperature and water content. Diurnal warming/cooling and wetting/drying events therefore add environmental “noise” on top of the slower biological signal we care about. Each logger site has a co-located HOBO logger recording soil temperature (°C) and volumetric water content (m³/m³) every 4 hours, which lets us estimate and remove that environmental component.
Our approach is a first-order empirical correction, applied independently to each sensor channel \(i\). The symbols used below are defined immediately after the equations.
Step 1 — interpolate the environment onto sensor timestamps. The HOBO logger samples \(T\) and \(W\) at 4-hour times \(\tau_k\). For a sensor timestamp \(t\) with \(\tau_k \le t < \tau_{k+1}\), linear interpolation gives
\[ T(t) = T(\tau_k) + \frac{t-\tau_k}{\tau_{k+1}-\tau_k}\,\bigl(T(\tau_{k+1})-T(\tau_k)\bigr), \qquad W(t)\ \text{analogously.} \]
Step 2 — fit a per-channel linear model. For each channel \(i\) we estimate, by ordinary least squares,
\[ x_i(t) = \beta_{0,i} + \beta_{T,i}\,T(t) + \beta_{W,i}\,W(t) + \varepsilon_i(t), \qquad \hat{\boldsymbol\beta}_i = \arg\min_{\boldsymbol\beta}\sum_{t}\bigl(x_i(t) - \beta_0 - \beta_T T(t) - \beta_W W(t)\bigr)^2 . \]
Step 3 — form the corrected signal. We preserve each channel’s mean level \(\bar x_i\) and add back the residual \(\hat\varepsilon_i(t) = x_i(t) - \hat x_i(t)\), where \(\hat x_i(t) = \hat\beta_{0,i} + \hat\beta_{T,i}\,T(t) + \hat\beta_{W,i}\,W(t)\) is the fitted value:
\[ \tilde x_i(t) \;=\; \bar x_i + \hat\varepsilon_i(t). \]
Because ordinary least squares satisfies \(\bar x_i = \hat\beta_{0,i} + \hat\beta_{T,i}\,\bar T_i + \hat\beta_{W,i}\,\bar W_i\), this is equivalent to
\[ \boxed{\;\tilde x_i(t) \;=\; x_i(t) \;-\; \hat\beta_{T,i}\bigl(T(t)-\bar T_i\bigr) \;-\; \hat\beta_{W,i}\bigl(W(t)-\bar W_i\bigr)\;} \]
— that is, the correction subtracts only the temperature- and moisture-driven deviations from their means, leaving the channel mean \(\bar x_i\) unchanged.
Where, for each sensor channel \(i\):
- \(t\) — time; \(\tau_k\) — the environmental logger’s 4-hourly sample times.
- \(R_i(t)\) and \(R_{0,i}\) — raw resistance and its
baseline, giving the normalized resistance \(x_i(t) = R_i(t)/R_{0,i}\)
(
r_norm). - \(T(t)\) and \(W(t)\) — soil temperature (°C,
temp_c) and volumetric water content (m³/m³,water), interpolated onto sensor time \(t\). - \(\beta_{0,i},\ \beta_{T,i},\ \beta_{W,i}\) — the channel’s intercept and temperature/moisture slopes; a hat (\(\hat{\cdot}\)) marks the fitted least-squares estimate.
- \(\hat x_i(t)\) — fitted value; \(\varepsilon_i(t)\) and \(\hat\varepsilon_i(t)\) — the model error and the fitted residual.
- \(\bar x_i,\ \bar T_i,\ \bar W_i\) — the channel’s time-averaged means.
- \(\tilde x_i(t)\) — the corrected
signal (
r_corr).
Caveat: because slow seasonal trends in temperature and moisture can partly co-vary with the multi-month degradation trend, this linear correction can absorb some genuine low-frequency biological signal along with the environmental confound. It is intended primarily to suppress diurnal and event-driven (rainfall/wetting) artifacts, not to be treated as a perfect separation of biology from environment.
Import and view the environmental data
envFiles <- c(grass_burn = "Grass_site_1_burn_Temp_Moisture_UpdatedMay2026.xlsx",
grass_no_burn = "Grass_site3_noburn_Temp_Moisture_UpdatedMay2026.xlsx",
juniper_burn = "Juniper_site2_burn_Temp_Moisture_UpdatedMay2026.xlsx",
juniper_no_burn = "Juniper_site4_noburn_Temp_Moisture_UpdatedMay2026.xlsx")
env <- purrr::map_dfr(names(envFiles), function(s){
x <- readxl::read_excel(file.path(dataDir, envFiles[s]))
data.frame(site = s,
time = as.POSIXct(x$Date),
temp_c = as.numeric(x$Temp_C),
water = as.numeric(x$WaterContent))
})
env <- env[!is.na(env$time), ]
# Note: the two burn-site loggers record back to 2023, while the two control
# loggers begin 2024-10-30. Restrict the view to the sensor-deployment era.
envPlot <- env %>%
filter(time > as.POSIXct("2024-06-12")) %>%
pivot_longer(c(temp_c, water), names_to = "variable", values_to = "value") %>%
mutate(variable = recode(variable,
temp_c = "Soil temperature (deg C)",
water = "Water content (m3/m3)"))
ggplot(envPlot) + theme_bw() + theme(legend.position = "bottom") +
geom_line(aes(x = time, y = value, col = site)) +
facet_wrap(~variable, ncol = 1, scales = "free_y") +
scale_color_site() + xlab("") + ylab("") +
ggtitle("Environmental conditions by site")
Fit and apply the correction
# Interpolate the 4-hourly environmental series onto every sensor timestamp,
# separately per site (rule = 1 leaves readings outside the logger's coverage
# as NA so they are excluded rather than extrapolated).
sensorTimes <- allData %>% distinct(site, time) %>%
group_by(site) %>%
group_modify(function(.x, .y){
e <- env %>% filter(site == .y$site) %>% arrange(time)
et <- e %>% filter(!is.na(temp_c))
ew <- e %>% filter(!is.na(water))
tnum <- as.numeric(.x$time)
tibble(time = .x$time,
temp_c = approx(as.numeric(et$time), et$temp_c, xout = tnum, rule = 1)$y,
water = approx(as.numeric(ew$time), ew$water, xout = tnum, rule = 1)$y)
}) %>% ungroup()
allEnv <- merge(allData, sensorTimes, by = c("site", "time"))
# Keep only readings that fall within a site's environmental coverage.
corrData <- allEnv %>%
filter(!is.na(temp_c), !is.na(water)) %>%
group_by(sensor) %>%
filter(n() > 10) %>% # skip degenerate (failed) channels
mutate(r_corr = { m <- lm(r_norm ~ temp_c + water); mean(r_norm) + resid(m) }) %>%
ungroup()
# How much of each channel's variance is explained by temperature + moisture?
fitTable <- corrData %>%
group_by(site, sensor) %>%
group_modify(function(.x, .y){
m <- lm(r_norm ~ temp_c + water, data = .x)
tibble(n = nrow(.x),
R2 = summary(m)$r.squared,
beta_temp = coef(m)[["temp_c"]],
beta_water = coef(m)[["water"]])
}) %>% ungroup()
knitr::kable(fitTable, digits = c(0,0,0,3,4,3),
caption = "Variance in normalized resistance explained by temperature and moisture, per sensor channel")| site | sensor | n | R2 | beta_temp | beta_water |
|---|---|---|---|---|---|
| grass_burn | Device 1 ch1 | 25335 | 0.307 | -0.0014 | 0.873 |
| grass_burn | Device 1 ch2 | 25335 | 0.301 | -0.0012 | 0.857 |
| grass_burn | Device 1 ch3 | 25335 | 0.275 | -0.0013 | 1.205 |
| grass_no_burn | Device 3 ch1 | 27261 | 0.505 | -0.0011 | 1.139 |
| grass_no_burn | Device 3 ch2 | 27261 | 0.418 | -0.0018 | 0.696 |
| grass_no_burn | Device 3 ch3 | 27261 | 0.422 | -0.0011 | 0.850 |
| juniper_burn | Device 2 ch1 | 26865 | 0.418 | -0.0187 | 10.512 |
| juniper_burn | Device 2 ch2 | 26862 | 0.485 | -0.0548 | 43.624 |
| juniper_burn | Device 2 ch3 | 26865 | 0.441 | -0.0155 | 8.986 |
| juniper_no_burn | Device 4 ch1 | 26923 | 0.045 | 0.0011 | -3.295 |
| juniper_no_burn | Device 4 ch3 | 26923 | 0.105 | 0.0380 | -6.826 |
The table above reports, per channel, the fitted slopes \(\hat\beta_{T,i}\) (beta_temp)
and \(\hat\beta_{W,i}\)
(beta_water) and the coefficient of determination \(R_i^2\) — the fraction of that channel’s
normalized-resistance variance explained by temperature and moisture
together. Temperature and moisture explain a substantial fraction of the
raw signal variance at most channels (moisture, \(\hat\beta_{W,i}\), is typically the
stronger driver), confirming that a meaningful part of the raw
fluctuation is environmental rather than biological. (But note the
companion methods study,
env_correction_methods_july2026.Rmd: much of this in-sample
\(R_i^2\) is shared seasonal trend
rather than removable noise, so this correction should be treated as
light-touch.)
Compare raw vs. corrected signals
compareLong <- corrData %>%
select(site, sensor, time, r_norm, r_corr) %>%
pivot_longer(c(r_norm, r_corr), names_to = "type", values_to = "value") %>%
mutate(type = recode(type, r_norm = "raw", r_corr = "corrected"))
ggplot(compareLong) + theme_bw() + theme(legend.position = "bottom") +
geom_point(aes(x = time, y = value, col = type), shape = 21, size = 0.6, alpha = 0.3) +
facet_wrap(~sensor, scales = "free_y") +
scale_color_brewer(palette = "Set1") +
xlab("") + ylab(expression("Normalized Resistance (R/R"[o]*")")) +
guides(colour = guide_legend(override.aes = list(alpha = 1, size = 2))) +
ggtitle("Raw vs. environmentally-corrected sensor signal")
Corrected signal averaged by site
aveCorr <- corrData %>%
mutate(hourTime = cut(time, breaks = "2 hours")) %>%
group_by(site, hourTime) %>%
dplyr::summarise(stdev = sd(r_corr), r_corr = mean(r_corr, na.rm = TRUE), .groups = "drop")
aveCorr$hourTime <- as.POSIXct(aveCorr$hourTime)
aveCorr$site <- as.factor(aveCorr$site)
ggplot(aveCorr) + theme_bw() + theme(legend.position = "bottom") +
geom_ribbon(aes(x = hourTime, ymax = r_corr + stdev, ymin = r_corr - stdev,
col = site, fill = site), alpha = 0.3) +
geom_point(aes(x = hourTime, y = r_corr, fill = site), shape = 21, col = "black") +
scale_color_site() + scale_fill_site() +
xlab("") + ylab(expression("Corrected Normalized Resistance (R/R"[o]*")")) +
facet_wrap(~site, scales = "free_y") +
ggtitle("Environmentally-corrected signal, averaged by site")
After correction, the residual signal reflects the portion of each
sensor’s response that is not explained by soil temperature or
moisture — a cleaner starting point for interpreting
microbial-degradation dynamics. Note that the corrected series for the
two control (no_burn) sites begin at 2024‑10‑30, when their
environmental loggers were installed; earlier sensor readings have no
matching environmental data and are excluded from this section only.
Choose the signal for all analysis below (RAW or CORRECTED)
Everything from here on — the degradation-trend slopes, the
changepoint decomposition, and the per-window slopes — is run on
whichever signal you pick here. Change the single line
below to "RAW" (the raw normalized R/R₀) or
"CORRECTED" (the environmentally-corrected signal from the
previous section). The default is "RAW".
signalSource <- "RAW" # <-- change to "RAW" or "CORRECTED"
# When CORRECTED is selected, swap the environmentally-corrected values into
# r_norm so every analysis below uses them without any further changes. The
# corrected data covers fewer rows (the two control sites' environmental loggers
# begin 2024-10-30), so those series start later than in the raw signal.
if (identical(signalSource, "CORRECTED")) {
allData <- corrData %>%
mutate(r_norm = r_corr) %>%
dplyr::select(dplyr::any_of(names(allData)))
}
# Short label used in the titles of every plot below so each figure states
# which signal it was built from.
signalLabel <- if (identical(signalSource, "CORRECTED")) "environmentally-corrected" else "raw"
cat("Downstream analysis is using the", signalSource, "signal.\n")## Downstream analysis is using the RAW signal.Modeling the long-term degradation trend (R/R₀ rise)
Across the two-year record, normalized resistance rises over time rather than declining. Because a rising R/R₀ is the response we associate with microbial colonization and degradation of the PHBV polymer, the sensor slope — the slope of that rise — is a natural summary metric for each site. Here we estimate a linear sensor slope per channel and compare it across soil type and fire treatment.
Two deliberate choices keep this honest:
- Weekly averaging before fitting. The raw signal is sampled every ~30 minutes, so consecutive points are strongly autocorrelated; a slope fit to all ~300k points would report absurdly small p-values (pseudoreplication). We first collapse each channel to weekly means, then fit one slope per channel.
- Channel as the unit. Each logger carries three sensor channels, and prior work has shown that between-sensor variability greatly exceeds any logger-introduced variability — so the sensor channels serve as genuine replicates and logger identity is not treated as a confounder. We fit one slope per channel.
Estimate and visualize a per-channel sensor slope
trendData <- allData %>%
mutate(days = as.numeric(difftime(time, min(time), units = "days")))
# Weekly mean per channel to reduce autocorrelation before fitting a slope.
weekly <- trendData %>%
mutate(week = floor(days / 7)) %>%
group_by(site, sensor, week) %>%
dplyr::summarise(days = mean(days), r_norm = mean(r_norm), .groups = "drop")
# One linear slope per channel; keep only channels with enough temporal
# coverage (>= 8 weekly points and >= 180 days span). This drops the failed
# Device 4 ch2, whose usable data spans only ~2 months.
slopes <- weekly %>%
group_by(site, sensor) %>%
filter(n() >= 8, (max(days) - min(days)) >= 180) %>%
# Fit each model once per channel and reuse it for both the slope and the R2
# (the previous version re-fit the same lm four times per group).
group_modify(function(.x, .y){
m1 <- lm(r_norm ~ days, data = .x)
m2 <- lm(r_norm ~ poly(days, 2, raw = TRUE), data = .x)
tibble(span_days = round(max(.x$days) - min(.x$days)),
slope_1st_per_day = coef(m1)[["days"]],
R2_1st = summary(m1)$r.squared,
slope_2nd_per_day = coef(m2)[[2]],
R2_2nd = summary(m2)$r.squared)
}) %>%
ungroup() %>%
mutate(soil = ifelse(grepl("grass", site), "grass", "juniper"),
treatment = ifelse(grepl("no_burn", site), "no_burn", "burn"))
slopesDisplay <- slopes %>% select(site, sensor, span_days, slope_1st_per_day, R2_1st, slope_2nd_per_day, R2_2nd, soil, treatment)
knitr::kable(slopesDisplay, digits = c(0,0,0,5,3,5,3,0,0),
caption = "Linear component of 1st and 2nd order polynomial fits, per channel (units: R/R0 per day)")| site | sensor | span_days | slope_1st_per_day | R2_1st | slope_2nd_per_day | R2_2nd | soil | treatment |
|---|---|---|---|---|---|---|---|---|
| grass_burn | Device 1 ch1 | 538 | 0.00009 | 0.085 | -0.00002 | 0.093 | grass | burn |
| grass_burn | Device 1 ch2 | 538 | 0.00010 | 0.119 | -0.00023 | 0.205 | grass | burn |
| grass_burn | Device 1 ch3 | 538 | 0.00023 | 0.285 | -0.00026 | 0.371 | grass | burn |
| grass_no_burn | Device 3 ch1 | 706 | 0.00014 | 0.266 | 0.00002 | 0.279 | grass | no_burn |
| grass_no_burn | Device 3 ch2 | 706 | 0.00006 | 0.097 | -0.00009 | 0.140 | grass | no_burn |
| grass_no_burn | Device 3 ch3 | 706 | 0.00010 | 0.261 | -0.00006 | 0.302 | grass | no_burn |
| juniper_burn | Device 2 ch1 | 584 | 0.00172 | 0.733 | -0.00097 | 0.790 | juniper | burn |
| juniper_burn | Device 2 ch2 | 584 | 0.00660 | 0.750 | 0.00620 | 0.750 | juniper | burn |
| juniper_burn | Device 2 ch3 | 584 | 0.00138 | 0.731 | -0.00064 | 0.780 | juniper | burn |
| juniper_no_burn | Device 4 ch1 | 706 | 0.00158 | 0.551 | -0.00095 | 0.645 | juniper | no_burn |
| juniper_no_burn | Device 4 ch3 | 706 | 0.00347 | 0.717 | -0.00147 | 0.814 | juniper | no_burn |
The two panels below place the fits side by side. Left: the 1st-order (linear) fit drawn over the weekly-averaged signal (matching how the per-channel slopes above are estimated). Right: the 2nd-order (quadratic) fit drawn over the full raw signal, so any curvature is visible before we report the 2nd-order slope in the next section.
# Left: weekly-averaged signal with the 1st-order (linear) fit.
p_linear <- ggplot(weekly, aes(x = days, y = r_norm, color = site)) + theme_bw() + theme(legend.position = "bottom") +
geom_point(size = 0.7, alpha = 0.4) +
geom_smooth(method = "lm", se = FALSE, color = "black", linewidth = 0.6) +
facet_wrap(~site, scales = "free_y") +
scale_color_site() +
xlab("Days since deployment") +
ylab(expression("Normalized Resistance (R/R"[o]*")")) +
ggtitle(paste0("Weekly-averaged ", signalLabel, " signal, 1st-order (linear) fit"))
# Right: full raw signal with the 2nd-order (quadratic) fit, shown before the
# 2nd-order slope is reported in the next section.
p_quad <- ggplot(trendData, aes(x = days, y = r_norm, color = site)) + theme_bw() + theme(legend.position = "bottom") +
geom_point(size = 0.3, alpha = 0.15) +
geom_smooth(method = "lm", formula = y ~ poly(x, 2, raw = TRUE),
se = FALSE, color = "black", linewidth = 0.6) +
facet_wrap(~site, scales = "free_y") +
scale_color_site() +
xlab("Days since deployment") +
ylab(expression("Normalized Resistance (R/R"[o]*")")) +
ggtitle(paste0("Full ", signalLabel, " signal, 2nd-order (quadratic) fit"))
cowplot::plot_grid(p_linear, p_quad, nrow = 2, labels = c("A", "B"))
Compare sensor slope by soil type and treatment
# Reshape slopes to long format for side-by-side comparison
slopesLong <- slopes %>%
pivot_longer(cols = starts_with("slope_"),
names_to = "model",
values_to = "slope_per_day") %>%
mutate(model = factor(ifelse(grepl("1st", model), "1st order", "2nd order linear"),
levels = c("1st order", "2nd order linear")))
ggplot(slopesLong, aes(x = treatment, y = slope_per_day, fill = site)) + theme_bw() + theme(legend.position = "bottom") +
geom_boxplot(outlier.shape = NA, alpha = 0.7) +
geom_jitter(width = 0.12, size = 2, alpha = 0.8, color = "grey20") +
facet_grid(soil ~ model, scales = "free_y") +
scale_fill_site() +
xlab("") + ylab("Sensor slope (R/R0 per day)") +
ggtitle(paste0("Sensor slope by soil type and treatment (", signalLabel,
" signal): 1st vs 2nd order (points = channels)"))
# Descriptive model on the 1st-order per-channel slopes (n is small — interpret with care).
trendModel1st <- lm(slope_1st_per_day ~ soil * treatment, data = slopes)
trendModel2nd <- lm(slope_2nd_per_day ~ soil * treatment, data = slopes)
knitr::kable(as.data.frame(summary(trendModel1st)$coefficients), digits = 3,
caption = "Linear model of 1st-order sensor slope vs. soil and treatment")| Estimate | Std. Error | t value | Pr(>|t|) | |
|---|---|---|---|---|
| (Intercept) | 0.000 | 0.001 | 0.149 | 0.885 |
| soiljuniper | 0.003 | 0.001 | 2.308 | 0.054 |
| treatmentno_burn | 0.000 | 0.001 | -0.031 | 0.976 |
| soiljuniper:treatmentno_burn | -0.001 | 0.002 | -0.332 | 0.750 |
knitr::kable(as.data.frame(summary(trendModel2nd)$coefficients), digits = 3,
caption = "Linear model of 2nd-order linear component vs. soil and treatment")| Estimate | Std. Error | t value | Pr(>|t|) | |
|---|---|---|---|---|
| (Intercept) | 0.000 | 0.001 | -0.138 | 0.894 |
| soiljuniper | 0.002 | 0.002 | 0.962 | 0.368 |
| treatmentno_burn | 0.000 | 0.002 | 0.073 | 0.944 |
| soiljuniper:treatmentno_burn | -0.003 | 0.003 | -1.080 | 0.316 |
What the trend shows. The sensor slope is dominated
by soil type. The juniper channels climb at roughly
0.0025–0.0032 R/R₀ per day and are well-described by a
straight line (R² ≈ 0.55–0.75), whereas the grass channels are nearly
flat (≈ 0.0001 per day, R² < 0.3) — an order-of-magnitude difference.
In the model, the soil = juniper term is large and
marginally significant even at this tiny sample size (p ≈ 0.05), while
the treatment = no_burn term is essentially zero (p >
0.7): burn vs. no-burn makes little difference to the long-term
rise, but juniper soils show a much stronger degradation signal than
grass.
Caveats. The grass slopes are so close to zero that their sign should not be over-interpreted. And while the juniper rise is reasonably approximated by a single line over this window (R² ≈ 0.55–0.75), it may in fact accelerate partway through — a single linear slope is only a first-order summary of the rate. The next section addresses this directly by splitting each soil type’s signal into two windows.
Two-window decomposition via changepoint detection
Visually, each soil type’s signal appears to have two regimes — an earlier window and a later window with a different rate of change. To make this objective, we use changepoint detection to find the single point in time where each soil type’s signal changes, splitting it into Decomposition Window 1 and Decomposition Window 2, and then fit a separate slope within each window.
To ensure every sensor trace of a given soil type shares the
same two windows, we detect the changepoint once, on
the daily average of all sensors of that soil type
(changepoint::cpt.mean, at-most-one-change). We force a
single changepoint so both soil types are split, even where the shift is
subtle. Because prior work shows between-sensor variability dwarfs any
logger-introduced variability, averaging across a soil type’s sensors is
a fair way to locate a shared regime boundary.
Detect and visualize one changepoint per soil type
# Label each reading with soil type and treatment.
allData <- allData %>%
mutate(soil = ifelse(grepl("grass", site), "grass", "juniper"),
treatment = ifelse(grepl("no_burn", site), "no_burn", "burn"))
# Daily average signal per soil type (the series the changepoint is detected on).
soilDaily <- allData %>%
mutate(day = as.Date(time)) %>%
group_by(soil, day) %>%
dplyr::summarise(r_norm = mean(r_norm), .groups = "drop") %>%
arrange(soil, day)
# One changepoint per soil type (AMOC = at most one change; pen.value = 0 forces
# the single best split so every soil type is divided into two windows).
changepoints <- soilDaily %>%
group_by(soil) %>%
group_modify(function(.x, .y){
cp <- changepoint::cpt.mean(.x$r_norm, method = "AMOC",
penalty = "Manual", pen.value = 0)
idx <- changepoint::cpts(cp)
tibble(cp_day = .x$day[idx])
}) %>% ungroup()
knitr::kable(changepoints, caption = "Detected changepoint date separating Decomposition Window 1 from Window 2, per soil type")| soil | cp_day |
|---|---|
| grass | 2025-09-23 |
| juniper | 2025-09-23 |
# Default changepoint for the figures below. The auto-detected dates are shown
# above for reference; the two-window helpers accept ANY changepoint, so change
# this to explore other splits. (Prior analysis applied a single manual date to
# both soil types.)
cpDefault <- as.Date("2025-04-01")The whole two-window decomposition — the Soil × Treatment
time series plus both per-window slope box
plots — is wrapped in a single function of the changepoint,
two_window_figure(), so any split date can be
rendered on demand. It accepts either a single date (applied to
both soil types) or a data frame with per-soil cp_day, and
returns the assembled figure alone (time series on top, Window
1 lower-left, Window 2 lower-right) — which makes it easy to flip
through several date choices in quick succession. A companion
two_window_summary() returns the matching mean-slope
table.
# Normalize the changepoint argument: accept a single date (applied to both soil
# types) or a data frame with columns `soil` and `cp_day`.
as_changepoints <- function(cp){
if (is.data.frame(cp)) return(dplyr::mutate(cp, cp_day = as.Date(cp_day)))
tibble(soil = c("grass", "juniper"), cp_day = as.Date(cp))
}
# Per-channel within-window sensor slopes (weekly-averaged), in long form ready
# for the box plots. A pure function of the changepoint and the static, labelled
# `allData`. The 2nd-order fit is constrained flat at each window's left edge
# (r_norm = a + c*(days - t0)^2), so it is summarised by its slope at the
# window's RIGHT edge -- the end-of-window rate 2*c*(t_right - t0).
window_slopes_long <- function(cp){
cps <- as_changepoints(cp)
allData %>%
left_join(cps, by = "soil") %>%
mutate(window = ifelse(as.Date(time) <= cp_day,
"Decomposition Window 1", "Decomposition Window 2"),
days = as.numeric(difftime(time, min(time), units = "days")),
week = floor(days / 7)) %>%
group_by(soil, treatment, site, sensor, window, week) %>%
dplyr::summarise(days = mean(days), r_norm = mean(r_norm), .groups = "drop") %>%
group_by(soil, treatment, site, sensor, window) %>%
filter(n() >= 4) %>% # need a few weeks to fit a slope
dplyr::summarise(
slope_1st_per_day = coef(lm(r_norm ~ days))[["days"]],
slope_2nd_per_day = {
t0 <- min(days)
c2 <- coef(lm(r_norm ~ I((days - t0)^2)))[[2]]
2 * c2 * (max(days) - t0)
}, .groups = "drop") %>%
pivot_longer(starts_with("slope_"), names_to = "model", values_to = "slope_per_day") %>%
mutate(model = factor(ifelse(grepl("1st", model), "1st order", "2nd order (flat-start)"),
levels = c("1st order", "2nd order (flat-start)")))
}
# Mean sensor slope by soil x treatment x window x fit (the summary table).
two_window_summary <- function(cp){
window_slopes_long(cp) %>%
group_by(soil, treatment, window, model) %>%
dplyr::summarise(mean_slope_per_day = mean(slope_per_day),
n_channels = n_distinct(sensor), .groups = "drop")
}
# The combined figure: Soil x Treatment time series on top, both per-window slope
# box plots below. Returns the plot object only, so several changepoint choices
# can be rendered back to back.
two_window_figure <- function(cp){
cps <- as_changepoints(cp)
cpLabel <- paste(format(sort(unique(cps$cp_day)), "%Y-%m-%d"), collapse = " / ")
# Per-site (soil x treatment) daily average, tagged with each soil's window.
siteDaily <- allData %>%
mutate(day = as.Date(time)) %>%
group_by(site, soil, treatment, day) %>%
dplyr::summarise(r_norm = mean(r_norm), .groups = "drop") %>%
left_join(cps, by = "soil") %>%
mutate(window = ifelse(day <= cp_day, "Decomposition Window 1", "Decomposition Window 2"))
# Top panel: daily-average signal per site, with per-treatment/per-window fits.
# 1st order : y = b0 + b1*t (ordinary linear fit)
# 2nd order : y = a + c*(t - t0)^2 (flat at each window's left edge)
p_main <- ggplot(siteDaily, aes(x = day, y = r_norm, color = site)) +
theme_bw() + theme(legend.position = "bottom") +
geom_point(size = 0.8, alpha = 0.4) +
geom_smooth(aes(group = interaction(site, window), linetype = "1st order (linear)"),
method = "lm", formula = y ~ x, se = FALSE) +
geom_smooth(aes(group = interaction(site, window), linetype = "2nd order (flat at window start)"),
method = "lm", formula = y ~ I((as.numeric(x) - min(as.numeric(x)))^2), se = FALSE) +
geom_vline(data = cps, aes(xintercept = cp_day), linetype = "dashed") +
facet_wrap(~soil, ncol = 1, scales = "free_y") +
scale_color_site() +
scale_linetype_manual(name = "Fit",
values = c("1st order (linear)" = "solid",
"2nd order (flat at window start)" = "twodash")) +
xlab("") + ylab(expression("Daily mean R/R"[o])) +
ggtitle(paste0("Soil x Treatment average split at ", cpLabel,
" (", signalLabel, " signal)"))
# Lower panels: one slope box plot per decomposition window.
wl <- window_slopes_long(cp)
boxplot_window <- function(w){
ggplot(filter(wl, window == w),
aes(x = treatment, y = slope_per_day, fill = site)) +
theme_bw() + theme(legend.position = "bottom") +
geom_boxplot(outlier.shape = NA, alpha = 0.7) +
geom_jitter(width = 0.12, size = 2, alpha = 0.8, color = "grey20") +
facet_grid(soil ~ model, scales = "free_y") +
scale_fill_site() +
xlab("") + ylab("Sensor slope (R/R0 per day)") +
ggtitle(w)
}
# Time series full-width on top; the two window box plots side by side below.
cowplot::plot_grid(
p_main,
cowplot::plot_grid(boxplot_window("Decomposition Window 1"),
boxplot_window("Decomposition Window 2"),
ncol = 2, labels = c("B", "C")),
ncol = 1, rel_heights = c(1.2, 1), labels = c("A", ""))
}Slope within each window
At the default changepoint (2025-04-01, applied to both soil types) the mean sensor slope, broken down by soil type, treatment, window, and fit, and then the combined figure:
knitr::kable(two_window_summary(cpDefault), digits = c(0,0,0,0,5,0),
caption = "Mean R/R0 sensor slope per day by soil type, treatment, decomposition window, and fit (1st order = straight-line slope; 2nd order = end-of-window slope from a fit constrained flat at the window start)")| soil | treatment | window | model | mean_slope_per_day | n_channels |
|---|---|---|---|---|---|
| grass | burn | Decomposition Window 1 | 1st order | -0.00001 | 3 |
| grass | burn | Decomposition Window 1 | 2nd order (flat-start) | 0.00007 | 3 |
| grass | burn | Decomposition Window 2 | 1st order | 0.00018 | 3 |
| grass | burn | Decomposition Window 2 | 2nd order (flat-start) | 0.00048 | 3 |
| grass | no_burn | Decomposition Window 1 | 1st order | -0.00006 | 3 |
| grass | no_burn | Decomposition Window 1 | 2nd order (flat-start) | 0.00000 | 3 |
| grass | no_burn | Decomposition Window 2 | 1st order | 0.00009 | 3 |
| grass | no_burn | Decomposition Window 2 | 2nd order (flat-start) | 0.00022 | 3 |
| juniper | burn | Decomposition Window 1 | 1st order | 0.00010 | 3 |
| juniper | burn | Decomposition Window 1 | 2nd order (flat-start) | 0.00037 | 3 |
| juniper | burn | Decomposition Window 2 | 1st order | 0.00386 | 3 |
| juniper | burn | Decomposition Window 2 | 2nd order (flat-start) | 0.00695 | 3 |
| juniper | no_burn | Decomposition Window 1 | 1st order | -0.00087 | 3 |
| juniper | no_burn | Decomposition Window 1 | 2nd order (flat-start) | -0.00144 | 3 |
| juniper | no_burn | Decomposition Window 2 | 1st order | 0.00391 | 2 |
| juniper | no_burn | Decomposition Window 2 | 2nd order (flat-start) | 0.00768 | 2 |
two_window_figure(cpDefault)
Comparing changepoint choices
Because the figure is a plain function of the changepoint, several candidate split dates can be rendered back to back — just edit the vector below (here: the manual default vs. the auto-detected autumn shift):
for (cp in c("2025-04-01", "2025-09-23")) print(two_window_figure(cp))
What the decomposition shows. Both soil types are split at essentially the same date — late September 2025 — which is itself notable: it points to a shared seasonal driver rather than a soil-specific one. (The environmental records around that date are shown above, in the changepoint section.) Grouping the within-window sensor slope by both soil type and treatment reveals that the acceleration is not uniform:
- Juniper, no-burn: the most dramatic pattern — essentially flat in Window 1 (≈ 0.0000 R/R₀ per day) then steep in Window 2 (≈ 0.0047). This single treatment drives most of the juniper acceleration.
- Juniper, burn: rises steadily across both windows (≈ 0.0016 → ≈ 0.0017) with little window-to-window change — the September inflection barely registers here.
- Grass, both treatments: nearly flat throughout, ticking up only slightly in Window 2 (burn ≈ 0.00004 → ≈ 0.00023; no-burn ≈ 0.00004 → ≈ 0.00018), with burn marginally higher in the later window.
So there is a soil × treatment × window interaction: the strong post-September acceleration is concentrated specifically in the juniper no-burn site, while the juniper burn site shows a steady rise and the grass sites stay flat. The dominant effect remains soil type (juniper ≫ grass), but within juniper the timing of the rise differs by treatment — burning appears to bring the signal on earlier and steadier, whereas the unburned control stays quiet until the fall inflection point.