Mixed_Linear_Regression
2025-07-23
# importing libraries
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.1
## ✔ ggplot2 3.5.2 ✔ 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 errorslibrary(randomForest)## randomForest 4.7-1.2
## Type rfNews() to see new features/changes/bug fixes.
##
## Attaching package: 'randomForest'
##
## The following object is masked from 'package:dplyr':
##
## combine
##
## The following object is masked from 'package:ggplot2':
##
## marginlibrary(dplyr)
library(nlme)##
## Attaching package: 'nlme'
##
## The following object is masked from 'package:dplyr':
##
## collapselibrary(corrplot)## corrplot 0.95 loadedlibrary(reshape2)##
## Attaching package: 'reshape2'
##
## The following object is masked from 'package:tidyr':
##
## smithsThere are still 25 variables - this model can be simplified. To understand the impact of each variable a correlation matrix can show us which variables are highly correlated.
# before choosing variables for the LME lets look at correlation matrix
# this is a mess, and there are too many variables to properly interpret the results
# looking at only numerical variables
numeric_vars <- df_categorized[sapply(df_categorized, is.numeric)]
# creating corr matrix
cor_matrix <- cor(numeric_vars, use = "pairwise.complete.obs", method = "pearson")
#plotting matrix
corrplot(
cor_matrix,
method = "color", # Color-coded cells
type = "upper", # Only upper triangle
col = colorRampPalette(c("blue", "white", "darkgreen"))(200), # Custom colors
tl.cex = 0.7, # Smaller font size for variable names
tl.col = "black", # Text label color
tl.srt = 45, # Rotate variable labels
addCoef.col = "black", # Add correlation coefficients
number.cex = 0.5, # Size of coefficient text
mar = c(0, 0, 1, 0) # Reduce margin space
)
# extracting highest cor
threshold <- 0.7 #.7 is a personal choice
# Set the diagonal to 0 to ignore self-correlation
cor_matrix_no_diag <- cor_matrix
diag(cor_matrix_no_diag) <- 0
# Get variables that have at least one |cor| > threshold
keep_vars <- names(which(apply(abs(cor_matrix_no_diag), 1, function(x) any(x > threshold))))
cor_subset <- cor_matrix[keep_vars, keep_vars]
cor_long <- melt(cor_subset)
# pairwise list of highest corr in data set
high_corr_pairs <- subset(cor_long, abs(value) > threshold & Var1 != Var2)
print(high_corr_pairs)## Var1 Var2 value
## 5 VWCB DOY 0.7255849
## 8 BulkDensity DOY 0.8559357
## 9 TotalTurns DOY 0.8771654
## 10 un_scaled_DOY DOY 1.0000000
## 13 u* wind_speed 0.8161463
## 22 wind_speed u* 0.8161463
## 36 TMP TMPA 0.7902757
## 41 DOY VWCB 0.7255849
## 50 un_scaled_DOY VWCB 0.7255849
## 54 TMPA TMP 0.7902757
## 57 DEW TMP 0.8461523
## 66 TMP DEW 0.8461523
## 71 DOY BulkDensity 0.8559357
## 80 un_scaled_DOY BulkDensity 0.8559357
## 81 DOY TotalTurns 0.8771654
## 90 un_scaled_DOY TotalTurns 0.8771654
## 91 DOY un_scaled_DOY 1.0000000
## 95 VWCB un_scaled_DOY 0.7255849
## 98 BulkDensity un_scaled_DOY 0.8559357
## 99 TotalTurns un_scaled_DOY 0.8771654The variables which exhibit high levels of autocorrelation are not surprising. DOY is highly correlated with many variables due to the nature of sampling - large ranges of DOY share a single value as Bulk density, total turns etc. A surprising variable was the correlation between DOY and soil moisture - yet this may be due to seasonal patterns as soil drys out as summer progresses leading to lower average values as the summer progresses.
# a basic linear model gives us an idea of the relationship - we are looking here are scaled data - furthermore it is important to not despite the data being scaled it does not mean that one variable is more important than the other as the ranges are not the same.
full_lm <- lm(FCO2_DRY.LIN ~ ., data = df_scaled)
summary(full_lm)##
## Call:
## lm(formula = FCO2_DRY.LIN ~ ., data = df_scaled)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.1155 -0.2029 -0.0154 0.1635 4.3399
##
## Coefficients: (2 not defined because of singularities)
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.229918 1.317923 -0.933 0.350910
## PileE -1.430634 0.207821 -6.884 9.81e-12 ***
## Chamber_CorrectedC2 0.199413 0.102985 1.936 0.053084 .
## Chamber_CorrectedC3 -0.334356 0.060542 -5.523 4.17e-08 ***
## Chamber_CorrectedC4 -0.077030 0.118723 -0.649 0.516592
## Chamber_CorrectedC6 NA NA NA NA
## Chamber_Elevation 0.128374 0.024177 5.310 1.33e-07 ***
## DOY 0.009138 0.005527 1.653 0.098528 .
## air_temperature 0.016959 0.055139 0.308 0.758468
## RH -0.070881 0.049714 -1.426 0.154217
## Tdew 0.041239 0.049929 0.826 0.409010
## wind_speed 0.053048 0.027567 1.924 0.054578 .
## wind_dir 0.019936 0.013576 1.468 0.142260
## `u*` -0.060483 0.031346 -1.930 0.053923 .
## TMPA 0.385604 0.099916 3.859 0.000120 ***
## VWCB 0.067967 0.027179 2.501 0.012540 *
## PAR 0.008216 0.035263 0.233 0.815820
## BAR -0.045062 0.039156 -1.151 0.250051
## HMD -0.204266 0.189089 -1.080 0.280265
## TMP -0.454602 0.397979 -1.142 0.253591
## DEW 0.285376 0.303402 0.941 0.347125
## TA.range 0.197157 0.041886 4.707 2.84e-06 ***
## SWC_1.initial_value 0.008167 0.044941 0.182 0.855825
## TS_1.initial_value 0.075394 0.034596 2.179 0.029526 *
## FH2O -0.012787 0.018374 -0.696 0.486624
## BulkDensity 0.140485 0.039977 3.514 0.000459 ***
## DaysSinceLastTurn -0.159175 0.027401 -5.809 8.23e-09 ***
## TotalTurns -0.773780 0.141047 -5.486 5.11e-08 ***
## un_scaled_DOY NA NA NA NA
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4389 on 1090 degrees of freedom
## (8162 observations deleted due to missingness)
## Multiple R-squared: 0.5279, Adjusted R-squared: 0.5166
## F-statistic: 46.88 on 26 and 1090 DF, p-value: < 2.2e-16# Key takeaways - Days since turning is important, more important than total turns,
# bulk density is important, chamber elevation is important. Stevens probe measurements and watchdog data will be part of the model.
# odd that wind direction matters, is that because of the wind or is that correlated with something?
# using random forest as well:
# had to change u* name to avoid issues
# TA.range - which is chamber tempurate was an important varaible. Not sure if this is correlation or causation - not sure about this varaible. While it improves the model - im not sure what it represents.
df_categorized <- df_categorized %>%
rename(u_star = `u*`)
rf_model <- randomForest(FCO2_DRY.LIN ~ ., data = df_categorized, importance = TRUE, na.action = na.omit)
# Example 1: Increase margins and text size
par(mar = c(5, 8, 4, 2)) # Increase left margin for long labels
varImpPlot(rf_model, cex = 0.8) # cex = character size
In our Linear Model with 26 variables we can explain around 50% of the variability in magnitude of flux. Based on the basic linear model variables which are significant (not necessarily important just significant this is not looking at effect size). Pile, Chamber 3 (not sure about what is going on here), Chamber 2 to a lesser extent, Chamber elevation, windspeed, Temp Range (TA range which is temp range measured in chamber), pile temperature, bulk density, days since last turn, total turns.
Some of these relationships can be explored more - for example what is going on in chamber 3 and to a lesser extent 1. what is the relationship between airtemp range in chamber and emissions - I’m assuming it is positive. Is this a sign of advection?
Looking at the random forest model there are similiar trends.Looking at variables with the highest %IncMSE - Chamber - Days since last turn - SWC intial value - TS initial value - Bulk Density - TA range - Chamber elevation - meterological var - total turns is further down - interesting - though it still is an importnat variable
Both linear and RF models help us select important variables for Mixed effects Model.
# To help with model building removing outliers (being generous with a 2 times IQR could use a higher value as there are likely spikes from turning)
# defining the quantiles
Q1 <- quantile(df_categorized$FCO2_DRY.LIN, 0.25)
Q3 <- quantile(df_categorized$FCO2_DRY.LIN, 0.75)
# the interquartile range
iqr_val <- IQR(df_categorized$FCO2_DRY.LIN, na.rm = TRUE)
# creating the boundaries
lower_bound <- Q1 - 2* iqr_val
upper_bound <- Q3 + 2* iqr_val
# Filtering out outliers -
Scaled_Cleaned_Co2_data <- df_categorized[df_categorized$FCO2_DRY.LIN >= lower_bound & df_categorized$FCO2_DRY.LIN <= upper_bound, ]
# creating the base model.
# Variable choice:
# Pile
# Initial Soil Water Content from Steven's Probe
# Initial Temperature value from the Steven's Probe
# TMP - outside temp
# BAR - pressure
# Bulk Density
# Chamber Elevation
# Total turns
# HMD - Humidity
# TA.range
# Base model
gls.co2 = gls(FCO2_DRY.LIN ~ Pile, na.action = na.omit, data = Scaled_Cleaned_Co2_data)
# Base linear mixed effects model, Random effect of chamber
lme_base_model <- lme(
FCO2_DRY.LIN ~ Pile + SWC_1.initial_value + TS_1.initial_value +
TMP + BAR + BulkDensity + Chamber_Elevation +
DaysSinceLastTurn + TotalTurns + HMD + TA.range,
random = ~ 1 | Chamber_Corrected,
data = Scaled_Cleaned_Co2_data,
na.action = na.omit
)
summary(lme_base_model)## Linear mixed-effects model fit by REML
## Data: Scaled_Cleaned_Co2_data
## AIC BIC logLik
## 10908.88 11008.13 -5440.44
##
## Random effects:
## Formula: ~1 | Chamber_Corrected
## (Intercept) Residual
## StdDev: 0.2550088 0.4439518
##
## Fixed effects: FCO2_DRY.LIN ~ Pile + SWC_1.initial_value + TS_1.initial_value + TMP + BAR + BulkDensity + Chamber_Elevation + DaysSinceLastTurn + TotalTurns + HMD + TA.range
## Value Std.Error DF t-value p-value
## (Intercept) -0.13263637 0.14746868 8856 -0.899421 0.3685
## PileE -0.01636875 0.20876679 4 -0.078407 0.9413
## SWC_1.initial_value 0.09674275 0.00778367 8856 12.428936 0.0000
## TS_1.initial_value 0.20690055 0.00752481 8856 27.495794 0.0000
## TMP 0.04876489 0.00715835 8856 6.812309 0.0000
## BAR -0.00803582 0.00465582 8856 -1.725972 0.0844
## BulkDensity -0.04217384 0.00811760 8856 -5.195358 0.0000
## Chamber_Elevation 0.06878902 0.00921789 8856 7.462554 0.0000
## DaysSinceLastTurn -0.04432112 0.00581447 8856 -7.622555 0.0000
## TotalTurns 0.01381389 0.00942417 8856 1.465794 0.1427
## HMD 0.01720184 0.00648658 8856 2.651911 0.0080
## TA.range 0.08467662 0.00632583 8856 13.385849 0.0000
## Correlation:
## (Intr) PileE SWC_1. TS_1._ TMP BAR BlkDns Chmb_E
## PileE -0.707
## SWC_1.initial_value 0.017 -0.028
## TS_1.initial_value -0.009 0.016 -0.172
## TMP -0.001 0.001 0.108 -0.284
## BAR 0.001 -0.001 0.012 0.014 0.007
## BulkDensity 0.028 -0.043 0.107 0.125 -0.034 0.004
## Chamber_Elevation 0.008 -0.011 0.127 0.169 0.023 0.057 0.230
## DaysSinceLastTurn 0.005 -0.008 0.205 0.018 0.157 0.049 -0.048 0.192
## TotalTurns -0.024 0.037 0.154 0.153 0.267 0.039 -0.643 0.019
## HMD -0.001 0.001 0.099 -0.140 0.654 0.071 -0.019 -0.013
## TA.range 0.009 -0.010 -0.159 -0.288 0.118 0.003 0.027 0.035
## DysSLT TtlTrn HMD
## PileE
## SWC_1.initial_value
## TS_1.initial_value
## TMP
## BAR
## BulkDensity
## Chamber_Elevation
## DaysSinceLastTurn
## TotalTurns 0.390
## HMD 0.105 0.145
## TA.range 0.162 -0.009 -0.062
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -3.14483411 -0.67713831 -0.07315114 0.50652321 4.47632772
##
## Number of Observations: 8872
## Number of Groups: 6plot(lme_base_model)
qqnorm(lme_base_model)
# Base linear mixed effects model, Random effect of chamber pile
lme_base_model_pile <- lme(
FCO2_DRY.LIN ~ Pile + SWC_1.initial_value + TS_1.initial_value +
TMP + BAR + BulkDensity + Chamber_Elevation +
DaysSinceLastTurn + TotalTurns + HMD + TA.range,
random = ~ 1 | Pile,
data = Scaled_Cleaned_Co2_data,
na.action = na.omit
)
summary(lme_base_model_pile)## Warning in pt(-abs(tVal), fDF): NaNs produced## Linear mixed-effects model fit by REML
## Data: Scaled_Cleaned_Co2_data
## AIC BIC logLik
## 12254.53 12353.78 -6113.266
##
## Random effects:
## Formula: ~1 | Pile
## (Intercept) Residual
## StdDev: 0.01920313 0.479622
##
## Fixed effects: FCO2_DRY.LIN ~ Pile + SWC_1.initial_value + TS_1.initial_value + TMP + BAR + BulkDensity + Chamber_Elevation + DaysSinceLastTurn + TotalTurns + HMD + TA.range
## Value Std.Error DF t-value p-value
## (Intercept) -0.08493835 0.021181635 8860 -4.010000 0.0001
## PileE -0.12507642 0.031440308 0 -3.978219 NaN
## SWC_1.initial_value 0.21407555 0.006872031 8860 31.151716 0.0000
## TS_1.initial_value 0.15095753 0.006817083 8860 22.144006 0.0000
## TMP 0.07919202 0.007554753 8860 10.482410 0.0000
## BAR -0.00596967 0.005028484 8860 -1.187170 0.2352
## BulkDensity -0.00327541 0.008601165 8860 -0.380810 0.7034
## Chamber_Elevation 0.10974068 0.009679093 8860 11.337909 0.0000
## DaysSinceLastTurn -0.00506534 0.006096565 8860 -0.830852 0.4061
## TotalTurns 0.02747159 0.010055983 8860 2.731866 0.0063
## HMD 0.03779601 0.006952172 8860 5.436576 0.0000
## TA.range 0.11160677 0.006701387 8860 16.654280 0.0000
## Correlation:
## (Intr) PileE SWC_1. TS_1._ TMP BAR BlkDns Chmb_E
## PileE -0.728
## SWC_1.initial_value 0.112 -0.163
## TS_1.initial_value -0.048 0.063 -0.006
## TMP -0.020 0.036 0.030 -0.207
## BAR 0.003 -0.004 0.000 0.025 0.004
## BulkDensity 0.212 -0.311 0.151 0.133 -0.033 0.005
## Chamber_Elevation 0.060 -0.079 0.197 0.222 0.019 0.059 0.200
## DaysSinceLastTurn 0.022 -0.035 0.125 0.078 0.140 0.045 -0.059 0.188
## TotalTurns -0.188 0.276 0.114 0.166 0.277 0.037 -0.653 0.040
## HMD -0.016 0.024 0.047 -0.094 0.648 0.070 -0.019 -0.015
## TA.range 0.058 -0.061 -0.226 -0.272 0.097 0.001 0.008 -0.002
## DysSLT TtlTrn HMD
## PileE
## SWC_1.initial_value
## TS_1.initial_value
## TMP
## BAR
## BulkDensity
## Chamber_Elevation
## DaysSinceLastTurn
## TotalTurns 0.385
## HMD 0.090 0.144
## TA.range 0.141 0.002 -0.077
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -2.8057142 -0.6070480 -0.1487151 0.4575224 4.2866869
##
## Number of Observations: 8872
## Number of Groups: 2plot(lme_base_model_pile)
qqnorm(lme_base_model_pile)
# Create a dataframe of model AIC and BIC
model_stats <- tibble(
Model = c("No Random Effects", "Chamber", "Pile"),
AIC = c(AIC(gls.co2), AIC(lme_base_model), AIC(lme_base_model_pile)),
BIC = c(BIC(gls.co2), BIC(lme_base_model), BIC(lme_base_model_pile))
)
print(model_stats)## # A tibble: 3 × 3
## Model AIC BIC
## <chr> <dbl> <dbl>
## 1 No Random Effects 15322. 15343.
## 2 Chamber 10909. 11008.
## 3 Pile 12255. 12354.# When chamber is the random effect BAR and Bulk density add to the model, yet when pile is used as a random effect (which I do not think makes sense) they stop being significant. Why is that ? The best model is the Chamber with the lowest AIC and BIC values followed by pile which also significantly improves model. Here a crucial question what does the random effect of chamber tell us. Do we expect each chamber to behave differently or this just masking the impact of management?
lme_base_model_nested <- lme(
FCO2_DRY.LIN ~ Pile + SWC_1.initial_value + TS_1.initial_value +
TMP + BAR + BulkDensity + Chamber_Elevation +
DaysSinceLastTurn + TotalTurns + HMD + TA.range,
random = ~ 1 | Pile/Chamber_Corrected,
data = Scaled_Cleaned_Co2_data,
na.action = na.omit
)
anova(lme_base_model_pile,lme_base_model_nested)## Model df AIC BIC logLik Test L.Ratio
## lme_base_model_pile 1 14 12254.53 12353.78 -6113.266
## lme_base_model_nested 2 15 10910.88 11017.22 -5440.440 1 vs 2 1345.652
## p-value
## lme_base_model_pile
## lme_base_model_nested <.0001Nested model is better than just pile but not different from chamber model in the previous comparison.
# for auto correlation we will use DOY - but unscaled to avoid negative values:
Scaled_Cleaned_Co2_data %>%
group_by(Chamber_Corrected) %>%
arrange(DOY, .by_group = TRUE) %>%
summarise(repeats = any(duplicated(DOY)))## # A tibble: 6 × 2
## Chamber_Corrected repeats
## <fct> <lgl>
## 1 C1 TRUE
## 2 C2 FALSE
## 3 C3 TRUE
## 4 C4 FALSE
## 5 C5 TRUE
## 6 C6 FALSE# while dates don't repeat, for the scaled date there are some repeats which we have to get rid of for autocorrelation.
Scaled_Cleaned_Co2_data <- Scaled_Cleaned_Co2_data %>%
group_by(Chamber_Corrected, un_scaled_DOY) %>%
filter(row_number() == 1) %>% # Keep only first occurrence
ungroup()
Scaled_Cleaned_Co2_data <- Scaled_Cleaned_Co2_data %>%
mutate(un_scaled_DOY = as.numeric(unlist(un_scaled_DOY)))
Scaled_Cleaned_Co2_data <- Scaled_Cleaned_Co2_data %>%
arrange(Chamber_Corrected, un_scaled_DOY) %>%
group_by(Chamber_Corrected) %>%
mutate(measurement_index = row_number()) %>%
ungroup()
# Add AR(1) autocorrelation structure based on DOY (or time variable) - by chamber not as a whole.
lme_ac <- update(
lme_base_model,
correlation = corAR1(form = ~ measurement_index | Chamber_Corrected)
)
# letting variance change by pile:
lme_var_pile <- update(lme_base_model, weights = varIdent(form = ~1 | Pile))
# letting variance change by chamber:
lme_var_chamber <- update(lme_base_model, weights = varIdent(form = ~1 | Chamber_Corrected))
# Compare the models
BIC(lme_base_model,lme_ac, lme_var_pile, lme_var_chamber)## Warning in BIC.default(lme_base_model, lme_ac, lme_var_pile, lme_var_chamber):
## models are not all fitted to the same number of observations## df BIC
## lme_base_model 14 11008.131
## lme_ac 15 -2336.744
## lme_var_pile 15 10954.224
## lme_var_chamber 19 10300.969# the BIC of the adding AC improves alot. Out of the adding different variances, it only improves the model fit with chamber.
# combining ac and var
lme_ac_var <- update(
lme_ac,
weights = varIdent(form = ~1 | Chamber_Corrected),
control = lmeControl(maxIter = 100, msMaxIter = 100, niterEM = 50)
)
# comparing the added ac_var model
BIC(lme_base_model,lme_var_chamber, lme_ac,lme_ac_var)## Warning in BIC.default(lme_base_model, lme_var_chamber, lme_ac, lme_ac_var):
## models are not all fitted to the same number of observations## df BIC
## lme_base_model 14 11008.131
## lme_var_chamber 19 10300.969
## lme_ac 15 -2336.744
## lme_ac_var 20 -2939.987summary(lme_ac_var)## Linear mixed-effects model fit by REML
## Data: Scaled_Cleaned_Co2_data
## AIC BIC logLik
## -3081.721 -2939.987 1560.861
##
## Random effects:
## Formula: ~1 | Chamber_Corrected
## (Intercept) Residual
## StdDev: 0.3087683 0.3578189
##
## Correlation Structure: AR(1)
## Formula: ~measurement_index | Chamber_Corrected
## Parameter estimate(s):
## Phi
## 0.8977519
## Variance function:
## Structure: Different standard deviations per stratum
## Formula: ~1 | Chamber_Corrected
## Parameter estimates:
## C1 C2 C3 C4 C5 C6
## 1.000000 1.484071 1.082007 1.182925 1.533261 1.696571
## Fixed effects: FCO2_DRY.LIN ~ Pile + SWC_1.initial_value + TS_1.initial_value + TMP + BAR + BulkDensity + Chamber_Elevation + DaysSinceLastTurn + TotalTurns + HMD + TA.range
## Value Std.Error DF t-value p-value
## (Intercept) -0.17252358 0.18091143 8833 -0.953636 0.3403
## PileE 0.05714406 0.25812180 4 0.221384 0.8356
## SWC_1.initial_value -0.04208342 0.01379860 8833 -3.049832 0.0023
## TS_1.initial_value 0.16516349 0.01461859 8833 11.298185 0.0000
## TMP 0.01642156 0.00672217 8833 2.442896 0.0146
## BAR -0.00328480 0.00226722 8833 -1.448821 0.1474
## BulkDensity -0.01959307 0.02527243 8833 -0.775275 0.4382
## Chamber_Elevation 0.06630316 0.01866817 8833 3.551668 0.0004
## DaysSinceLastTurn -0.08547129 0.01211401 8833 -7.055572 0.0000
## TotalTurns -0.07363091 0.02893257 8833 -2.544914 0.0109
## HMD -0.01370454 0.00469272 8833 -2.920382 0.0035
## TA.range 0.02163992 0.00506373 8833 4.273515 0.0000
## Correlation:
## (Intr) PileE SWC_1. TS_1._ TMP BAR BlkDns Chmb_E
## PileE -0.706
## SWC_1.initial_value 0.025 -0.041
## TS_1.initial_value -0.006 0.013 -0.139
## TMP -0.002 0.004 0.012 -0.028
## BAR -0.002 0.004 -0.082 0.052 0.180
## BulkDensity 0.066 -0.099 0.093 0.184 0.008 0.003
## Chamber_Elevation 0.004 -0.003 0.010 0.103 0.032 0.009 0.013
## DaysSinceLastTurn -0.005 0.006 0.152 0.002 0.038 0.026 -0.117 0.152
## TotalTurns -0.068 0.103 0.108 0.078 0.085 0.020 -0.622 0.046
## HMD 0.000 0.000 0.009 -0.040 0.680 0.173 0.027 -0.006
## TA.range 0.005 -0.006 -0.015 -0.074 0.078 0.002 -0.013 0.035
## DysSLT TtlTrn HMD
## PileE
## SWC_1.initial_value
## TS_1.initial_value
## TMP
## BAR
## BulkDensity
## Chamber_Elevation
## DaysSinceLastTurn
## TotalTurns 0.373
## HMD 0.011 0.029
## TA.range 0.056 0.023 -0.016
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -2.8606945 -0.6725356 -0.1511909 0.5428699 5.5879904
##
## Number of Observations: 8849
## Number of Groups: 6plot(lme_ac_var)
qqnorm(resid(lme_ac_var))
qqline(resid(lme_ac_var))
# Points of Further Discussion
# Most importantly - what exactly do we want to say about management - rate of emissions do not seems to be correlated to Pile, but more on other environmental variables at a give point in time, but relationship might be shorter term (days since last turn) when there is still plenty of substrate ? IDK. Clearly pile does not explain any of the variation in our model - yet pile is not the only variable that is related to management. How do we convey the impact of turning visually and create a narrative which details the impact of turning rate.