PSU 2022 Spatial Analysis: SpATS Spatial Correction for Both P Treatments
Fausto Rodriguez
April 21, 2026
Overview
This notebook fits a SpATS PSANOVA surface per phenotype and
treatment and returns spatially-corrected plot-level values for
each phenotype. The approach uses statgenHTP::fitModels()
(which wraps SpATS when provided a TimePoints object) and
getCorrected() to extract corrected plot-level
observations. The workflow:
- Load and clean data
- Create
plot_data(plot-level means) - Build TimePoints objects for both +P and -P treatments
- Loop over usable phenotypes and apply
fitModels()→getCorrected()for each treatment - Merge corrected values from both treatments into a combined dataset
- Perform second-stage treatment × genotype analysis
- Save corrected table for downstream analyses
Warning: fitting many SpATS models can be computationally intensive. Run interactively on a machine with sufficient RAM/CPU.
1. Load libraries
library(tidyverse)
library(dplyr)
library(ggplot2)
library(ggtext)
library(ggfx)
library(nlme) # Mixed-effects models (if needed)
library(gstat) # Variograms
library(rstatix)
library(ggpubr)
library(VIM)
library(knitr)
library(SpATS) # SpATS package (used internally by statgenHTP)
library(statgenHTP) # Provides fitModels() and getCorrected()2. Load and prepare data
plant_csv <- file.path(paths$data, "22_NCS_PSU_LANGEBIO_FIELDS_PSU_P_field.csv")
ear_csv <- file.path(paths$data, "22_NCS_PSU_LANGEBIO_FIELDS_PSU_P_field_ear_pheno.csv")
field_data_raw <- read.csv(plant_csv, na.strings = c("", "#N/A", "NA"))
ear_data_raw <- read.csv(ear_csv, na.strings = c("", "n/a", "NA"), skip = 1)
field_data <- field_data_raw %>%
filter(P22. >= 3004, P22. <= 4192) %>%
rename(
rowid = P22.,
Genotype = Who.What,
PH = Height_Anthesis,
STW40 = X40_DAP_dw,
STW50 = X50_DAP_dw,
STW60 = X60_DAP_dw,
STWHV = harvest_dw
) %>%
mutate(
# Set +P as reference level (first level)
Treatment = factor(Treatment, levels = c("HighP", "LowP")),
Genotype = factor(Genotype),
Rep = as.factor(Rep)
) %>%
mutate(
Treatment = fct_recode(Treatment, "+P" = "HighP", "-P" = "LowP"),
# Set CTRL as reference genotype
Genotype = fct_recode(Genotype, "Inv4m" = "INV4M"),
Genotype = fct_relevel(Genotype, "CTRL")
)
ear_data <- ear_data_raw %>%
select(-description, -RK, -CC, -NIR, ear_rep = rep) %>%
rename(rowid = row) %>%
arrange(rowid) %>%
group_by(rowid) %>%
select(-ear_rep) %>%
summarise_all(mean, na.rm = TRUE) %>%
droplevels()
field_data_complete <- field_data %>%
left_join(ear_data, by = "rowid") %>%
# Bender 2013-style harvest index: grain biomass / total above-ground
# biomass = TKW / (TKW + STWHV). Bounded [0, 1], matches the mineral HI
# convention (grain content / total content) used downstream.
mutate(HI = TKW / (TKW + STWHV)) %>%
mutate(
Plot_Column = case_when(
Plot == "PlotVIII" ~ Plot_Column + 10,
Plot == "PlotVI" ~ Plot_Column,
TRUE ~ Plot_Column
)
) %>%
filter(!is.na(PH)) %>%
droplevels()
# Load raw ICP-OES ionome data (one row per tissue x rowid) and reshape to
# one row per plot with {mineral}_grain and {mineral}_stover columns. Apply
# the same concentration outlier masks used in PSU2022_ionome.Rmd so that
# content derivations stay consistent with the concentration analyses.
ionome_csv <- file.path(paths$data, "PSU_inv4m_ionome_all.csv")
reliable_minerals <- c("Ca", "Fe", "K", "Mg", "Mn", "P", "S", "Zn")
ionome_wide <- read.csv(ionome_csv, na.strings = c("", "NA")) %>%
mutate(tissue = recode(tissue, "stalk" = "stover", "seed" = "grain")) %>%
select(tissue, rowid, all_of(reliable_minerals)) %>%
pivot_wider(
id_cols = rowid,
names_from = tissue,
values_from = all_of(reliable_minerals),
names_glue = "{.value}_{tissue}"
)
ionome_wide$Fe_stover[ionome_wide$Fe_stover > 800] <- NA
ionome_wide$Mn_stover[ionome_wide$Mn_stover > 75] <- NA
field_data_complete <- field_data_complete %>%
left_join(ionome_wide, by = "rowid")
# Derive per-plant mineral content and mineral harvest index for each
# reliable mineral at the raw plot level. These flow through the SpATS
# loop below as first-class traits (same way HI does), so that the
# corrected content variables are detrended as a single trait rather than
# reconstructed from two independently-corrected inputs.
for (ion in reliable_minerals) {
grain_col <- paste0(ion, "_grain")
stover_col <- paste0(ion, "_stover")
field_data_complete[[paste0(ion, "_grain_mg")]] <-
field_data_complete$TKW * field_data_complete[[grain_col]] / 1000
field_data_complete[[paste0(ion, "_stover_mg")]] <-
field_data_complete$STWHV * field_data_complete[[stover_col]] / 1000
field_data_complete[[paste0(ion, "_total_mg")]] <-
field_data_complete[[paste0(ion, "_grain_mg")]] +
field_data_complete[[paste0(ion, "_stover_mg")]]
field_data_complete[[paste0(ion, "_HI")]] <-
field_data_complete[[paste0(ion, "_grain_mg")]] /
field_data_complete[[paste0(ion, "_total_mg")]]
field_data_complete[[paste0(ion, "_ratio")]] <-
field_data_complete[[grain_col]] /
field_data_complete[[stover_col]]
}
# Define phenotypes
base_phenotypes <- c("PH", "DTA", "DTS", "STW40", "STW50", "STW60", "STWHV")
ear_phenotypes <- c("EW", "EL", "ED", "KRN", "CD", "TKW", "KW50", "TKN", "HI")
ionome_conc_phenotypes <- c(
paste0(reliable_minerals, "_grain"),
paste0(reliable_minerals, "_stover")
)
ionome_content_phenotypes <- c(
paste0(reliable_minerals, "_grain_mg"),
paste0(reliable_minerals, "_stover_mg"),
paste0(reliable_minerals, "_total_mg"),
paste0(reliable_minerals, "_HI"),
paste0(reliable_minerals, "_ratio")
)
all_phenotypes <- unique(c(base_phenotypes, ear_phenotypes,
ionome_conc_phenotypes,
ionome_content_phenotypes))
plot_data <- field_data_complete %>%
group_by(rowid, Rep, Plot_Row, Plot_Column, Genotype, Treatment) %>%
summarise(
across(all_of(all_phenotypes), ~ mean(.x, na.rm = TRUE)),
n_plants = n(),
.groups = 'drop'
) %>%
mutate(across(all_of(all_phenotypes), ~ ifelse(is.nan(.x), NA, .x)))
plot_data$Plot_Row <- as.numeric(plot_data$Plot_Row)
plot_data$Plot_Column <- as.numeric(plot_data$Plot_Column)
plot_data$Rep <- as.factor(plot_data$Rep)
cat("Plot-level summary dimensions:", dim(plot_data), "\n")## Plot-level summary dimensions: 127 79cat("Average plants per plot:", round(mean(plot_data$n_plants, na.rm = TRUE), 1), "\n")## Average plants per plot: 13. Identify usable phenotypes (<50% missing)
missing_summary <- plot_data %>%
select(all_of(all_phenotypes)) %>%
summarise(across(everything(), ~ sum(is.na(.)))) %>%
pivot_longer(everything(), names_to = "phenotype", values_to = "missing_count") %>%
mutate(
total_obs = nrow(plot_data),
missing_pct = round(missing_count / total_obs * 100, 1),
available_n = total_obs - missing_count
) %>%
arrange(missing_pct)
print(missing_summary)## # A tibble: 72 × 5
## phenotype missing_count total_obs missing_pct available_n
## <chr> <int> <int> <dbl> <int>
## 1 PH 0 127 0 127
## 2 DTA 0 127 0 127
## 3 DTS 0 127 0 127
## 4 STW40 2 127 1.6 125
## 5 STWHV 2 127 1.6 125
## 6 Ca_stover 2 127 1.6 125
## 7 K_stover 2 127 1.6 125
## 8 Mg_stover 2 127 1.6 125
## 9 P_stover 2 127 1.6 125
## 10 S_stover 2 127 1.6 125
## # ℹ 62 more rowsusable_phenotypes <- missing_summary %>%
filter(missing_pct < 50) %>%
pull(phenotype)
cat("Usable phenotypes (<50% missing):", paste(usable_phenotypes, collapse = ", "), "\n")## Usable phenotypes (<50% missing): PH, DTA, DTS, STW40, STWHV, Ca_stover, K_stover, Mg_stover, P_stover, S_stover, Zn_stover, Ca_stover_mg, K_stover_mg, Mg_stover_mg, P_stover_mg, S_stover_mg, Zn_stover_mg, Fe_stover, Fe_stover_mg, Mn_stover, Mn_stover_mg, STW50, STW60, EW, EL, ED, KRN, CD, TKW, KW50, TKN, HI, Ca_grain, Fe_grain, K_grain, Mg_grain, Mn_grain, P_grain, S_grain, Zn_grain, Ca_grain_mg, Fe_grain_mg, K_grain_mg, Mg_grain_mg, Mn_grain_mg, P_grain_mg, S_grain_mg, Zn_grain_mg, Ca_total_mg, Fe_total_mg, K_total_mg, Mg_total_mg, P_total_mg, S_total_mg, Zn_total_mg, Ca_HI, Fe_HI, K_HI, Mg_HI, P_HI, S_HI, Zn_HI, Ca_ratio, Fe_ratio, K_ratio, Mg_ratio, P_ratio, S_ratio, Zn_ratio, Mn_total_mg, Mn_HI, Mn_ratio4. Spatial correction for +P treatment
# Prepare data for TimePoints
single_time <- plot_data %>%
mutate(
Genotype = as.character(Genotype),
Treatment = as.character(Treatment),
Rep = as.character(Rep),
Plot_Row = as.numeric(Plot_Row),
Plot_Column = as.numeric(Plot_Column)
) %>%
mutate(time = as.POSIXct("2022-05-22 10:00:00"))
# Create TimePoints object for +P treatment
tp_data_plus_p <- createTimePoints(
dat = single_time %>% filter(Treatment == "+P"),
experimentName = "PSU2022_PlusP",
genotype = "Genotype",
timePoint = "time",
plotId = "rowid",
repId = "Rep",
rowNum = "Plot_Row",
colNum = "Plot_Column",
checkGenotypes = c("NUE", "B73"),
addCheck = TRUE
)
summary(tp_data_plus_p)
# Initialize lists for +P treatment
corrected_list_plus_p <- list()
fit_info_plus_p <- list()
all_models_plus_p <- list()
for (phen in usable_phenotypes) {
message("Fitting SpATS model for +P treatment, phenotype: ", phen)
# Check data availability
no_nas_plus_p <- single_time %>%
filter(Treatment == "+P", !is.na(.data[[phen]]))
if (nrow(no_nas_plus_p) < 10) {
message("Insufficient data for ", phen, " in +P treatment (n=", nrow(no_nas_plus_p), ")")
next
}
# Fit model
fit_try <- try(
fitModels(
TP = tp_data_plus_p,
trait = phen,
timePoints = 1,
what = "fixed"
),
silent = TRUE
)
if (inherits(fit_try, "try-error")) {
message("Fit failed for ", phen, " in +P treatment: ", fit_try)
next
}
all_models_plus_p[[phen]] <- fit_try
# Extract corrected values
corr <- try(getCorrected(fit_try), silent = TRUE)
if (inherits(corr, "try-error") || is.null(corr)) {
message("getCorrected() failed for ", phen, " in +P treatment")
next
}
corr <- corr %>% select(-phen)
corrected_list_plus_p[[phen]] <- corr
fit_info_plus_p[[phen]] <- list(n_obs = nrow(no_nas_plus_p), status = "success")
}
# Combine all corrected data for +P treatment
sp_corrected_plus_p <- lapply(
corrected_list_plus_p,
function(corr_df) {
corr_df %>%
select(timeNumber, timePoint, genotype, rowId, colId, plotId)
}
) %>%
bind_rows() %>%
arrange(plotId) %>%
distinct()
# Join each corrected phenotype for +P
for (phen in names(corrected_list_plus_p)) {
phen_corr <- paste0(phen, "_corr")
corr_df <- corrected_list_plus_p[[phen]] %>%
select(plotId, all_of(phen_corr))
sp_corrected_plus_p <- left_join(
sp_corrected_plus_p,
corr_df,
by = "plotId"
)
}
# Clean column names (remove _corr suffix)
colnames(sp_corrected_plus_p) <- gsub(
"_corr", "",
colnames(sp_corrected_plus_p)
)
# Add treatment information
sp_corrected_plus_p <- sp_corrected_plus_p %>%
inner_join(
tp_data_plus_p[[1]] %>% select(plotId, Treatment),
by = "plotId"
) %>%
select(timeNumber:genotype, Treatment, everything())
cat("Final sp_corrected_plus_p dimensions:", dim(sp_corrected_plus_p), "\n")5. Spatial correction for -P treatment
# Create TimePoints object for -P treatment
tp_data_minus_p <- createTimePoints(
dat = single_time %>% filter(Treatment == "-P"),
experimentName = "PSU2022_MinusP",
genotype = "Genotype",
timePoint = "time",
plotId = "rowid",
repId = "Rep",
rowNum = "Plot_Row",
colNum = "Plot_Column",
checkGenotypes = c("NUE", "B73"),
addCheck = TRUE
)
summary(tp_data_minus_p)
# Initialize lists for -P treatment
corrected_list_minus_p <- list()
fit_info_minus_p <- list()
all_models_minus_p <- list()
for (phen in usable_phenotypes) {
message("Fitting SpATS model for -P treatment, phenotype: ", phen)
# Check if we have enough data for this phenotype in -P
no_nas_minus_p <- single_time %>%
filter(Treatment == "-P", !is.na(.data[[phen]]))
if (nrow(no_nas_minus_p) < 10) {
message("Insufficient data for ", phen, " in -P treatment (n=", nrow(no_nas_minus_p), ")")
next
}
# Fit model
fit_try <- try(
fitModels(
TP = tp_data_minus_p,
trait = phen,
timePoints = 1,
what = "fixed"
),
silent = TRUE
)
if (inherits(fit_try, "try-error")) {
message("Fit failed for ", phen, " in -P treatment: ", fit_try)
next
}
all_models_minus_p[[phen]] <- fit_try
# Extract corrected values
corr <- try(getCorrected(fit_try), silent = TRUE)
if (inherits(corr, "try-error") || is.null(corr)) {
message("getCorrected() failed for ", phen, " in -P treatment")
next
}
corr <- corr %>% select(-phen)
corrected_list_minus_p[[phen]] <- corr
fit_info_minus_p[[phen]] <- list(n_obs = nrow(no_nas_minus_p), status = "success")
}
# Combine all corrected data for -P treatment
sp_corrected_minus_p <- lapply(
corrected_list_minus_p,
function(corr_df) {
corr_df %>%
select(timeNumber, timePoint, genotype, rowId, colId, plotId)
}
) %>%
bind_rows() %>%
arrange(plotId) %>%
distinct()
# Join each corrected phenotype for -P
for (phen in names(corrected_list_minus_p)) {
phen_corr <- paste0(phen, "_corr")
corr_df <- corrected_list_minus_p[[phen]] %>%
select(plotId, all_of(phen_corr))
sp_corrected_minus_p <- left_join(
sp_corrected_minus_p,
corr_df,
by = "plotId"
)
}
# Clean column names (remove _corr suffix)
colnames(sp_corrected_minus_p) <- gsub(
"_corr", "",
colnames(sp_corrected_minus_p)
)
# Add treatment information
sp_corrected_minus_p <- sp_corrected_minus_p %>%
inner_join(
tp_data_minus_p[[1]] %>% select(plotId, Treatment),
by = "plotId"
) %>%
select(timeNumber:genotype, Treatment, everything())
cat("Final sp_corrected_minus_p dimensions:", dim(sp_corrected_minus_p), "\n")6. Merge both spatially corrected datasets
# Combine +P and -P corrected data
sp_corrected_combined <- bind_rows(
sp_corrected_plus_p, # +P treatment
sp_corrected_minus_p # -P treatment
) %>%
arrange(plotId)
cat("Combined sp_corrected dimensions:", dim(sp_corrected_combined), "\n")## Combined sp_corrected dimensions: 127 79cat("Treatment distribution:\n")## Treatment distribution:table(sp_corrected_combined$Treatment)##
## -P +P
## 63 64# Verify we have data for both treatments
treatment_summary <- sp_corrected_combined %>%
group_by(Treatment) %>%
summarise(
n_plots = n(),
n_genotypes = n_distinct(genotype),
.groups = "drop"
)
print(treatment_summary)## # A tibble: 2 × 3
## Treatment n_plots n_genotypes
## <chr> <int> <int>
## 1 +P 64 4
## 2 -P 63 4# Check for overlapping phenotypes between treatments
plus_p_phenos <- names(sp_corrected_plus_p)[!names(sp_corrected_plus_p) %in%
c("timeNumber", "timePoint", "genotype",
"rowId", "colId", "plotId", "Treatment")]
minus_p_phenos <- names(sp_corrected_minus_p)[!names(sp_corrected_minus_p) %in%
c("timeNumber", "timePoint", "genotype",
"rowId", "colId", "plotId", "Treatment")]
common_phenos <- intersect(plus_p_phenos, minus_p_phenos)
cat("Phenotypes corrected in both treatments:", paste(common_phenos, collapse = ", "), "\n")## Phenotypes corrected in both treatments: PH, DTA, DTS, STW40, STWHV, Ca_stover, K_stover, Mg_stover, P_stover, S_stover, Zn_stover, Ca_stover_mg, K_stover_mg, Mg_stover_mg, P_stover_mg, S_stover_mg, Zn_stover_mg, Fe_stover, Fe_stover_mg, Mn_stover, Mn_stover_mg, STW50, STW60, EW, EL, ED, KRN, CD, TKW, KW50, TKN, HI, Ca_grain, Fe_grain, K_grain, Mg_grain, Mn_grain, P_grain, S_grain, Zn_grain, Ca_grain_mg, Fe_grain_mg, K_grain_mg, Mg_grain_mg, Mn_grain_mg, P_grain_mg, S_grain_mg, Zn_grain_mg, Ca_total_mg, Fe_total_mg, K_total_mg, Mg_total_mg, P_total_mg, S_total_mg, Zn_total_mg, Ca_HI, Fe_HI, K_HI, Mg_HI, P_HI, S_HI, Zn_HI, Ca_ratio, Fe_ratio, K_ratio, Mg_ratio, P_ratio, S_ratio, Zn_ratio, Mn_total_mg, Mn_HI, Mn_ratio7. Spatial correction diagnostic plots
# Example spatial diagnostic plots for key phenotypes
key_phenotypes <- intersect(c("PH", "HI", "DTS", "TKW"), common_phenos)
for (phen in key_phenotypes[1:2]) { # Limit to first 2 to avoid too many plots
if (phen %in% names(all_models_plus_p)) {
cat("\n### +P Treatment:", phen, "\n")
plot(all_models_plus_p[[phen]], timePoints = 1)
plot(all_models_plus_p[[phen]], timePoints = 1, plotType = "spatial", spaTrend = "raw")
}
if (phen %in% names(all_models_minus_p)) {
cat("\n### -P Treatment:", phen, "\n")
plot(all_models_minus_p[[phen]], timePoints = 1)
plot(all_models_minus_p[[phen]], timePoints = 1, plotType = "spatial", spaTrend = "raw")
}
}##
## ### +P Treatment: PH
##
## ### -P Treatment: PH
##
## ### +P Treatment: HI
##
## ### -P Treatment: HI
8. Second-stage analysis: Treatment and Genotype × Treatment effects
# Filter to NILs for treatment × genotype analysis
nil_corrected_combined <- sp_corrected_combined %>%
filter(genotype %in% c("Inv4m", "CTRL")) %>%
mutate(
# Ensure proper factor levels with +P and CTRL as references
Treatment = factor(Treatment, levels = c("+P", "-P")),
genotype = factor(genotype, levels = c("CTRL", "Inv4m"))
) %>%
droplevels()
cat("NIL data dimensions:", dim(nil_corrected_combined), "\n")## NIL data dimensions: 64 79cat("Treatment × Genotype combinations:\n")## Treatment × Genotype combinations:table(nil_corrected_combined$Treatment, nil_corrected_combined$genotype)##
## CTRL Inv4m
## +P 16 16
## -P 16 16# Extract model coefficients for all common phenotypes
model_coefficients <- list()
cat("\n=== Treatment and Genotype × Treatment Models ===\n")##
## === Treatment and Genotype × Treatment Models ===for (phen in common_phenos) {
if (sum(!is.na(nil_corrected_combined[[phen]])) < 20) {
message("Skipping ", phen, " - insufficient data")
next
}
model_formula <- as.formula(paste(phen, "~ genotype * Treatment"))
lm_result <- try(
lm(model_formula, data = nil_corrected_combined),
silent = TRUE
)
if (!inherits(lm_result, "try-error")) {
# Extract coefficients
coef_summary <- summary(lm_result)$coefficients
coef_df <- data.frame(
trait = phen,
term = rownames(coef_summary),
estimate = coef_summary[, "Estimate"],
std_error = coef_summary[, "Std. Error"],
t_value = coef_summary[, "t value"],
p_value = coef_summary[, "Pr(>|t|)"],
stringsAsFactors = FALSE
)
model_coefficients[[phen]] <- coef_df
cat("\n=== ", phen, " ===\n")
print(summary(lm_result))
}
}##
## === PH ===
##
## Call:
## lm(formula = model_formula, data = nil_corrected_combined)
##
## Residuals:
## Min 1Q Median 3Q Max
## -9.8349 -2.7792 0.4735 2.8230 9.6399
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 161.658 1.047 154.406 < 2e-16 ***
## genotypeInv4m 5.570 1.481 3.762 0.000385 ***
## Treatment-P 2.588 1.481 1.748 0.085629 .
## genotypeInv4m:Treatment-P 1.677 2.094 0.801 0.426346
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.188 on 60 degrees of freedom
## Multiple R-squared: 0.4486, Adjusted R-squared: 0.421
## F-statistic: 16.27 on 3 and 60 DF, p-value: 7.49e-08
##
##
## === DTA ===
##
## Call:
## lm(formula = model_formula, data = nil_corrected_combined)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.35736 -0.70516 -0.01931 0.60844 2.88214
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 72.1200 0.2599 277.458 < 2e-16 ***
## genotypeInv4m -1.4778 0.3676 -4.020 0.000165 ***
## Treatment-P 3.4329 0.3676 9.339 2.67e-13 ***
## genotypeInv4m:Treatment-P 0.3374 0.5199 0.649 0.518774
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.04 on 60 degrees of freedom
## Multiple R-squared: 0.784, Adjusted R-squared: 0.7732
## F-statistic: 72.59 on 3 and 60 DF, p-value: < 2.2e-16
##
##
## === DTS ===
##
## Call:
## lm(formula = model_formula, data = nil_corrected_combined)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.30539 -0.58763 0.05908 0.52957 2.51458
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 73.6544 0.2277 323.494 < 2e-16 ***
## genotypeInv4m -1.2842 0.3220 -3.988 0.000183 ***
## Treatment-P 3.0621 0.3220 9.510 1.38e-13 ***
## genotypeInv4m:Treatment-P 0.7087 0.4554 1.556 0.124886
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.9107 on 60 degrees of freedom
## Multiple R-squared: 0.8028, Adjusted R-squared: 0.7929
## F-statistic: 81.42 on 3 and 60 DF, p-value: < 2.2e-16
##
##
## === STW40 ===
##
## Call:
## lm(formula = model_formula, data = nil_corrected_combined)
##
## Residuals:
## Min 1Q Median 3Q Max
## -8.338 -1.303 -0.032 1.450 6.476
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 12.0737 0.7239 16.679 < 2e-16 ***
## genotypeInv4m -2.9345 1.0076 -2.912 0.00506 **
## Treatment-P -6.5170 1.0076 -6.468 2.14e-08 ***
## genotypeInv4m:Treatment-P 2.8619 1.4134 2.025 0.04741 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.804 on 59 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.5015, Adjusted R-squared: 0.4761
## F-statistic: 19.78 on 3 and 59 DF, p-value: 5.389e-09
##
##
## === STWHV ===
##
## Call:
## lm(formula = model_formula, data = nil_corrected_combined)
##
## Residuals:
## Min 1Q Median 3Q Max
## -27.370 -7.546 1.087 5.429 43.449
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 106.989 3.433 31.163 < 2e-16 ***
## genotypeInv4m -11.655 4.701 -2.479 0.0161 *
## Treatment-P -22.407 4.701 -4.766 1.3e-05 ***
## genotypeInv4m:Treatment-P 7.421 6.537 1.135 0.2609
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 12.85 on 58 degrees of freedom
## (2 observations deleted due to missingness)
## Multiple R-squared: 0.3988, Adjusted R-squared: 0.3677
## F-statistic: 12.82 on 3 and 58 DF, p-value: 1.557e-06
##
##
## === Ca_stover ===
##
## Call:
## lm(formula = model_formula, data = nil_corrected_combined)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1072.41 -163.08 -42.64 313.41 809.12
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4534.32 103.04 44.006 < 2e-16 ***
## genotypeInv4m -411.42 141.09 -2.916 0.00503 **
## Treatment-P 65.12 141.09 0.462 0.64612
## genotypeInv4m:Treatment-P -152.75 196.18 -0.779 0.43936
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 385.5 on 58 degrees of freedom
## (2 observations deleted due to missingness)
## Multiple R-squared: 0.3065, Adjusted R-squared: 0.2706
## F-statistic: 8.545 on 3 and 58 DF, p-value: 8.674e-05
##
##
## === K_stover ===
##
## Call:
## lm(formula = model_formula, data = nil_corrected_combined)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4079.4 -1175.1 -89.2 1185.2 4493.6
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 11126.9 468.8 23.734 <2e-16 ***
## genotypeInv4m -577.1 641.9 -0.899 0.372
## Treatment-P -526.2 641.9 -0.820 0.416
## genotypeInv4m:Treatment-P 101.2 892.6 0.113 0.910
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1754 on 58 degrees of freedom
## (2 observations deleted due to missingness)
## Multiple R-squared: 0.04045, Adjusted R-squared: -0.009183
## F-statistic: 0.815 on 3 and 58 DF, p-value: 0.4908
##
##
## === Mg_stover ===
##
## Call:
## lm(formula = model_formula, data = nil_corrected_combined)
##
## Residuals:
## Min 1Q Median 3Q Max
## -283.30 -80.86 4.04 76.83 354.14
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1560.4620 35.8895 43.480 <2e-16 ***
## genotypeInv4m 28.6282 49.1437 0.583 0.562
## Treatment-P 0.5654 49.1437 0.012 0.991
## genotypeInv4m:Treatment-P -46.5051 68.3316 -0.681 0.499
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 134.3 on 58 degrees of freedom
## (2 observations deleted due to missingness)
## Multiple R-squared: 0.0163, Adjusted R-squared: -0.03458
## F-statistic: 0.3203 on 3 and 58 DF, p-value: 0.8106
##
##
## === P_stover ===
##
## Call:
## lm(formula = model_formula, data = nil_corrected_combined)
##
## Residuals:
## Min 1Q Median 3Q Max
## -542.93 -121.52 -9.96 95.48 610.64
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2313.06 61.77 37.444 <2e-16 ***
## genotypeInv4m 115.05 84.59 1.360 0.179
## Treatment-P -1592.10 84.59 -18.822 <2e-16 ***
## genotypeInv4m:Treatment-P -80.49 117.62 -0.684 0.496
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 231.1 on 58 degrees of freedom
## (2 observations deleted due to missingness)
## Multiple R-squared: 0.9306, Adjusted R-squared: 0.927
## F-statistic: 259.3 on 3 and 58 DF, p-value: < 2.2e-16
##
##
## === S_stover ===
##
## Call:
## lm(formula = model_formula, data = nil_corrected_combined)
##
## Residuals:
## Min 1Q Median 3Q Max
## -126.052 -36.988 1.733 43.026 111.986
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 938.94 15.35 61.169 < 2e-16 ***
## genotypeInv4m 18.24 21.02 0.868 0.38919
## Treatment-P 113.06 21.02 5.379 1.4e-06 ***
## genotypeInv4m:Treatment-P -93.81 29.23 -3.210 0.00217 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 57.43 on 58 degrees of freedom
## (2 observations deleted due to missingness)
## Multiple R-squared: 0.3748, Adjusted R-squared: 0.3425
## F-statistic: 11.59 on 3 and 58 DF, p-value: 4.7e-06
##
##
## === Zn_stover ===
##
## Call:
## lm(formula = model_formula, data = nil_corrected_combined)
##
## Residuals:
## Min 1Q Median 3Q Max
## -7.1408 -1.5863 -0.0824 2.1290 5.8840
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 23.5778 0.8185 28.805 < 2e-16 ***
## genotypeInv4m -1.5686 1.1208 -1.399 0.167
## Treatment-P 6.8546 1.1208 6.116 8.78e-08 ***
## genotypeInv4m:Treatment-P -1.8365 1.5585 -1.178 0.243
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.063 on 58 degrees of freedom
## (2 observations deleted due to missingness)
## Multiple R-squared: 0.5505, Adjusted R-squared: 0.5273
## F-statistic: 23.68 on 3 and 58 DF, p-value: 3.919e-10
##
##
## === Ca_stover_mg ===
##
## Call:
## lm(formula = model_formula, data = nil_corrected_combined)
##
## Residuals:
## Min 1Q Median 3Q Max
## -176.714 -26.418 0.996 35.070 161.180
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 485.40 15.04 32.268 < 2e-16 ***
## genotypeInv4m -93.64 20.60 -4.546 2.83e-05 ***
## Treatment-P -97.38 20.60 -4.727 1.49e-05 ***
## genotypeInv4m:Treatment-P 30.30 28.64 1.058 0.295
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 56.29 on 58 degrees of freedom
## (2 observations deleted due to missingness)
## Multiple R-squared: 0.5142, Adjusted R-squared: 0.4891
## F-statistic: 20.46 on 3 and 58 DF, p-value: 3.623e-09
##
##
## === K_stover_mg ===
##
## Call:
## lm(formula = model_formula, data = nil_corrected_combined)
##
## Residuals:
## Min 1Q Median 3Q Max
## -498.14 -128.76 -20.39 95.79 693.62
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1202.25 61.36 19.594 < 2e-16 ***
## genotypeInv4m -213.04 84.02 -2.536 0.013940 *
## Treatment-P -298.85 84.02 -3.557 0.000755 ***
## genotypeInv4m:Treatment-P 116.59 116.82 0.998 0.322405
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 229.6 on 58 degrees of freedom
## (2 observations deleted due to missingness)
## Multiple R-squared: 0.2914, Adjusted R-squared: 0.2548
## F-statistic: 7.951 on 3 and 58 DF, p-value: 0.0001583
##
##
## === Mg_stover_mg ===
##
## Call:
## lm(formula = model_formula, data = nil_corrected_combined)
##
## Residuals:
## Min 1Q Median 3Q Max
## -42.745 -17.934 -1.323 11.758 67.189
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 167.221 6.111 27.364 < 2e-16 ***
## genotypeInv4m -15.315 8.368 -1.830 0.072359 .
## Treatment-P -34.842 8.368 -4.164 0.000105 ***
## genotypeInv4m:Treatment-P 7.330 11.635 0.630 0.531173
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 22.87 on 58 degrees of freedom
## (2 observations deleted due to missingness)
## Multiple R-squared: 0.3568, Adjusted R-squared: 0.3236
## F-statistic: 10.73 on 3 and 58 DF, p-value: 1.046e-05
##
##
## === P_stover_mg ===
##
## Call:
## lm(formula = model_formula, data = nil_corrected_combined)
##
## Residuals:
## Min 1Q Median 3Q Max
## -111.824 -15.521 -3.535 6.321 151.105
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 249.25 11.49 21.692 <2e-16 ***
## genotypeInv4m -19.23 15.73 -1.222 0.227
## Treatment-P -186.00 15.73 -11.821 <2e-16 ***
## genotypeInv4m:Treatment-P 16.54 21.88 0.756 0.453
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 42.99 on 58 degrees of freedom
## (2 observations deleted due to missingness)
## Multiple R-squared: 0.82, Adjusted R-squared: 0.8107
## F-statistic: 88.07 on 3 and 58 DF, p-value: < 2.2e-16
##
##
## === S_stover_mg ===
##
## Call:
## lm(formula = model_formula, data = nil_corrected_combined)
##
## Residuals:
## Min 1Q Median 3Q Max
## -34.941 -10.262 -0.297 7.540 41.479
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 101.054 4.144 24.386 <2e-16 ***
## genotypeInv4m -10.005 5.674 -1.763 0.0831 .
## Treatment-P -11.417 5.674 -2.012 0.0489 *
## genotypeInv4m:Treatment-P -1.352 7.890 -0.171 0.8645
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 15.51 on 58 degrees of freedom
## (2 observations deleted due to missingness)
## Multiple R-squared: 0.2195, Adjusted R-squared: 0.1791
## F-statistic: 5.436 on 3 and 58 DF, p-value: 0.002307
##
##
## === Zn_stover_mg ===
##
## Call:
## lm(formula = model_formula, data = nil_corrected_combined)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.11447 -0.18265 0.01989 0.29016 0.71211
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.50263 0.11512 21.740 < 2e-16 ***
## genotypeInv4m -0.44031 0.15763 -2.793 0.00706 **
## Treatment-P 0.04069 0.15763 0.258 0.79719
## genotypeInv4m:Treatment-P 0.06470 0.21917 0.295 0.76889
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4307 on 58 degrees of freedom
## (2 observations deleted due to missingness)
## Multiple R-squared: 0.2003, Adjusted R-squared: 0.159
## F-statistic: 4.843 on 3 and 58 DF, p-value: 0.004478
##
##
## === Fe_stover ===
##
## Call:
## lm(formula = model_formula, data = nil_corrected_combined)
##
## Residuals:
## Min 1Q Median 3Q Max
## -117.393 -38.014 -4.325 23.388 245.859
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 115.860 17.008 6.812 6.54e-09 ***
## genotypeInv4m -4.349 22.898 -0.190 0.850
## Treatment-P 30.994 22.898 1.354 0.181
## genotypeInv4m:Treatment-P -6.242 31.534 -0.198 0.844
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 61.32 on 57 degrees of freedom
## (3 observations deleted due to missingness)
## Multiple R-squared: 0.05736, Adjusted R-squared: 0.007751
## F-statistic: 1.156 on 3 and 57 DF, p-value: 0.3344
##
##
## === Fe_stover_mg ===
##
## Call:
## lm(formula = model_formula, data = nil_corrected_combined)
##
## Residuals:
## Min 1Q Median 3Q Max
## -11.2378 -3.3238 -0.7895 2.6821 23.1808
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 13.4675 1.5150 8.890 2.34e-12 ***
## genotypeInv4m -3.7180 2.0396 -1.823 0.0736 .
## Treatment-P -0.8055 2.0396 -0.395 0.6944
## genotypeInv4m:Treatment-P 1.7596 2.8089 0.626 0.5335
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.462 on 57 degrees of freedom
## (3 observations deleted due to missingness)
## Multiple R-squared: 0.07147, Adjusted R-squared: 0.0226
## F-statistic: 1.462 on 3 and 57 DF, p-value: 0.2344
##
##
## === Mn_stover ===
##
## Call:
## lm(formula = model_formula, data = nil_corrected_combined)
##
## Residuals:
## Min 1Q Median 3Q Max
## -13.5811 -3.5144 -0.4248 3.6682 17.6171
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 45.7021 1.7523 26.082 <2e-16 ***
## genotypeInv4m 0.7173 2.3940 0.300 0.7656
## Treatment-P 4.2071 2.3590 1.783 0.0799 .
## genotypeInv4m:Treatment-P -3.4645 3.2743 -1.058 0.2946
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 6.318 on 56 degrees of freedom
## (4 observations deleted due to missingness)
## Multiple R-squared: 0.06428, Adjusted R-squared: 0.01415
## F-statistic: 1.282 on 3 and 56 DF, p-value: 0.2894
##
##
## === Mn_stover_mg ===
##
## Call:
## lm(formula = model_formula, data = nil_corrected_combined)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.33317 -0.45456 -0.00399 0.44304 1.38830
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.85544 0.18475 26.281 <2e-16 ***
## genotypeInv4m -0.38530 0.25242 -1.526 0.1325
## Treatment-P -0.66245 0.24873 -2.663 0.0101 *
## genotypeInv4m:Treatment-P -0.01322 0.34522 -0.038 0.9696
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.6661 on 56 degrees of freedom
## (4 observations deleted due to missingness)
## Multiple R-squared: 0.2597, Adjusted R-squared: 0.22
## F-statistic: 6.548 on 3 and 56 DF, p-value: 0.0007118
##
##
## === STW50 ===
##
## Call:
## lm(formula = model_formula, data = nil_corrected_combined)
##
## Residuals:
## Min 1Q Median 3Q Max
## -12.2540 -2.2394 0.1365 2.5199 11.9398
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 35.984 1.159 31.060 <2e-16 ***
## genotypeInv4m -4.113 1.638 -2.510 0.0152 *
## Treatment-P -18.759 1.583 -11.851 <2e-16 ***
## genotypeInv4m:Treatment-P 4.620 2.222 2.079 0.0425 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.177 on 53 degrees of freedom
## (7 observations deleted due to missingness)
## Multiple R-squared: 0.8097, Adjusted R-squared: 0.7989
## F-statistic: 75.15 on 3 and 53 DF, p-value: < 2.2e-16
##
##
## === STW60 ===
##
## Call:
## lm(formula = model_formula, data = nil_corrected_combined)
##
## Residuals:
## Min 1Q Median 3Q Max
## -31.086 -7.884 -0.080 4.937 35.809
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 77.003 3.517 21.896 < 2e-16 ***
## genotypeInv4m -2.906 5.076 -0.572 0.570
## Treatment-P -33.343 4.973 -6.704 1.74e-08 ***
## genotypeInv4m:Treatment-P 10.807 6.941 1.557 0.126
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 12.68 on 50 degrees of freedom
## (10 observations deleted due to missingness)
## Multiple R-squared: 0.5715, Adjusted R-squared: 0.5458
## F-statistic: 22.23 on 3 and 50 DF, p-value: 2.768e-09
##
##
## === EW ===
##
## Call:
## lm(formula = model_formula, data = nil_corrected_combined)
##
## Residuals:
## Min 1Q Median 3Q Max
## -54.605 -10.703 1.946 12.888 21.022
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 86.460 4.788 18.056 <2e-16 ***
## genotypeInv4m 12.388 6.399 1.936 0.0587 .
## Treatment-P -2.487 6.304 -0.395 0.6949
## genotypeInv4m:Treatment-P -7.283 8.784 -0.829 0.4111
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 15.88 on 49 degrees of freedom
## (11 observations deleted due to missingness)
## Multiple R-squared: 0.1261, Adjusted R-squared: 0.07264
## F-statistic: 2.358 on 3 and 49 DF, p-value: 0.08307
##
##
## === EL ===
##
## Call:
## lm(formula = model_formula, data = nil_corrected_combined)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.6282 -0.2277 0.0289 0.4594 1.3351
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 14.7579 0.2905 50.804 <2e-16 ***
## genotypeInv4m 0.8936 0.3882 2.302 0.0256 *
## Treatment-P 0.8272 0.3824 2.163 0.0354 *
## genotypeInv4m:Treatment-P -1.0603 0.5329 -1.990 0.0522 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.9634 on 49 degrees of freedom
## (11 observations deleted due to missingness)
## Multiple R-squared: 0.1154, Adjusted R-squared: 0.06126
## F-statistic: 2.131 on 3 and 49 DF, p-value: 0.1083
##
##
## === ED ===
##
## Call:
## lm(formula = model_formula, data = nil_corrected_combined)
##
## Residuals:
## Min 1Q Median 3Q Max
## -15.3132 -0.6363 0.1322 1.1502 4.9629
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 39.04209 0.84882 45.996 <2e-16 ***
## genotypeInv4m -0.13902 1.13429 -0.123 0.903
## Treatment-P 0.07826 1.11753 0.070 0.944
## genotypeInv4m:Treatment-P -2.13106 1.55712 -1.369 0.177
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.815 on 49 degrees of freedom
## (11 observations deleted due to missingness)
## Multiple R-squared: 0.1073, Adjusted R-squared: 0.05259
## F-statistic: 1.962 on 3 and 49 DF, p-value: 0.1319
##
##
## === KRN ===
##
## Call:
## lm(formula = model_formula, data = nil_corrected_combined)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.3323 -0.3699 0.0803 0.6031 1.7838
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 16.73541 0.28661 58.390 <2e-16 ***
## genotypeInv4m -0.23275 0.38300 -0.608 0.546
## Treatment-P 0.09206 0.37734 0.244 0.808
## genotypeInv4m:Treatment-P -0.56062 0.52578 -1.066 0.292
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.9506 on 49 degrees of freedom
## (11 observations deleted due to missingness)
## Multiple R-squared: 0.1014, Adjusted R-squared: 0.04643
## F-statistic: 1.844 on 3 and 49 DF, p-value: 0.1515
##
##
## === CD ===
##
## Call:
## lm(formula = model_formula, data = nil_corrected_combined)
##
## Residuals:
## Min 1Q Median 3Q Max
## -7.6747 -0.7966 0.0991 0.8650 2.9659
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 26.1623 0.5325 49.135 <2e-16 ***
## genotypeInv4m 0.1583 0.7115 0.222 0.825
## Treatment-P -0.1906 0.7010 -0.272 0.787
## genotypeInv4m:Treatment-P -2.6177 0.9768 -2.680 0.010 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.766 on 49 degrees of freedom
## (11 observations deleted due to missingness)
## Multiple R-squared: 0.3109, Adjusted R-squared: 0.2687
## F-statistic: 7.37 on 3 and 49 DF, p-value: 0.0003591
##
##
## === TKW ===
##
## Call:
## lm(formula = model_formula, data = nil_corrected_combined)
##
## Residuals:
## Min 1Q Median 3Q Max
## -46.299 -8.307 2.268 9.737 19.859
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 66.6977 4.0227 16.580 <2e-16 ***
## genotypeInv4m 11.8710 5.3755 2.208 0.0319 *
## Treatment-P -0.3746 5.2961 -0.071 0.9439
## genotypeInv4m:Treatment-P -8.0556 7.3794 -1.092 0.2803
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 13.34 on 49 degrees of freedom
## (11 observations deleted due to missingness)
## Multiple R-squared: 0.1326, Adjusted R-squared: 0.07952
## F-statistic: 2.497 on 3 and 49 DF, p-value: 0.07059
##
##
## === KW50 ===
##
## Call:
## lm(formula = model_formula, data = nil_corrected_combined)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.2755 -0.9140 -0.2116 0.7555 4.7007
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 10.4346 0.4073 25.617 < 2e-16 ***
## genotypeInv4m 0.2029 0.5443 0.373 0.71094
## Treatment-P -1.8014 0.5363 -3.359 0.00152 **
## genotypeInv4m:Treatment-P -0.1211 0.7472 -0.162 0.87193
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.351 on 49 degrees of freedom
## (11 observations deleted due to missingness)
## Multiple R-squared: 0.3437, Adjusted R-squared: 0.3035
## F-statistic: 8.553 on 3 and 49 DF, p-value: 0.0001139
##
##
## === TKN ===
##
## Call:
## lm(formula = model_formula, data = nil_corrected_combined)
##
## Residuals:
## Min 1Q Median 3Q Max
## -255.83 -52.03 11.53 64.12 150.58
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 338.75 26.40 12.833 <2e-16 ***
## genotypeInv4m 39.63 35.27 1.123 0.267
## Treatment-P 69.22 34.75 1.992 0.052 .
## genotypeInv4m:Treatment-P -21.07 48.42 -0.435 0.665
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 87.55 on 49 degrees of freedom
## (11 observations deleted due to missingness)
## Multiple R-squared: 0.1236, Adjusted R-squared: 0.06995
## F-statistic: 2.304 on 3 and 49 DF, p-value: 0.0885
##
##
## === HI ===
##
## Call:
## lm(formula = model_formula, data = nil_corrected_combined)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.19409 -0.03374 0.01116 0.04029 0.12389
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.37796 0.01855 20.381 < 2e-16 ***
## genotypeInv4m 0.07149 0.02478 2.885 0.00581 **
## Treatment-P 0.05926 0.02442 2.427 0.01895 *
## genotypeInv4m:Treatment-P -0.04856 0.03402 -1.427 0.15981
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.06151 on 49 degrees of freedom
## (11 observations deleted due to missingness)
## Multiple R-squared: 0.202, Adjusted R-squared: 0.1532
## F-statistic: 4.135 on 3 and 49 DF, p-value: 0.01085
##
##
## === Ca_grain ===
##
## Call:
## lm(formula = model_formula, data = nil_corrected_combined)
##
## Residuals:
## Min 1Q Median 3Q Max
## -22.1547 -5.0408 0.3704 5.5123 19.0990
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 18.1609 2.6056 6.970 8.18e-09 ***
## genotypeInv4m -0.1958 3.5404 -0.055 0.956
## Treatment-P 18.9580 3.4305 5.526 1.32e-06 ***
## genotypeInv4m:Treatment-P -2.2921 4.8226 -0.475 0.637
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 8.642 on 48 degrees of freedom
## (12 observations deleted due to missingness)
## Multiple R-squared: 0.5388, Adjusted R-squared: 0.51
## F-statistic: 18.69 on 3 and 48 DF, p-value: 3.592e-08
##
##
## === Fe_grain ===
##
## Call:
## lm(formula = model_formula, data = nil_corrected_combined)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.7174 -0.8281 -0.1443 0.8420 3.6400
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 14.5813 0.4654 31.331 <2e-16 ***
## genotypeInv4m 0.4481 0.6323 0.709 0.482
## Treatment-P 0.8421 0.6127 1.374 0.176
## genotypeInv4m:Treatment-P -0.4310 0.8614 -0.500 0.619
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.544 on 48 degrees of freedom
## (12 observations deleted due to missingness)
## Multiple R-squared: 0.04958, Adjusted R-squared: -0.009824
## F-statistic: 0.8346 on 3 and 48 DF, p-value: 0.4815
##
##
## === K_grain ===
##
## Call:
## lm(formula = model_formula, data = nil_corrected_combined)
##
## Residuals:
## Min 1Q Median 3Q Max
## -556.66 -121.01 -27.16 118.69 837.01
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2978.51 81.84 36.392 <2e-16 ***
## genotypeInv4m -143.57 111.20 -1.291 0.203
## Treatment-P 97.92 107.75 0.909 0.368
## genotypeInv4m:Treatment-P -232.26 151.48 -1.533 0.132
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 271.4 on 48 degrees of freedom
## (12 observations deleted due to missingness)
## Multiple R-squared: 0.2383, Adjusted R-squared: 0.1907
## F-statistic: 5.006 on 3 and 48 DF, p-value: 0.004228
##
##
## === Mg_grain ===
##
## Call:
## lm(formula = model_formula, data = nil_corrected_combined)
##
## Residuals:
## Min 1Q Median 3Q Max
## -150.73 -39.67 -8.67 40.35 173.25
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1064.59 23.15 45.982 < 2e-16 ***
## genotypeInv4m -95.55 31.46 -3.037 0.00385 **
## Treatment-P -97.45 30.48 -3.197 0.00246 **
## genotypeInv4m:Treatment-P 25.94 42.85 0.605 0.54781
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 76.79 on 48 degrees of freedom
## (12 observations deleted due to missingness)
## Multiple R-squared: 0.3708, Adjusted R-squared: 0.3315
## F-statistic: 9.43 on 3 and 48 DF, p-value: 5.231e-05
##
##
## === Mn_grain ===
##
## Call:
## lm(formula = model_formula, data = nil_corrected_combined)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.00135 -0.24494 -0.04101 0.15373 1.69538
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.5375 0.1309 34.651 <2e-16 ***
## genotypeInv4m -0.3001 0.1779 -1.687 0.0982 .
## Treatment-P -0.2255 0.1724 -1.308 0.1971
## genotypeInv4m:Treatment-P -0.1563 0.2424 -0.645 0.5219
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4343 on 48 degrees of freedom
## (12 observations deleted due to missingness)
## Multiple R-squared: 0.2466, Adjusted R-squared: 0.1995
## F-statistic: 5.237 on 3 and 48 DF, p-value: 0.003298
##
##
## === P_grain ===
##
## Call:
## lm(formula = model_formula, data = nil_corrected_combined)
##
## Residuals:
## Min 1Q Median 3Q Max
## -431.90 -124.92 -28.28 138.92 710.52
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3050.03 71.15 42.865 < 2e-16 ***
## genotypeInv4m -188.93 96.68 -1.954 0.0565 .
## Treatment-P -671.82 93.68 -7.172 4.02e-09 ***
## genotypeInv4m:Treatment-P -11.80 131.70 -0.090 0.9290
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 236 on 48 degrees of freedom
## (12 observations deleted due to missingness)
## Multiple R-squared: 0.6977, Adjusted R-squared: 0.6788
## F-statistic: 36.92 on 3 and 48 DF, p-value: 1.61e-12
##
##
## === S_grain ===
##
## Call:
## lm(formula = model_formula, data = nil_corrected_combined)
##
## Residuals:
## Min 1Q Median 3Q Max
## -151.336 -47.028 4.235 37.206 191.770
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1088.815 22.736 47.888 <2e-16 ***
## genotypeInv4m 1.035 30.893 0.033 0.9734
## Treatment-P 79.095 29.934 2.642 0.0111 *
## genotypeInv4m:Treatment-P -93.291 42.082 -2.217 0.0314 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 75.41 on 48 degrees of freedom
## (12 observations deleted due to missingness)
## Multiple R-squared: 0.2172, Adjusted R-squared: 0.1683
## F-statistic: 4.441 on 3 and 48 DF, p-value: 0.007821
##
##
## === Zn_grain ===
##
## Call:
## lm(formula = model_formula, data = nil_corrected_combined)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.8328 -1.2633 0.0212 1.0526 6.1201
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 21.1051 0.5734 36.804 <2e-16 ***
## genotypeInv4m -0.2654 0.7792 -0.341 0.735
## Treatment-P 0.1264 0.7550 0.167 0.868
## genotypeInv4m:Treatment-P -0.3019 1.0614 -0.284 0.777
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.902 on 48 degrees of freedom
## (12 observations deleted due to missingness)
## Multiple R-squared: 0.0151, Adjusted R-squared: -0.04646
## F-statistic: 0.2453 on 3 and 48 DF, p-value: 0.8643
##
##
## === Ca_grain_mg ===
##
## Call:
## lm(formula = model_formula, data = nil_corrected_combined)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.66264 -0.40687 0.02316 0.48866 1.17384
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.16155 0.21855 5.315 2.74e-06 ***
## genotypeInv4m 0.14057 0.29696 0.473 0.638
## Treatment-P 1.26432 0.28774 4.394 6.12e-05 ***
## genotypeInv4m:Treatment-P -0.07703 0.40451 -0.190 0.850
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.7249 on 48 degrees of freedom
## (12 observations deleted due to missingness)
## Multiple R-squared: 0.4336, Adjusted R-squared: 0.3982
## F-statistic: 12.25 on 3 and 48 DF, p-value: 4.504e-06
##
##
## === Fe_grain_mg ===
##
## Call:
## lm(formula = model_formula, data = nil_corrected_combined)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.49837 -0.09672 0.00531 0.14628 0.33946
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.93456 0.05453 17.139 < 2e-16 ***
## genotypeInv4m 0.20558 0.07409 2.775 0.00785 **
## Treatment-P 0.08108 0.07179 1.129 0.26433
## genotypeInv4m:Treatment-P -0.14592 0.10093 -1.446 0.15472
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1809 on 48 degrees of freedom
## (12 observations deleted due to missingness)
## Multiple R-squared: 0.1498, Adjusted R-squared: 0.0967
## F-statistic: 2.82 on 3 and 48 DF, p-value: 0.04876
##
##
## === K_grain_mg ===
##
## Call:
## lm(formula = model_formula, data = nil_corrected_combined)
##
## Residuals:
## Min 1Q Median 3Q Max
## -139.075 -19.475 5.616 28.014 102.639
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 195.50 13.51 14.466 <2e-16 ***
## genotypeInv4m 21.04 18.36 1.146 0.258
## Treatment-P 7.55 17.79 0.424 0.673
## genotypeInv4m:Treatment-P -34.81 25.01 -1.392 0.170
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 44.82 on 48 degrees of freedom
## (12 observations deleted due to missingness)
## Multiple R-squared: 0.05223, Adjusted R-squared: -0.007009
## F-statistic: 0.8817 on 3 and 48 DF, p-value: 0.4573
##
##
## === Mg_grain_mg ===
##
## Call:
## lm(formula = model_formula, data = nil_corrected_combined)
##
## Residuals:
## Min 1Q Median 3Q Max
## -40.634 -5.939 0.914 9.759 25.526
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 69.419 3.949 17.577 <2e-16 ***
## genotypeInv4m 4.701 5.366 0.876 0.385
## Treatment-P -6.132 5.200 -1.179 0.244
## genotypeInv4m:Treatment-P -5.404 7.310 -0.739 0.463
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 13.1 on 48 degrees of freedom
## (12 observations deleted due to missingness)
## Multiple R-squared: 0.1256, Adjusted R-squared: 0.07097
## F-statistic: 2.299 on 3 and 48 DF, p-value: 0.08929
##
##
## === Mn_grain_mg ===
##
## Call:
## lm(formula = model_formula, data = nil_corrected_combined)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.175076 -0.025671 0.005726 0.038737 0.083618
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.29010 0.01736 16.707 <2e-16 ***
## genotypeInv4m 0.03828 0.02359 1.622 0.111
## Treatment-P -0.01076 0.02286 -0.471 0.640
## genotypeInv4m:Treatment-P -0.04991 0.03214 -1.553 0.127
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.05759 on 48 degrees of freedom
## (12 observations deleted due to missingness)
## Multiple R-squared: 0.1462, Adjusted R-squared: 0.09285
## F-statistic: 2.74 on 3 and 48 DF, p-value: 0.05347
##
##
## === P_grain_mg ===
##
## Call:
## lm(formula = model_formula, data = nil_corrected_combined)
##
## Residuals:
## Min 1Q Median 3Q Max
## -118.332 -15.488 4.903 21.948 53.931
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 196.43 11.07 17.744 < 2e-16 ***
## genotypeInv4m 22.75 15.04 1.512 0.13698
## Treatment-P -41.14 14.57 -2.823 0.00691 **
## genotypeInv4m:Treatment-P -26.34 20.49 -1.286 0.20469
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 36.72 on 48 degrees of freedom
## (12 observations deleted due to missingness)
## Multiple R-squared: 0.3962, Adjusted R-squared: 0.3584
## F-statistic: 10.5 on 3 and 48 DF, p-value: 2.007e-05
##
##
## === S_grain_mg ===
##
## Call:
## lm(formula = model_formula, data = nil_corrected_combined)
##
## Residuals:
## Min 1Q Median 3Q Max
## -47.963 -8.093 2.278 12.340 22.014
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 70.392 4.559 15.440 <2e-16 ***
## genotypeInv4m 13.226 6.194 2.135 0.0379 *
## Treatment-P 6.176 6.002 1.029 0.3086
## genotypeInv4m:Treatment-P -14.342 8.438 -1.700 0.0957 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 15.12 on 48 degrees of freedom
## (12 observations deleted due to missingness)
## Multiple R-squared: 0.08962, Adjusted R-squared: 0.03272
## F-statistic: 1.575 on 3 and 48 DF, p-value: 0.2076
##
##
## === Zn_grain_mg ===
##
## Call:
## lm(formula = model_formula, data = nil_corrected_combined)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.79638 -0.19906 0.00456 0.23822 0.55141
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.37452 0.09013 15.251 <2e-16 ***
## genotypeInv4m 0.21675 0.12246 1.770 0.0831 .
## Treatment-P 0.02343 0.11866 0.198 0.8443
## genotypeInv4m:Treatment-P -0.16307 0.16681 -0.978 0.3332
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2989 on 48 degrees of freedom
## (12 observations deleted due to missingness)
## Multiple R-squared: 0.07776, Adjusted R-squared: 0.02012
## F-statistic: 1.349 on 3 and 48 DF, p-value: 0.2696
##
##
## === Ca_total_mg ===
##
## Call:
## lm(formula = model_formula, data = nil_corrected_combined)
##
## Residuals:
## Min 1Q Median 3Q Max
## -157.377 -23.970 4.685 28.329 121.318
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 456.24 15.88 28.727 < 2e-16 ***
## genotypeInv4m -82.60 21.58 -3.828 0.000374 ***
## Treatment-P -71.96 20.91 -3.442 0.001208 **
## genotypeInv4m:Treatment-P 30.86 29.40 1.050 0.299002
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 52.67 on 48 degrees of freedom
## (12 observations deleted due to missingness)
## Multiple R-squared: 0.4118, Adjusted R-squared: 0.375
## F-statistic: 11.2 on 3 and 48 DF, p-value: 1.091e-05
##
##
## === Fe_total_mg ===
##
## Call:
## lm(formula = model_formula, data = nil_corrected_combined)
##
## Residuals:
## Min 1Q Median 3Q Max
## -10.398 -3.260 -0.813 1.301 24.175
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 13.6021 1.7005 7.999 2.21e-10 ***
## genotypeInv4m -4.5130 2.3105 -1.953 0.0566 .
## Treatment-P -0.3847 2.2388 -0.172 0.8643
## genotypeInv4m:Treatment-P 3.1100 3.1473 0.988 0.3280
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.64 on 48 degrees of freedom
## (12 observations deleted due to missingness)
## Multiple R-squared: 0.09523, Adjusted R-squared: 0.03868
## F-statistic: 1.684 on 3 and 48 DF, p-value: 0.1829
##
##
## === K_total_mg ===
##
## Call:
## lm(formula = model_formula, data = nil_corrected_combined)
##
## Residuals:
## Min 1Q Median 3Q Max
## -334.13 -153.24 5.93 100.85 743.40
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1338.22 63.93 20.934 < 2e-16 ***
## genotypeInv4m -172.30 86.86 -1.984 0.05303 .
## Treatment-P -227.00 84.16 -2.697 0.00962 **
## genotypeInv4m:Treatment-P 91.88 118.32 0.777 0.44124
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 212 on 48 degrees of freedom
## (12 observations deleted due to missingness)
## Multiple R-squared: 0.2175, Adjusted R-squared: 0.1686
## F-statistic: 4.448 on 3 and 48 DF, p-value: 0.007755
##
##
## === Mg_total_mg ===
##
## Call:
## lm(formula = model_formula, data = nil_corrected_combined)
##
## Residuals:
## Min 1Q Median 3Q Max
## -56.79 -15.27 -0.39 14.73 49.51
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 233.001 7.132 32.670 <2e-16 ***
## genotypeInv4m -10.869 9.691 -1.122 0.2676
## Treatment-P -37.812 9.390 -4.027 0.0002 ***
## genotypeInv4m:Treatment-P 2.042 13.200 0.155 0.8777
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 23.65 on 48 degrees of freedom
## (12 observations deleted due to missingness)
## Multiple R-squared: 0.4015, Adjusted R-squared: 0.3641
## F-statistic: 10.73 on 3 and 48 DF, p-value: 1.633e-05
##
##
## === P_total_mg ===
##
## Call:
## lm(formula = model_formula, data = nil_corrected_combined)
##
## Residuals:
## Min 1Q Median 3Q Max
## -161.423 -16.988 3.266 24.689 138.575
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 441.42 14.49 30.473 < 2e-16 ***
## genotypeInv4m 14.48 19.68 0.736 0.465
## Treatment-P -224.23 19.07 -11.758 9.74e-16 ***
## genotypeInv4m:Treatment-P -21.80 26.81 -0.813 0.420
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 48.04 on 48 degrees of freedom
## (12 observations deleted due to missingness)
## Multiple R-squared: 0.8663, Adjusted R-squared: 0.858
## F-statistic: 103.7 on 3 and 48 DF, p-value: < 2.2e-16
##
##
## === S_total_mg ===
##
## Call:
## lm(formula = model_formula, data = nil_corrected_combined)
##
## Residuals:
## Min 1Q Median 3Q Max
## -55.966 -9.698 -0.140 10.543 39.551
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 166.663 6.044 27.575 <2e-16 ***
## genotypeInv4m 8.486 8.212 1.033 0.3066
## Treatment-P -1.895 7.957 -0.238 0.8127
## genotypeInv4m:Treatment-P -19.388 11.186 -1.733 0.0895 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 20.05 on 48 degrees of freedom
## (12 observations deleted due to missingness)
## Multiple R-squared: 0.1339, Adjusted R-squared: 0.07977
## F-statistic: 2.474 on 3 and 48 DF, p-value: 0.07284
##
##
## === Zn_total_mg ===
##
## Call:
## lm(formula = model_formula, data = nil_corrected_combined)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.38126 -0.28524 -0.00275 0.40330 1.12508
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.62291 0.16188 22.381 <2e-16 ***
## genotypeInv4m 0.03828 0.21995 0.174 0.863
## Treatment-P 0.30537 0.21312 1.433 0.158
## genotypeInv4m:Treatment-P -0.33378 0.29961 -1.114 0.271
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.5369 on 48 degrees of freedom
## (12 observations deleted due to missingness)
## Multiple R-squared: 0.06093, Adjusted R-squared: 0.00224
## F-statistic: 1.038 on 3 and 48 DF, p-value: 0.3842
##
##
## === Ca_HI ===
##
## Call:
## lm(formula = model_formula, data = nil_corrected_combined)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.005313 -0.001034 0.000071 0.001365 0.005104
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.0026785 0.0005876 4.558 3.56e-05 ***
## genotypeInv4m 0.0006395 0.0007984 0.801 0.427
## Treatment-P 0.0037292 0.0007736 4.821 1.48e-05 ***
## genotypeInv4m:Treatment-P 0.0003748 0.0010875 0.345 0.732
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.001949 on 48 degrees of freedom
## (12 observations deleted due to missingness)
## Multiple R-squared: 0.525, Adjusted R-squared: 0.4954
## F-statistic: 17.69 on 3 and 48 DF, p-value: 7.191e-08
##
##
## === Fe_HI ===
##
## Call:
## lm(formula = model_formula, data = nil_corrected_combined)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.076353 -0.020799 0.005299 0.023567 0.051187
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.094503 0.009593 9.852 4.12e-13 ***
## genotypeInv4m 0.040262 0.013034 3.089 0.00333 **
## Treatment-P -0.002256 0.012629 -0.179 0.85900
## genotypeInv4m:Treatment-P -0.020068 0.017754 -1.130 0.26397
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.03181 on 48 degrees of freedom
## (12 observations deleted due to missingness)
## Multiple R-squared: 0.2393, Adjusted R-squared: 0.1918
## F-statistic: 5.034 on 3 and 48 DF, p-value: 0.0041
##
##
## === K_HI ===
##
## Call:
## lm(formula = model_formula, data = nil_corrected_combined)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.102111 -0.020879 0.005362 0.022461 0.074217
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.15334 0.01109 13.821 <2e-16 ***
## genotypeInv4m 0.03574 0.01507 2.371 0.0218 *
## Treatment-P 0.03447 0.01461 2.360 0.0224 *
## genotypeInv4m:Treatment-P -0.03345 0.02053 -1.629 0.1098
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.0368 on 48 degrees of freedom
## (12 observations deleted due to missingness)
## Multiple R-squared: 0.1451, Adjusted R-squared: 0.09162
## F-statistic: 2.715 on 3 and 48 DF, p-value: 0.05507
##
##
## === Mg_HI ===
##
## Call:
## lm(formula = model_formula, data = nil_corrected_combined)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.176810 -0.022169 0.005663 0.026794 0.091461
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.29561 0.01605 18.420 <2e-16 ***
## genotypeInv4m 0.04456 0.02181 2.043 0.0465 *
## Treatment-P 0.03196 0.02113 1.513 0.1369
## genotypeInv4m:Treatment-P -0.03401 0.02970 -1.145 0.2580
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.05323 on 48 degrees of freedom
## (12 observations deleted due to missingness)
## Multiple R-squared: 0.09755, Adjusted R-squared: 0.04114
## F-statistic: 1.729 on 3 and 48 DF, p-value: 0.1735
##
##
## === P_HI ===
##
## Call:
## lm(formula = model_formula, data = nil_corrected_combined)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.23753 -0.03165 0.01261 0.05037 0.13722
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.44406 0.02321 19.133 < 2e-16 ***
## genotypeInv4m 0.03600 0.03154 1.141 0.259
## Treatment-P 0.27995 0.03056 9.162 4.1e-12 ***
## genotypeInv4m:Treatment-P -0.04176 0.04296 -0.972 0.336
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.07698 on 48 degrees of freedom
## (12 observations deleted due to missingness)
## Multiple R-squared: 0.7529, Adjusted R-squared: 0.7375
## F-statistic: 48.76 on 3 and 48 DF, p-value: 1.314e-14
##
##
## === S_HI ===
##
## Call:
## lm(formula = model_formula, data = nil_corrected_combined)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.22205 -0.02288 0.01229 0.04099 0.10331
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.41515 0.02088 19.885 <2e-16 ***
## genotypeInv4m 0.05935 0.02837 2.092 0.0417 *
## Treatment-P 0.04777 0.02749 1.738 0.0886 .
## genotypeInv4m:Treatment-P -0.03614 0.03864 -0.935 0.3543
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.06924 on 48 degrees of freedom
## (12 observations deleted due to missingness)
## Multiple R-squared: 0.1279, Adjusted R-squared: 0.07337
## F-statistic: 2.346 on 3 and 48 DF, p-value: 0.08451
##
##
## === Zn_HI ===
##
## Call:
## lm(formula = model_formula, data = nil_corrected_combined)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.185848 -0.022055 0.009532 0.034694 0.088363
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.37367 0.01829 20.431 <2e-16 ***
## genotypeInv4m 0.06406 0.02485 2.578 0.0131 *
## Treatment-P -0.01651 0.02408 -0.685 0.4963
## genotypeInv4m:Treatment-P -0.01719 0.03385 -0.508 0.6139
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.06066 on 48 degrees of freedom
## (12 observations deleted due to missingness)
## Multiple R-squared: 0.2239, Adjusted R-squared: 0.1754
## F-statistic: 4.616 on 3 and 48 DF, p-value: 0.006453
##
##
## === Ca_ratio ===
##
## Call:
## lm(formula = model_formula, data = nil_corrected_combined)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.0055506 -0.0010635 -0.0002012 0.0011311 0.0061766
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.0039226 0.0006634 5.913 3.41e-07 ***
## genotypeInv4m 0.0007933 0.0009014 0.880 0.383
## Treatment-P 0.0041559 0.0008734 4.758 1.83e-05 ***
## genotypeInv4m:Treatment-P -0.0004932 0.0012279 -0.402 0.690
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.0022 on 48 degrees of freedom
## (12 observations deleted due to missingness)
## Multiple R-squared: 0.4594, Adjusted R-squared: 0.4256
## F-statistic: 13.6 on 3 and 48 DF, p-value: 1.511e-06
##
##
## === Fe_ratio ===
##
## Call:
## lm(formula = model_formula, data = nil_corrected_combined)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.150452 -0.033592 0.009758 0.030325 0.124006
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.14877 0.01387 10.729 2.41e-14 ***
## genotypeInv4m 0.05307 0.01884 2.817 0.00702 **
## Treatment-P -0.01874 0.01826 -1.027 0.30977
## genotypeInv4m:Treatment-P -0.03193 0.02567 -1.244 0.21946
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.04599 on 48 degrees of freedom
## (12 observations deleted due to missingness)
## Multiple R-squared: 0.2736, Adjusted R-squared: 0.2282
## F-statistic: 6.025 on 3 and 48 DF, p-value: 0.001439
##
##
## === K_ratio ===
##
## Call:
## lm(formula = model_formula, data = nil_corrected_combined)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.10033 -0.04145 -0.00832 0.02327 0.13151
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.28760 0.01677 17.152 <2e-16 ***
## genotypeInv4m 0.01377 0.02278 0.604 0.549
## Treatment-P 0.01485 0.02208 0.673 0.504
## genotypeInv4m:Treatment-P -0.02976 0.03103 -0.959 0.342
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.05561 on 48 degrees of freedom
## (12 observations deleted due to missingness)
## Multiple R-squared: 0.01923, Adjusted R-squared: -0.04207
## F-statistic: 0.3137 on 3 and 48 DF, p-value: 0.8154
##
##
## === Mg_ratio ===
##
## Call:
## lm(formula = model_formula, data = nil_corrected_combined)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.140810 -0.030710 0.002109 0.031276 0.167363
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.68729 0.01877 36.626 < 2e-16 ***
## genotypeInv4m -0.05176 0.02550 -2.030 0.04790 *
## Treatment-P -0.06732 0.02471 -2.725 0.00895 **
## genotypeInv4m:Treatment-P 0.02646 0.03473 0.762 0.44991
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.06224 on 48 degrees of freedom
## (12 observations deleted due to missingness)
## Multiple R-squared: 0.2253, Adjusted R-squared: 0.1769
## F-statistic: 4.653 on 3 and 48 DF, p-value: 0.006197
##
##
## === P_ratio ===
##
## Call:
## lm(formula = model_formula, data = nil_corrected_combined)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.89967 -0.16127 0.04248 0.15302 1.52704
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.3798 0.1234 11.184 5.74e-15 ***
## genotypeInv4m -0.2311 0.1676 -1.378 0.174
## Treatment-P 2.0921 0.1624 12.879 < 2e-16 ***
## genotypeInv4m:Treatment-P -0.0899 0.2284 -0.394 0.696
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4092 on 48 degrees of freedom
## (12 observations deleted due to missingness)
## Multiple R-squared: 0.8751, Adjusted R-squared: 0.8673
## F-statistic: 112.1 on 3 and 48 DF, p-value: < 2.2e-16
##
##
## === S_ratio ===
##
## Call:
## lm(formula = model_formula, data = nil_corrected_combined)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.290178 -0.042631 -0.006513 0.051470 0.237904
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.19186 0.02762 43.153 <2e-16 ***
## genotypeInv4m -0.05989 0.03753 -1.596 0.1171
## Treatment-P -0.07777 0.03636 -2.139 0.0376 *
## genotypeInv4m:Treatment-P 0.07122 0.05112 1.393 0.1699
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.0916 on 48 degrees of freedom
## (12 observations deleted due to missingness)
## Multiple R-squared: 0.09648, Adjusted R-squared: 0.04001
## F-statistic: 1.709 on 3 and 48 DF, p-value: 0.1778
##
##
## === Zn_ratio ===
##
## Call:
## lm(formula = model_formula, data = nil_corrected_combined)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.23475 -0.06879 -0.00414 0.05575 0.35589
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.01920 0.03757 27.127 < 2e-16 ***
## genotypeInv4m -0.03704 0.05105 -0.725 0.4717
## Treatment-P -0.30571 0.04946 -6.180 1.33e-07 ***
## genotypeInv4m:Treatment-P 0.12814 0.06954 1.843 0.0715 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1246 on 48 degrees of freedom
## (12 observations deleted due to missingness)
## Multiple R-squared: 0.5273, Adjusted R-squared: 0.4977
## F-statistic: 17.85 on 3 and 48 DF, p-value: 6.438e-08
##
##
## === Mn_total_mg ===
##
## Call:
## lm(formula = model_formula, data = nil_corrected_combined)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.1634 -0.3300 -0.0206 0.2955 1.8196
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.8404 0.1929 25.090 <2e-16 ***
## genotypeInv4m -0.3836 0.2671 -1.436 0.158
## Treatment-P -0.4232 0.2540 -1.666 0.102
## genotypeInv4m:Treatment-P 0.1225 0.3607 0.340 0.736
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.6399 on 47 degrees of freedom
## (13 observations deleted due to missingness)
## Multiple R-squared: 0.1277, Adjusted R-squared: 0.07199
## F-statistic: 2.293 on 3 and 47 DF, p-value: 0.09019
##
##
## === Mn_HI ===
##
## Call:
## lm(formula = model_formula, data = nil_corrected_combined)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.033580 -0.006102 0.001129 0.009401 0.027214
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.061723 0.004016 15.370 <2e-16 ***
## genotypeInv4m 0.010494 0.005560 1.888 0.0653 .
## Treatment-P 0.002244 0.005287 0.425 0.6731
## genotypeInv4m:Treatment-P -0.005273 0.007509 -0.702 0.4860
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.01332 on 47 degrees of freedom
## (13 observations deleted due to missingness)
## Multiple R-squared: 0.09054, Adjusted R-squared: 0.03249
## F-statistic: 1.56 on 3 and 47 DF, p-value: 0.2116
##
##
## === Mn_ratio ===
##
## Call:
## lm(formula = model_formula, data = nil_corrected_combined)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.0206469 -0.0094358 -0.0003677 0.0076911 0.0314576
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.103021 0.003890 26.481 < 2e-16 ***
## genotypeInv4m -0.004689 0.005386 -0.871 0.38841
## Treatment-P -0.014929 0.005122 -2.915 0.00544 **
## genotypeInv4m:Treatment-P 0.002685 0.007274 0.369 0.71365
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.0129 on 47 degrees of freedom
## (13 observations deleted due to missingness)
## Multiple R-squared: 0.2366, Adjusted R-squared: 0.1879
## F-statistic: 4.857 on 3 and 47 DF, p-value: 0.005038# Combine all coefficients into single data frame
if (length(model_coefficients) > 0) {
all_coefficients <- bind_rows(model_coefficients) %>%
mutate(
term = case_when(
term == "(Intercept)" ~ "Intercept (CTRL, +P)",
term == "genotypeInv4m" ~ "Genotype Effect (Inv4m vs CTRL)",
term == "Treatment-P" ~ "Treatment Effect (-P vs +P)",
term == "genotypeInv4m:Treatment-P" ~ "Inv4m:-P",
TRUE ~ term
)
)
# Apply FDR correction within each term type
all_coefficients <- all_coefficients %>%
group_by(term) %>%
mutate(
p_adj_fdr = p.adjust(p_value, method = "fdr")
) %>%
ungroup()
cat("\n=== All Model Coefficients ===\n")
print(all_coefficients %>% filter(grepl("Interaction",term)) %>% arrange(p_value))
}##
## === All Model Coefficients ===
## # A tibble: 0 × 7
## # ℹ 7 variables: trait <chr>, term <chr>, estimate <dbl>, std_error <dbl>,
## # t_value <dbl>, p_value <dbl>, p_adj_fdr <dbl>9. Effect tables: Treatment effects and Genotype × Treatment interactions
if (exists("all_coefficients") && nrow(all_coefficients) > 0) {
# Treatment effects table (main effect of -P vs +P)
treatment_effects <- all_coefficients %>%
filter(term == "Treatment Effect (-P vs +P)") %>%
arrange(p_value) %>%
mutate(
significance = case_when(
p_adj_fdr < 0.001 ~ "***",
p_adj_fdr < 0.01 ~ "**",
p_adj_fdr < 0.05 ~ "*",
p_adj_fdr < 0.1 ~ ".",
TRUE ~ ""
)
)
cat("\n=== Treatment Effects: -P vs +P (Reference) ===\n")
cat("Negative estimates indicate reduction under -P treatment\n")
cat("Positive estimates indicate increase under -P treatment\n\n")
print(treatment_effects %>%
select(trait, estimate, std_error, t_value, p_value, p_adj_fdr, significance) %>%
mutate(across(where(is.numeric), ~ round(.x, 4))))
# Genotype × Treatment interaction effects table
interaction_effects <- all_coefficients %>%
filter(term == "Inv4m:-P") %>%
arrange(p_value) %>%
mutate(
significance = case_when(
p_adj_fdr < 0.001 ~ "***",
p_adj_fdr < 0.01 ~ "**",
p_adj_fdr < 0.05 ~ "*",
p_adj_fdr < 0.1 ~ ".",
TRUE ~ ""
)
)
cat("\n=== Inv4m:-P Effects ===\n")
cat("Tests whether Inv4m responds differently to -P vs +P compared to CTRL\n")
cat("Negative estimates: Inv4m shows greater reduction under -P than CTRL\n")
cat("Positive estimates: Inv4m shows less reduction (or greater increase) under -P than CTRL\n\n")
print(interaction_effects %>%
select(trait, estimate, std_error, t_value, p_value, p_adj_fdr, significance) %>%
mutate(across(where(is.numeric), ~ round(.x, 4))))
# Genotype main effects table
genotype_effects <- all_coefficients %>%
filter(term == "Genotype Effect (Inv4m vs CTRL)") %>%
arrange(p_value) %>%
mutate(
significance = case_when(
p_adj_fdr < 0.001 ~ "***",
p_adj_fdr < 0.01 ~ "**",
p_adj_fdr < 0.05 ~ "*",
p_adj_fdr < 0.1 ~ ".",
TRUE ~ ""
)
)
cat("\n=== Genotype Effects: Inv4m vs CTRL (at +P reference level) ===\n")
print(genotype_effects %>%
select(trait, estimate, std_error, t_value, p_value, p_adj_fdr, significance) %>%
mutate(across(where(is.numeric), ~ round(.x, 4))))
# Summary counts
cat("\n=== Summary of Significant Effects (FDR < 0.05) ===\n")
cat("Treatment effects:", sum(treatment_effects$p_adj_fdr < 0.05, na.rm = TRUE), "of", nrow(treatment_effects), "\n")
cat("Interaction effects:", sum(interaction_effects$p_adj_fdr < 0.05, na.rm = TRUE), "of", nrow(interaction_effects), "\n")
cat("Genotype effects:", sum(genotype_effects$p_adj_fdr < 0.05, na.rm = TRUE), "of", nrow(genotype_effects), "\n")
}##
## === Treatment Effects: -P vs +P (Reference) ===
## Negative estimates indicate reduction under -P treatment
## Positive estimates indicate increase under -P treatment
##
## # A tibble: 72 × 7
## trait estimate std_error t_value p_value p_adj_fdr significance
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <chr>
## 1 P_stover -1592. 84.6 -18.8 0 0 ***
## 2 P_ratio 2.09 0.162 12.9 0 0 ***
## 3 P_stover_mg -186. 15.7 -11.8 0 0 ***
## 4 STW50 -18.8 1.58 -11.9 0 0 ***
## 5 P_total_mg -224. 19.1 -11.8 0 0 ***
## 6 DTS 3.06 0.322 9.51 0 0 ***
## 7 DTA 3.43 0.368 9.34 0 0 ***
## 8 P_HI 0.28 0.0306 9.16 0 0 ***
## 9 P_grain -672. 93.7 -7.17 0 0 ***
## 10 STW60 -33.3 4.97 -6.70 0 0 ***
## # ℹ 62 more rows
##
## === Inv4m:-P Effects ===
## Tests whether Inv4m responds differently to -P vs +P compared to CTRL
## Negative estimates: Inv4m shows greater reduction under -P than CTRL
## Positive estimates: Inv4m shows less reduction (or greater increase) under -P than CTRL
##
## # A tibble: 72 × 7
## trait estimate std_error t_value p_value p_adj_fdr significance
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <chr>
## 1 S_stover -93.8 29.2 -3.21 0.0022 0.156 ""
## 2 CD -2.62 0.977 -2.68 0.01 0.36 ""
## 3 S_grain -93.3 42.1 -2.22 0.0314 0.627 ""
## 4 STW50 4.62 2.22 2.08 0.0425 0.627 ""
## 5 STW40 2.86 1.41 2.02 0.0474 0.627 ""
## 6 EL -1.06 0.533 -1.99 0.0522 0.627 ""
## 7 Zn_ratio 0.128 0.0695 1.84 0.0715 0.672 ""
## 8 S_total_mg -19.4 11.2 -1.73 0.0895 0.672 ""
## 9 S_grain_mg -14.3 8.44 -1.70 0.0957 0.672 ""
## 10 K_HI -0.0335 0.0205 -1.63 0.110 0.672 ""
## # ℹ 62 more rows
##
## === Genotype Effects: Inv4m vs CTRL (at +P reference level) ===
## # A tibble: 72 × 7
## trait estimate std_error t_value p_value p_adj_fdr significance
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <chr>
## 1 Ca_stover_mg -93.6 20.6 -4.55 0 0.002 **
## 2 DTA -1.48 0.368 -4.02 0.0002 0.0044 **
## 3 DTS -1.28 0.322 -3.99 0.0002 0.0044 **
## 4 Ca_total_mg -82.6 21.6 -3.83 0.0004 0.0055 **
## 5 PH 5.57 1.48 3.76 0.0004 0.0055 **
## 6 Fe_HI 0.0403 0.013 3.09 0.0033 0.0396 *
## 7 Mg_grain -95.6 31.5 -3.04 0.0039 0.0396 *
## 8 Ca_stover -411. 141. -2.92 0.005 0.0405 *
## 9 STW40 -2.93 1.01 -2.91 0.0051 0.0405 *
## 10 HI 0.0715 0.0248 2.88 0.0058 0.0418 *
## # ℹ 62 more rows
##
## === Summary of Significant Effects (FDR < 0.05) ===
## Treatment effects: 34 of 72
## Interaction effects: 0 of 72
## Genotype effects: 13 of 7210. Visualization: Treatment and Genotype × Treatment effects
# Define plotting theme
pheno_theme2 <- ggpubr::theme_classic2(base_size = 16) +
ggplot2::theme(
plot.title = ggtext::element_markdown(hjust = 0.5, face = "bold"),
axis.title.y = ggtext::element_markdown(),
axis.title.x = element_text(face = "bold"),
axis.text.x = element_text(size = 14, face = "bold", color = "black"),
strip.background = element_blank(),
strip.text = ggtext::element_markdown(size = 14),
legend.position = "bottom",
legend.title = element_text(face = "bold"),
legend.text = element_text(size = 12)
)
if (length(common_phenos) > 0 && exists("all_coefficients")) {
# Prepare data for plotting
plot_data_combined <- nil_corrected_combined %>%
pivot_longer(
cols = all_of(common_phenos),
names_to = "trait",
values_to = "value"
) %>%
filter(!is.na(value)) %>%
mutate(trait = as.character(trait))
# Treatment × Genotype interaction plot (emphasizing treatment effect first)
interaction_plot <- ggplot(
plot_data_combined,
aes(x = Treatment, y = value, color = genotype, group = genotype)
) +
stat_summary(
fun = mean,
geom = "point",
size = 4
) +
stat_summary(
fun = mean,
geom = "line",
linewidth = 1.5
) +
stat_summary(
fun.data = mean_se,
geom = "errorbar",
width = 0.1,
linewidth = 1
) +
scale_color_manual(values = c("CTRL" = "gold", "Inv4m" = "purple4")) +
scale_x_discrete(labels = c("+P" = "+P\n(Reference)", "-P" = "-P\n(Low P)")) +
labs(
title = "Treatment Effects and Inv4m:-Ps",
subtitle = "+P is reference treatment; Lines show genotype-specific responses",
y = "Spatially Corrected Value",
x = "Phosphorus Treatment",
color = "Genotype"
) +
facet_wrap(~ trait, scales = "free_y", ncol = 4) +
pheno_theme2
print(interaction_plot)
# Focus on significant treatment effects
if (exists("treatment_effects")) {
sig_treatment_traits <- treatment_effects %>%
filter(p_adj_fdr < 0.05) %>%
pull(trait)
if (length(sig_treatment_traits) > 0) {
treatment_plot_data <- plot_data_combined %>%
filter(trait %in% sig_treatment_traits)
treatment_plot <- ggplot(
treatment_plot_data,
aes(x = Treatment, y = value, fill = Treatment)
) +
geom_boxplot(width = 0.6, alpha = 0.7, linewidth = 0.8) +
geom_point(
position = position_jitter(width = 0.2),
size = 2,
alpha = 0.6
) +
scale_fill_manual(
values = c("+P" = "lightblue", "-P" = "lightcoral"),
labels = c("+P" = "+P (Reference)", "-P" = "-P (Low P)")
) +
labs(
title = "Significant Treatment Effects (FDR < 0.05)",
subtitle = "Traits showing significant response to phosphorus limitation",
y = "Spatially Corrected Value",
x = "Phosphorus Treatment",
fill = "Treatment"
) +
facet_wrap(~ trait, scales = "free_y", ncol = 3) +
pheno_theme2
print(treatment_plot)
}
}
# Focus on significant interaction effects
if (exists("interaction_effects")) {
sig_interaction_traits <- interaction_effects %>%
filter(p_adj_fdr < 0.05) %>%
pull(trait)
if (length(sig_interaction_traits) > 0) {
interaction_plot_data <- plot_data_combined %>%
filter(trait %in% sig_interaction_traits)
interaction_detailed_plot <- ggplot(
interaction_plot_data,
aes(x = genotype, y = value, fill = Treatment)
) +
geom_boxplot(
width = 0.6,
alpha = 0.7,
position = position_dodge(width = 0.8),
linewidth = 0.8
) +
geom_point(
aes(color = Treatment),
position = position_jitterdodge(
dodge.width = 0.8,
jitter.width = 0.2
),
size = 1.5,
alpha = 0.6
) +
scale_fill_manual(
values = c("+P" = "lightblue", "-P" = "lightcoral"),
labels = c("+P" = "+P (Reference)", "-P" = "-P (Low P)")
) +
scale_color_manual(
values = c("+P" = "blue", "-P" = "red"),
labels = c("+P" = "+P (Reference)", "-P" = "-P (Low P)")
) +
labs(
title = "Significant Inv4m:-Ps (FDR < 0.05)",
subtitle = "Genotypes respond differently to phosphorus limitation",
y = "Spatially Corrected Value",
x = "Genotype",
fill = "P Treatment",
color = "P Treatment"
) +
facet_wrap(~ trait, scales = "free_y", ncol = 3) +
pheno_theme2 +
guides(color = "none") # Remove duplicate legend
print(interaction_detailed_plot)
}
}
}
11. Summary statistics and effect sizes
# Calculate effect sizes and summary statistics
effect_summary <- nil_corrected_combined %>%
pivot_longer(
cols = all_of(common_phenos),
names_to = "trait",
values_to = "value"
) %>%
filter(!is.na(value)) %>%
group_by(trait, Treatment, genotype) %>%
summarise(
n = n(),
mean = mean(value),
sd = sd(value),
se = sd / sqrt(n),
.groups = "drop"
) %>%
pivot_wider(
names_from = genotype,
values_from = c(n, mean, sd, se),
names_sep = "_"
) %>%
mutate(
diff_mean = mean_Inv4m - mean_CTRL,
percent_change = (diff_mean / mean_CTRL) * 100,
pooled_sd = sqrt(((n_Inv4m - 1) * sd_Inv4m^2 + (n_CTRL - 1) * sd_CTRL^2) /
(n_Inv4m + n_CTRL - 2)),
cohens_d = diff_mean / pooled_sd
) %>%
arrange(Treatment, desc(abs(cohens_d)))
cat("\n=== Effect Size Summary ===\n")##
## === Effect Size Summary ===print(effect_summary %>%
select(trait, Treatment, diff_mean, percent_change, cohens_d) %>%
mutate(across(where(is.numeric), ~ round(.x, 3))))## # A tibble: 144 × 5
## trait Treatment diff_mean percent_change cohens_d
## <chr> <fct> <dbl> <dbl> <dbl>
## 1 DTS +P -1.28 -1.74 -1.56
## 2 P_ratio +P -0.231 -16.7 -1.47
## 3 DTA +P -1.48 -2.05 -1.43
## 4 Ca_total_mg +P -82.6 -18.1 -1.43
## 5 Ca_stover_mg +P -93.6 -19.3 -1.42
## 6 PH +P 5.57 3.44 1.29
## 7 Fe_HI +P 0.04 42.6 1.24
## 8 Mg_grain +P -95.6 -8.98 -1.05
## 9 Fe_grain_mg +P 0.206 22.0 0.991
## 10 Ca_stover +P -411. -9.07 -0.974
## # ℹ 134 more rows12. Export results
# Export the combined spatially corrected dataset
combined_output_file <- file.path(paths$intermediate,
"PSU2022_spatially_corrected_both_treatments_with_ionome.csv")
write.csv(sp_corrected_combined, combined_output_file, row.names = FALSE)
# Export all model coefficients
if (exists("all_coefficients")) {
coefficients_output_file <- file.path(paths$intermediate, "PSU2022_model_coefficients_spatially_corrected.csv")
write.csv(all_coefficients, coefficients_output_file, row.names = FALSE)
# Export treatment effects table
if (exists("treatment_effects")) {
treatment_output_file <- file.path(paths$intermediate, "PSU2022_treatment_effects_spatially_corrected.csv")
write.csv(treatment_effects, treatment_output_file, row.names = FALSE)
}
# Export interaction effects table
if (exists("interaction_effects")) {
interaction_output_file <- file.path(paths$intermediate, "PSU2022_interaction_effects_spatially_corrected.csv")
write.csv(interaction_effects, interaction_output_file, row.names = FALSE)
}
# Export genotype effects table
if (exists("genotype_effects")) {
genotype_output_file <- file.path(paths$intermediate, "PSU2022_genotype_effects_spatially_corrected.csv")
write.csv(genotype_effects, genotype_output_file, row.names = FALSE)
}
}
# Export summary statistics with treatment contrasts
effect_summary_treatment <- nil_corrected_combined %>%
pivot_longer(
cols = all_of(common_phenos),
names_to = "trait",
values_to = "value"
) %>%
filter(!is.na(value)) %>%
group_by(trait, Treatment, genotype) %>%
summarise(
n = n(),
mean = mean(value),
sd = sd(value),
se = sd / sqrt(n),
.groups = "drop"
) %>%
# Calculate treatment effects (relative to +P reference)
pivot_wider(
names_from = Treatment,
values_from = c(n, mean, sd, se),
names_sep = "_"
) %>%
mutate(
treatment_effect = `mean_-P` - `mean_+P`,
treatment_percent_change = (treatment_effect / `mean_+P`) * 100
) %>%
# Add genotype comparison within each treatment
group_by(trait) %>%
mutate(
genotype_effect_plus_p = case_when(
genotype == "Inv4m" ~ `mean_+P` - `mean_+P`[genotype == "CTRL"],
TRUE ~ NA_real_
),
genotype_effect_minus_p = case_when(
genotype == "Inv4m" ~ `mean_-P` - `mean_-P`[genotype == "CTRL"],
TRUE ~ NA_real_
)
) %>%
ungroup()
summary_output_file <- file.path(paths$intermediate, "PSU2022_effect_summary_with_treatment_reference.csv")
write.csv(effect_summary_treatment, summary_output_file, row.names = FALSE)
cat("=== Files Exported ===\n")## === Files Exported ===cat("Spatially corrected data:", combined_output_file, "\n")## Spatially corrected data: /Users/fvrodriguez/Library/CloudStorage/GoogleDrive-frodrig4@ncsu.edu/My Drive/repos/inv4m/results/phosphorus_paper/intermediate/PSU2022_spatially_corrected_both_treatments_with_ionome.csvif (exists("all_coefficients")) {
cat("All model coefficients:", coefficients_output_file, "\n")
if (exists("treatment_effects")) cat("Treatment effects:", treatment_output_file, "\n")
if (exists("interaction_effects")) cat("Interaction effects:", interaction_output_file, "\n")
if (exists("genotype_effects")) cat("Genotype effects:", genotype_output_file, "\n")
}## All model coefficients: /Users/fvrodriguez/Library/CloudStorage/GoogleDrive-frodrig4@ncsu.edu/My Drive/repos/inv4m/results/phosphorus_paper/intermediate/PSU2022_model_coefficients_spatially_corrected.csv
## Treatment effects: /Users/fvrodriguez/Library/CloudStorage/GoogleDrive-frodrig4@ncsu.edu/My Drive/repos/inv4m/results/phosphorus_paper/intermediate/PSU2022_treatment_effects_spatially_corrected.csv
## Interaction effects: /Users/fvrodriguez/Library/CloudStorage/GoogleDrive-frodrig4@ncsu.edu/My Drive/repos/inv4m/results/phosphorus_paper/intermediate/PSU2022_interaction_effects_spatially_corrected.csv
## Genotype effects: /Users/fvrodriguez/Library/CloudStorage/GoogleDrive-frodrig4@ncsu.edu/My Drive/repos/inv4m/results/phosphorus_paper/intermediate/PSU2022_genotype_effects_spatially_corrected.csvcat("Effect summary (treatment reference):", summary_output_file, "\n")## Effect summary (treatment reference): /Users/fvrodriguez/Library/CloudStorage/GoogleDrive-frodrig4@ncsu.edu/My Drive/repos/inv4m/results/phosphorus_paper/intermediate/PSU2022_effect_summary_with_treatment_reference.csv# Final analysis summary
cat("\n=== Final Analysis Summary ===\n")##
## === Final Analysis Summary ===cat("Total phenotypes evaluated:", length(all_phenotypes), "\n")## Total phenotypes evaluated: 72cat("Usable phenotypes (<50% missing):", length(usable_phenotypes), "\n")## Usable phenotypes (<50% missing): 72cat("Successfully fitted models (+P):", length(all_models_plus_p), "\n")## Successfully fitted models (+P): 72cat("Successfully fitted models (-P):", length(all_models_minus_p), "\n")## Successfully fitted models (-P): 72cat("Common phenotypes between treatments:", length(common_phenos), "\n")## Common phenotypes between treatments: 72cat("Total plots with corrected data:", nrow(sp_corrected_combined), "\n")## Total plots with corrected data: 127cat("Significant traits (by treatment):", length(sig_treatment_traits), "\n")## Significant traits (by treatment): 34# Print model fit information for both treatments
if (length(fit_info_plus_p) > 0) {
cat("\n=== +P Model Fit Details ===\n")
for (trait in names(fit_info_plus_p)) {
info <- fit_info_plus_p[[trait]]
cat(sprintf("%s: n=%d, status=%s\n", trait, info$n_obs, info$status))
}
}##
## === +P Model Fit Details ===
## PH: n=64, status=success
## DTA: n=64, status=success
## DTS: n=64, status=success
## STW40: n=62, status=success
## STWHV: n=62, status=success
## Ca_stover: n=62, status=success
## K_stover: n=62, status=success
## Mg_stover: n=62, status=success
## P_stover: n=62, status=success
## S_stover: n=62, status=success
## Zn_stover: n=62, status=success
## Ca_stover_mg: n=62, status=success
## K_stover_mg: n=62, status=success
## Mg_stover_mg: n=62, status=success
## P_stover_mg: n=62, status=success
## S_stover_mg: n=62, status=success
## Zn_stover_mg: n=62, status=success
## Fe_stover: n=61, status=success
## Fe_stover_mg: n=61, status=success
## Mn_stover: n=60, status=success
## Mn_stover_mg: n=60, status=success
## STW50: n=51, status=success
## STW60: n=50, status=success
## EW: n=49, status=success
## EL: n=49, status=success
## ED: n=49, status=success
## KRN: n=49, status=success
## CD: n=49, status=success
## TKW: n=49, status=success
## KW50: n=49, status=success
## TKN: n=49, status=success
## HI: n=49, status=success
## Ca_grain: n=48, status=success
## Fe_grain: n=48, status=success
## K_grain: n=48, status=success
## Mg_grain: n=48, status=success
## Mn_grain: n=48, status=success
## P_grain: n=48, status=success
## S_grain: n=48, status=success
## Zn_grain: n=48, status=success
## Ca_grain_mg: n=48, status=success
## Fe_grain_mg: n=48, status=success
## K_grain_mg: n=48, status=success
## Mg_grain_mg: n=48, status=success
## Mn_grain_mg: n=48, status=success
## P_grain_mg: n=48, status=success
## S_grain_mg: n=48, status=success
## Zn_grain_mg: n=48, status=success
## Ca_total_mg: n=48, status=success
## Fe_total_mg: n=48, status=success
## K_total_mg: n=48, status=success
## Mg_total_mg: n=48, status=success
## P_total_mg: n=48, status=success
## S_total_mg: n=48, status=success
## Zn_total_mg: n=48, status=success
## Ca_HI: n=48, status=success
## Fe_HI: n=48, status=success
## K_HI: n=48, status=success
## Mg_HI: n=48, status=success
## P_HI: n=48, status=success
## S_HI: n=48, status=success
## Zn_HI: n=48, status=success
## Ca_ratio: n=48, status=success
## Fe_ratio: n=48, status=success
## K_ratio: n=48, status=success
## Mg_ratio: n=48, status=success
## P_ratio: n=48, status=success
## S_ratio: n=48, status=success
## Zn_ratio: n=48, status=success
## Mn_total_mg: n=47, status=success
## Mn_HI: n=47, status=success
## Mn_ratio: n=47, status=success# Final analysis summary
cat("\n=== Final Analysis Summary ===\n")##
## === Final Analysis Summary ===cat("Total phenotypes evaluated:", length(all_phenotypes), "\n")## Total phenotypes evaluated: 72cat("Usable phenotypes (<50% missing):", length(usable_phenotypes), "\n")## Usable phenotypes (<50% missing): 72cat("Successfully fitted models (+P):", length(all_models_plus_p), "\n")## Successfully fitted models (+P): 72cat("Successfully fitted models (-P):", length(all_models_minus_p), "\n")## Successfully fitted models (-P): 72cat("Common phenotypes between treatments:", length(common_phenos), "\n")## Common phenotypes between treatments: 72cat("Total plots with corrected data:", nrow(sp_corrected_combined), "\n")## Total plots with corrected data: 127if (exists("treatment_effects") && exists("interaction_effects") && exists("genotype_effects")) {
cat("Significant treatment effects (FDR < 0.05):", sum(treatment_effects$p_adj_fdr < 0.05, na.rm = TRUE), "of", nrow(treatment_effects), "\n")
cat("Significant interaction effects (FDR < 0.05):", sum(interaction_effects$p_adj_fdr < 0.05, na.rm = TRUE), "of", nrow(interaction_effects), "\n")
cat("Significant genotype effects (FDR < 0.05):", sum(genotype_effects$p_adj_fdr < 0.05, na.rm = TRUE), "of", nrow(genotype_effects), "\n")
}## Significant treatment effects (FDR < 0.05): 34 of 72
## Significant interaction effects (FDR < 0.05): 0 of 72
## Significant genotype effects (FDR < 0.05): 13 of 72# Print model fit information for both treatments
if (length(fit_info_plus_p) > 0) {
cat("\n=== +P Model Fit Details ===\n")
for (trait in names(fit_info_plus_p)) {
info <- fit_info_plus_p[[trait]]
cat(sprintf("%s: n=%d, status=%s\n", trait, info$n_obs, info$status))
}
}##
## === +P Model Fit Details ===
## PH: n=64, status=success
## DTA: n=64, status=success
## DTS: n=64, status=success
## STW40: n=62, status=success
## STWHV: n=62, status=success
## Ca_stover: n=62, status=success
## K_stover: n=62, status=success
## Mg_stover: n=62, status=success
## P_stover: n=62, status=success
## S_stover: n=62, status=success
## Zn_stover: n=62, status=success
## Ca_stover_mg: n=62, status=success
## K_stover_mg: n=62, status=success
## Mg_stover_mg: n=62, status=success
## P_stover_mg: n=62, status=success
## S_stover_mg: n=62, status=success
## Zn_stover_mg: n=62, status=success
## Fe_stover: n=61, status=success
## Fe_stover_mg: n=61, status=success
## Mn_stover: n=60, status=success
## Mn_stover_mg: n=60, status=success
## STW50: n=51, status=success
## STW60: n=50, status=success
## EW: n=49, status=success
## EL: n=49, status=success
## ED: n=49, status=success
## KRN: n=49, status=success
## CD: n=49, status=success
## TKW: n=49, status=success
## KW50: n=49, status=success
## TKN: n=49, status=success
## HI: n=49, status=success
## Ca_grain: n=48, status=success
## Fe_grain: n=48, status=success
## K_grain: n=48, status=success
## Mg_grain: n=48, status=success
## Mn_grain: n=48, status=success
## P_grain: n=48, status=success
## S_grain: n=48, status=success
## Zn_grain: n=48, status=success
## Ca_grain_mg: n=48, status=success
## Fe_grain_mg: n=48, status=success
## K_grain_mg: n=48, status=success
## Mg_grain_mg: n=48, status=success
## Mn_grain_mg: n=48, status=success
## P_grain_mg: n=48, status=success
## S_grain_mg: n=48, status=success
## Zn_grain_mg: n=48, status=success
## Ca_total_mg: n=48, status=success
## Fe_total_mg: n=48, status=success
## K_total_mg: n=48, status=success
## Mg_total_mg: n=48, status=success
## P_total_mg: n=48, status=success
## S_total_mg: n=48, status=success
## Zn_total_mg: n=48, status=success
## Ca_HI: n=48, status=success
## Fe_HI: n=48, status=success
## K_HI: n=48, status=success
## Mg_HI: n=48, status=success
## P_HI: n=48, status=success
## S_HI: n=48, status=success
## Zn_HI: n=48, status=success
## Ca_ratio: n=48, status=success
## Fe_ratio: n=48, status=success
## K_ratio: n=48, status=success
## Mg_ratio: n=48, status=success
## P_ratio: n=48, status=success
## S_ratio: n=48, status=success
## Zn_ratio: n=48, status=success
## Mn_total_mg: n=47, status=success
## Mn_HI: n=47, status=success
## Mn_ratio: n=47, status=successif (length(fit_info_minus_p) > 0) {
cat("\n=== -P Model Fit Details ===\n")
for (trait in names(fit_info_minus_p)) {
info <- fit_info_minus_p[[trait]]
cat(sprintf("%s: n=%d, status=%s\n", trait, info$n_obs, info$status))
}
}##
## === -P Model Fit Details ===
## PH: n=63, status=success
## DTA: n=63, status=success
## DTS: n=63, status=success
## STW40: n=63, status=success
## STWHV: n=63, status=success
## Ca_stover: n=63, status=success
## K_stover: n=63, status=success
## Mg_stover: n=63, status=success
## P_stover: n=63, status=success
## S_stover: n=63, status=success
## Zn_stover: n=63, status=success
## Ca_stover_mg: n=63, status=success
## K_stover_mg: n=63, status=success
## Mg_stover_mg: n=63, status=success
## P_stover_mg: n=63, status=success
## S_stover_mg: n=63, status=success
## Zn_stover_mg: n=63, status=success
## Fe_stover: n=63, status=success
## Fe_stover_mg: n=63, status=success
## Mn_stover: n=61, status=success
## Mn_stover_mg: n=61, status=success
## STW50: n=61, status=success
## STW60: n=57, status=success
## EW: n=56, status=success
## EL: n=56, status=success
## ED: n=56, status=success
## KRN: n=56, status=success
## CD: n=56, status=success
## TKW: n=56, status=success
## KW50: n=56, status=success
## TKN: n=56, status=success
## HI: n=56, status=success
## Ca_grain: n=56, status=success
## Fe_grain: n=56, status=success
## K_grain: n=56, status=success
## Mg_grain: n=56, status=success
## Mn_grain: n=56, status=success
## P_grain: n=56, status=success
## S_grain: n=56, status=success
## Zn_grain: n=56, status=success
## Ca_grain_mg: n=56, status=success
## Fe_grain_mg: n=56, status=success
## K_grain_mg: n=56, status=success
## Mg_grain_mg: n=56, status=success
## Mn_grain_mg: n=56, status=success
## P_grain_mg: n=56, status=success
## S_grain_mg: n=56, status=success
## Zn_grain_mg: n=56, status=success
## Ca_total_mg: n=56, status=success
## Fe_total_mg: n=56, status=success
## K_total_mg: n=56, status=success
## Mg_total_mg: n=56, status=success
## P_total_mg: n=56, status=success
## S_total_mg: n=56, status=success
## Zn_total_mg: n=56, status=success
## Ca_HI: n=56, status=success
## Fe_HI: n=56, status=success
## K_HI: n=56, status=success
## Mg_HI: n=56, status=success
## P_HI: n=56, status=success
## S_HI: n=56, status=success
## Zn_HI: n=56, status=success
## Ca_ratio: n=56, status=success
## Fe_ratio: n=56, status=success
## K_ratio: n=56, status=success
## Mg_ratio: n=56, status=success
## P_ratio: n=56, status=success
## S_ratio: n=56, status=success
## Zn_ratio: n=56, status=success
## Mn_total_mg: n=54, status=success
## Mn_HI: n=54, status=success
## Mn_ratio: n=54, status=success# Check for missing data patterns in final dataset
missing_final <- sp_corrected_combined %>%
select(all_of(common_phenos)) %>%
summarise(across(everything(), ~ sum(is.na(.)))) %>%
pivot_longer(everything(), names_to = "phenotype", values_to = "missing_count") %>%
mutate(
total_plots = nrow(sp_corrected_combined),
missing_pct = round(missing_count / total_plots * 100, 1)
) %>%
arrange(missing_pct)
cat("\n=== Missing Data in Final Dataset ===\n")##
## === Missing Data in Final Dataset ===print(missing_final)## # A tibble: 72 × 4
## phenotype missing_count total_plots missing_pct
## <chr> <int> <int> <dbl>
## 1 PH 0 127 0
## 2 DTA 0 127 0
## 3 DTS 0 127 0
## 4 STW40 2 127 1.6
## 5 STWHV 2 127 1.6
## 6 Ca_stover 2 127 1.6
## 7 K_stover 2 127 1.6
## 8 Mg_stover 2 127 1.6
## 9 P_stover 2 127 1.6
## 10 S_stover 2 127 1.6
## # ℹ 62 more rows# Key findings summary
if (exists("treatment_effects") && exists("interaction_effects")) {
cat("\n=== Key Findings Summary ===\n")
# Most significant treatment effects
top_treatment <- treatment_effects %>%
filter(p_adj_fdr < 0.05) %>%
arrange(p_adj_fdr) %>%
head(3)
if (nrow(top_treatment) > 0) {
cat("Top treatment effects (-P vs +P):\n")
for (i in 1:nrow(top_treatment)) {
trait <- top_treatment$trait[i]
est <- round(top_treatment$estimate[i], 3)
pval <- format(top_treatment$p_adj_fdr[i], scientific = TRUE, digits = 2)
direction <- ifelse(est < 0, "decreased", "increased")
cat(sprintf(" %s: %s by %s units (FDR = %s)\n", trait, direction, abs(est), pval))
}
}
# Most significant interaction effects
top_interactions <- interaction_effects %>%
filter(p_adj_fdr < 0.05) %>%
arrange(p_adj_fdr) %>%
head(3)
if (nrow(top_interactions) > 0) {
cat("\nTop genotype × treatment interactions:\n")
for (i in 1:nrow(top_interactions)) {
trait <- top_interactions$trait[i]
est <- round(top_interactions$estimate[i], 3)
pval <- format(top_interactions$p_adj_fdr[i], scientific = TRUE, digits = 2)
direction <- ifelse(est < 0, "more sensitive", "less sensitive")
cat(sprintf(" %s: Inv4m is %s to -P than CTRL (interaction = %s, FDR = %s)\n",
trait, direction, est, pval))
}
}
}##
## === Key Findings Summary ===
## Top treatment effects (-P vs +P):
## P_stover: decreased by 1592.099 units (FDR = 1.6e-24)
## P_ratio: increased by 2.092 units (FDR = 1.1e-15)
## P_stover_mg: decreased by 185.998 units (FDR = 1.1e-15)