library(openxlsx)
library(MASS)
library(plotrix)
library(agricolae)
library(patchwork)
library(tidyverse)
library(magrittr)
library(performance)
library(corrplot)E11 Cypripedium - Analyse signifikanter Umweltparameter
1. Vorbereitende Schritte
Genutze R-Packages
Selbsterstellte Funktionen
Nutzung mit freundlicher Genehmigung von Markus Bernhardt-Römermann.
invlogit <- function(x) return(round(exp(x) / (1 + exp(x)), 3))
rename.letters <- function(x) {
x$groups <- as.character(x$groups)
if (grepl("a", x$groups[1])) {
new.letters <- x$groups
} else {
which.letters <- unique(strsplit(paste(x$groups, collapse = ""), "")[[1]])
which.letters <- which.letters[order(which.letters)]
tab.translate <- data.frame(orig = unique(strsplit(paste(x$groups, collapse = ""), "")[[1]]), new = NA)
tab.translate$orig <- as.character(tab.translate$orig)
tab.translate$new[1] <- "a"
which.letters <- which.letters[-(which.letters %in% tab.translate$new[1])]
for (j in 2:dim(tab.translate)[1]) { # j=3
if (!(tab.translate[j, "orig"] %in% which.letters)) {
tab.translate[j, "new"] <- which.letters[order(which.letters)][1]
} else {
tab.translate[j, "new"] <- tab.translate[j, "orig"]
}
which.letters <- which.letters[-(which.letters %in% tab.translate$new[j])]
}
new.letters <- rep(NA, dim(x)[1])
for (k in 1:dim(x)[1]) { # k = 1
new.letters[which(x$groups %in% tab.translate$orig[k])] <- tab.translate$new[k]
}
if (any(is.na(new.letters))) {
for (l in which(is.na(new.letters))) { # l=3
multi.letter <- unique(strsplit(paste(x$groups[l], collapse = ""), "")[[1]])
for (m in 1:length(multi.letter)) { # m = 1
multi.letter[m] <- tab.translate$new[which(tab.translate$orig %in% multi.letter[m])]
}
multi.letter <- multi.letter[order(multi.letter)]
new.letters[l] <- paste0(multi.letter, collapse = "")
}
}
}
x$groups <- new.letters
return(x)
}
pplot_reg.lines.single.plots <- function(model_final, df.year, year) {
expl.vars <- data.frame(orig.name = colnames(model_final$model), stringsAsFactors = FALSE)
expl.vars$new.name <- var.names[match(expl.vars$orig.name, vars$parameters)]
if (any(sapply(model_final$model, class) == "factor")) {
factor.vars <- colnames(model_final$model)[which(sapply(model_final$model, class) == "factor")]
expl.vars <- expl.vars[(expl.vars$orig.name != factor.vars), ]
}
for (i in 2:dim(expl.vars)[1]) { # i=2
plots_reg.lines.single <- ggplot(data = df.year, aes_string(y = expl.vars$orig.name[1], x = expl.vars$orig.name[i])) +
geom_point(colour = "black", alpha = 0.25, shape = 16, size = 0.6) +
geom_smooth(method = MASS::glm.nb, formula = y ~ x, color = "darkgreen", se = TRUE, fill = "green", alpha = 0.35) +
labs(x = expl.vars$new.name[i], y = expl.vars$new.name[1]) +
theme(panel.background = element_blank(),
panel.border = element_rect(colour = "black", fill = NA),
plot.margin = unit(c(0, 0.05, 0.2, 0), "cm"),
legend.position = "none")
#### ploting
file_name <- paste0("Plots/regline_", expl.vars$orig.name[1], "_", expl.vars$orig.name[i], "_", year, "_alternative", ".png")
ggsave(filename = file_name, plot = plots_reg.lines.single, width = 7.4, height = 7.4, units = "cm")
}
}
write_model_table <- function(model.result = NULL, file.name = NULL) {
if (file.name %in% list.files()) {
wb.model.tables <- loadWorkbook(xlsxFile = file.name)
} else {
wb.model.tables <- createWorkbook()
}
dep_var <- names(attributes(summary(model.result)$term)$dataClasses[1])
if (grepl(pattern = "(", dep_var, fixed = TRUE)) {
start.p <- regexpr(pattern = "(", text = dep_var, fixed = TRUE)[1] + 1
if (grepl(pattern = ",", dep_var, fixed = TRUE)) {
end.p <- regexpr(pattern = ",", text = dep_var, fixed = TRUE)[1] - 1
} else {
end.p <- regexpr(pattern = ")", text = dep_var, fixed = TRUE)[1] - 1
}
dep_var <- str_sub(dep_var, start = start.p, end = end.p)
}
coef_table <- as.data.frame(summary(model.result)$coefficients)
coef_table[, 1] <- round(coef_table[, 1], 3)
coef_table[, 2] <- round(coef_table[, 2], 3)
coef_table[, 3] <- round(coef_table[, 3], 3)
coef_table[, 4] <- round(coef_table[, 4], 4)
style_border <- createStyle(border = "TopBottomLeftRight ", borderStyle = "thin", borderColour = "black")
style_header <- createStyle(border = "TopBottomLeftRight", fgFill = "lightgrey")
if (length(names(wb.model.tables) != 0)) {
if (dep_var %in% names(wb.model.tables))
removeWorksheet(wb.model.tables, sheet = dep_var)
}
addWorksheet(wb.model.tables, sheetName = dep_var)
writeData(wb.model.tables, sheet = dep_var, headerStyle = style_header, rowNames = TRUE, x = coef_table)
addStyle(wb.model.tables, sheet = dep_var, style = style_border, rows = 1:(dim(coef_table)[1] + 1), cols = 1:(dim(coef_table)[2] + 1), stack = TRUE, gridExpand = TRUE)
# save table with model results
saveWorkbook(wb.model.tables, file = file.name, overwrite = TRUE)
}Einlesen der Daten von 2019 und 2024
data_2024 <- read.xlsx("data_2024.xlsx", sheet = 1)
str(data_2024)'data.frame': 80 obs. of 26 variables:
$ Labnumber : num 1 2 3 4 5 6 7 8 9 10 ...
$ YourID : chr "WI1" "WI2" "WI3" "WI4" ...
$ Zone.UTM : chr "32U" "32U" "32U" "32U" ...
$ East : num 11.5 11.5 11.5 11.5 11.5 ...
$ North : num 50.9 50.9 50.9 50.9 50.9 ...
$ nb.stem : num 8 8 5 4 23 15 2 8 13 3 ...
$ stem.per.sqm: num 55.1 63.8 31.7 74.1 115 ...
$ area.bunch : num 1452 1254 1575 540 2000 ...
$ nb.flower : num 3 6 4 3 21 11 1 3 9 0 ...
$ prop.flower : num 37.5 75 80 75 91.3 ...
$ stem.height : num 37 34.8 31.9 32.2 34.2 ...
$ l.leaves : num 13.5 12 13.5 12.2 14.6 ...
$ w.leaves : num 10.75 8.25 8.55 7.9 9.25 ...
$ management2 : chr "Wald" "Wald" "Wald" "Wald" ...
$ exposition : num 176 152 133 162 166 184 177 165 178 171 ...
$ slope : num 60 72 59 74 60 57 65 60 68 65 ...
$ soil_depth : num 25 15.7 22.3 15.3 24.3 ...
$ soil_water : num 22.8 25.6 29.4 25.8 31.1 ...
$ PAR : num 21.4 26.8 15.8 21.6 20.1 ...
$ HL_cover : num 40 60 85 63 80 45 80 75 75 20 ...
$ SL_cover : num 20 70 40 60 30 40 3 40 40 90 ...
$ TL_cover : num 20 15 40 40 10 0 0 20 20 20 ...
$ soil_cover : num 5 40 5 30 20 5 10 10 5 20 ...
$ moss_cover : num 90 40 40 30 40 70 80 80 60 40 ...
$ vh.max : num 43 45 77 27 40 50 52 48 46 40 ...
$ vh.90 : num 15 16 21 12 20 20 10 25 21 15 ...data_2019 <- read.xlsx("data_2019.xlsx", sheet = 1)
str(data_2019)'data.frame': 91 obs. of 27 variables:
$ Labnumber : num 60 57 59 56 53 58 55 54 51 52 ...
$ YourID : chr "01F" "02F" "03F" "04F" ...
$ Zone.UTM : chr "32U" "32U" "32U" "32U" ...
$ East : num 677065 677072 677067 677066 677069 ...
$ North : num 5643900 5643904 5643900 5643903 5643909 ...
$ nb.stem : num 19 11 11 14 8 22 3 26 2 6 ...
$ stem.per.sqm: num 43.2 149.7 201.5 40.7 124.2 ...
$ area.bunch : num 4399 735 546 3440 644 ...
$ nb.flower : num 9 11 8 7 8 8 1 0 2 2 ...
$ prop.flower : num 47.4 100 72.7 50 100 ...
$ stem.height : num 46 43.5 43.5 47.5 54 47 47.5 27 55 43.5 ...
$ FvFm : num 0.755 0.77 0.785 0.765 0.745 0.81 0.755 0.78 0.755 0.76 ...
$ PI : num 0.88 0.855 1.69 1.315 0.93 ...
$ SLA : num 35.5 32.2 34.6 28.7 27.9 ...
$ LDMC : num 162 161 163 175 181 ...
$ management : chr "CC" "CC" "CC" "CC" ...
$ management2 : chr "CC" "CC" "CC" "CC" ...
$ exposition : num 169 170 163 167 162 180 175 156 190 173 ...
$ slope : num 35 24 19 28 30 33 25 29 26 25 ...
$ soil_depth : num 14.2 8.2 10.2 13.2 14.3 14.2 14 12.3 11.9 13.3 ...
$ soil_water1 : num 12.6 11.3 18.1 18.8 14.1 20.8 12.8 20.5 11.3 15.1 ...
$ soil_water2 : num 6.83 3.1 7.03 17.2 5.4 ...
$ PAR : num 30 29.2 23.3 30.3 34.4 ...
$ HL_cover : num 70 60 30 80 90 60 30 40 40 50 ...
$ SL_cover : num 0 0 20 0 0 0 0 0 0 0 ...
$ soil_cover : num 0 10 40 10 3 10 0 10 10 5 ...
$ moss_cover : num 30 30 30 10 5 30 70 30 50 40 ...Vorbereitung der Daten
data_2024 <- data_2024 |>
mutate(leaf.area = l.leaves*w.leaves)
data_2024$management2 <- factor(data_2024$management2,
levels = c("Wald", "Hang", "Buche", "Fichte"),
labels = c("Kiefer (n. entb.)", "Kiefer (entb.)", "Buche", "Fichte"))
data_2019$management2 <- factor(data_2019$management2,
levels = c("CC", "5yr", "beech_forest", "spruce_forest"),
labels = c("Kiefer (n. entb.)", "Kiefer (entb.)", "Buche", "Fichte"))
data_2019 <- data_2019 %>%
mutate(soil_water = rowMeans(select(., soil_water1, soil_water2)))
# area.bunch from cm² to m²
data_2024$area.bunch <- data_2024$area.bunch / 10000
data_2019$area.bunch <- data_2019$area.bunch / 10000
parameters <- colnames(data_2024[, c(6:13, 15:27)])
var.names <- c("Sprossanzahl pro Horst",
"Sprossanzahl pro m² Horstgröße",
"Horstgröße [m²]",
"Anzahl blühender Sprosse",
"Anteil blühender Sprosse [%]",
"Sprosshöhe [cm]",
"Blattlänge [cm]",
"Blattbreite [cm]",
"Exposition [°]",
"Neigung [%]",
"Bodentiefe [cm]",
"Bodenfeuchte [%]",
"PAR [%]",
"Deckung Krautschicht [%]",
"Deckung Strauchschicht [%]",
"Deckung Baumschicht [%]",
"Deckung offener Boden [%]",
"Deckung Moosschicht [%]",
"Maximale Vegetationshöhe [cm]",
"Höhe 90% der Vegetation [cm]",
"Blattfläche [cm²]")
vars <- data.frame(parameters, var.names)2. Ökologische Modellierung 2024
Korrelationsmatrix
subsetcorr <- subset(data_2024, select=c(nb.stem, stem.per.sqm, prop.flower, area.bunch, nb.flower, stem.height, exposition, slope, soil_depth, soil_water, PAR, HL_cover, SL_cover, soil_cover, moss_cover, TL_cover, vh.max, vh.90, l.leaves, w.leaves))
#add leaf area
subsetcorr <- subsetcorr %>%
mutate(leaf.area = w.leaves * l.leaves)
#make matrix
cormatrix <- cor(subsetcorr)
corpvalues <- cor.mtest(subsetcorr)
#matrix with numbers
png("Plots/corrmatrix_numbers_2024.png", width = 25, height = 25, units = "cm", res = 300)
corrplot(cormatrix, p.mat=corpvalues$p, type="upper", method="number")
dev.off()
#matrix with circles
png("Plots/corrmatrix_circle_2024.png", width = 25, height = 25, units = "cm", res = 300)
corrplot(cormatrix, p.mat=corpvalues$p, type="upper", method="circle")
dev.off()corrplot(cormatrix, p.mat=corpvalues$p, type="upper", method="number")
corrplot(cormatrix, p.mat=corpvalues$p, type="upper", method="circle")
Modell zu Sprossanzahl pro Horst (nb.stem)
Erstellen und Simplifikation des Modells. Außerdem Exportieren der Ergebnisse in Excel-Datei.
glm1 <- glm.nb(nb.stem ~ management2 + exposition + slope + soil_depth + soil_water + PAR + HL_cover + SL_cover + TL_cover + soil_cover + moss_cover + vh.max + vh.90, data = data_2024)
summary(glm1)
glm2 <- stepAIC(glm1)
summary(glm2)
glm3 <- update(glm2, .~. -soil_water)
anova(glm2, glm3)
summary(glm3)
model_final <- glm3write_model_table(model.result = model_final, file.name = "Model.result.tables_2024.xlsx")summary(model_final)
Call:
glm.nb(formula = nb.stem ~ management2 + exposition + HL_cover +
TL_cover, data = data_2024, init.theta = 3.433064531, link = log)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -1.120490 1.331376 -0.842 0.40001
management2Kiefer (entb.) 0.373457 0.223904 1.668 0.09533 .
management2Buche -2.838777 1.084882 -2.617 0.00888 **
management2Fichte -0.255322 0.330233 -0.773 0.43943
exposition 0.015240 0.007419 2.054 0.03996 *
HL_cover 0.010604 0.004755 2.230 0.02573 *
TL_cover -0.011846 0.003772 -3.140 0.00169 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for Negative Binomial(3.4331) family taken to be 1)
Null deviance: 128.217 on 79 degrees of freedom
Residual deviance: 70.219 on 73 degrees of freedom
AIC: 386.66
Number of Fisher Scoring iterations: 1
Theta: 3.433
Std. Err.: 0.909
2 x log-likelihood: -370.657 Testen der Modell-Voraussetzungen.
par(mfrow = c(2, 2))
plot(model_final)
par(mfrow = c(1, 1))
check_model(model_final)
Visualisierungen
Sprosszahl ~ Alle signifikanten Umweltvariablen
# leicht bearbeitete Funktion aus der Analyse von 2019 für nichtlineare Regression
pplot_reg.lines.single.plots(model_final, data_2024, 2024)Modell zu Horstgröße (area.bunch)
Erstellen und Simplifikation des Modells. Außerdem Exportieren der Ergebnisse in Excel-Datei.
lm1 <- lm(log(area.bunch) ~ management2 + exposition + slope + soil_depth + soil_water + PAR + HL_cover + SL_cover + TL_cover + soil_cover + moss_cover + vh.max + vh.90, data = data_2024)
#plot(lm2)
summary(lm1)
lm2 <- stepAIC(lm1)
summary(lm2)
lm3 <- update(lm2, .~. -TL_cover)
anova(lm2, lm3)
summary(lm3)
lm4 <- update(lm3, .~. -moss_cover)
anova(lm3, lm4)
summary(lm4)
model_final <- lm4write_model_table(model.result = model_final, file.name = "Model.result.tables_2024.xlsx")summary(model_final)
Call:
lm(formula = log(area.bunch) ~ soil_depth + HL_cover, data = data_2024)
Residuals:
Min 1Q Median 3Q Max
-2.2718 -0.6489 -0.1654 0.6651 3.4218
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -2.664149 0.495133 -5.381 7.74e-07 ***
soil_depth -0.041884 0.019634 -2.133 0.0361 *
HL_cover 0.012596 0.005497 2.291 0.0247 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.125 on 77 degrees of freedom
Multiple R-squared: 0.1016, Adjusted R-squared: 0.07829
F-statistic: 4.355 on 2 and 77 DF, p-value: 0.01615Testen der Modell-Voraussetzungen.
par(mfrow = c(2, 2))
plot(model_final)
par(mfrow = c(1, 1))
check_model(model_final)
Visualisierungen
Horstgröße ~ Bodentiefe
lm.soil_depth <- lm(log(area.bunch) ~ soil_depth, data = data_2024)
range.x <- range(data_2024$soil_depth)
new.x <- data.frame(soil_depth = seq(from = range.x[1], to = range.x[2], 0.1))
new.y <- predict(lm.soil_depth,
newdata = new.x,
se.fit = TRUE) %>%
as.data.frame() %>%
mutate(soil_depth = seq(from = range.x[1], to = range.x[2], 0.1),
lower = (fit - 1.96*se.fit),
point.estimate = (fit),
upper = (fit + 1.96*se.fit))
plot <- ggplot(data_2024, aes(x = soil_depth, y = log(area.bunch))) +
geom_point(colour = "black",
alpha = 0.25,
shape = 16,
size = 0.6) +
geom_line(aes(x = soil_depth, y = point.estimate),
data = new.y,
colour = "darkgreen",
linewidth = 1) +
geom_ribbon(aes(x = soil_depth, ymin = lower, ymax = upper),
data = new.y,
color = NA, fill = "green",
alpha = 0.35,
inherit.aes = FALSE) +
labs(x = "Bodentiefe [cm]",
y = "Horstgröße [m²]") +
theme(panel.background = element_blank(),
panel.border = element_rect(colour = "black", fill = NA),
legend.position = "none")
plot
ggsave(filename = "Plots/regline_log(area.bunch)~soil_depth_2024.png", plot = plot, width = 7.4, height = 7.4, units = "cm")Horstgröße ~ Deckung Krautschicht
lm.HL_cover <- lm(log(area.bunch) ~ HL_cover, data = data_2024)
range.x <- range(data_2024$HL_cover)
new.x <- data.frame(HL_cover = seq(from = range.x[1], to = range.x[2], 0.1))
new.y <- predict(lm.HL_cover,
newdata = new.x,
se.fit = TRUE) %>%
as.data.frame() %>%
mutate(HL_cover = seq(from = range.x[1], to = range.x[2], 0.1),
lower = (fit - 1.96*se.fit),
point.estimate = (fit),
upper = (fit + 1.96*se.fit))
plot <- ggplot(data_2024, aes(x = HL_cover, y = log(area.bunch))) +
geom_point(colour = "black",
alpha = 0.25,
shape = 16,
size = 0.6) +
geom_line(aes(x = HL_cover, y = point.estimate),
data = new.y,
colour = "darkgreen",
linewidth = 1) +
geom_ribbon(aes(x = HL_cover, ymin = lower, ymax = upper),
data = new.y,
color = NA, fill = "green",
alpha = 0.35,
inherit.aes = FALSE) +
labs(x = "Deckung Krautschicht [%]",
y = "Horstgröße [m²]") +
theme(panel.background = element_blank(),
panel.border = element_rect(colour = "black", fill = NA),
legend.position = "none")
plot
ggsave(filename = "Plots/regline_log(area.bunch)~HL_cover_2024.png", plot = plot, width = 7.4, height = 7.4, units = "cm")Modell zu Blüherfolg (cbind(nb.flower, nb.stem))
Erstellen und Simplifikation des Modells. Außerdem Exportieren der Ergebnisse in Excel-Datei.
glm1 <- glm(cbind(nb.flower, nb.stem) ~ management2 + exposition + slope + soil_depth + soil_water + PAR + HL_cover + SL_cover + TL_cover + soil_cover + moss_cover + vh.max + vh.90, data = data_2024, family = binomial)
summary(glm1)
glm2 <- stepAIC(glm1)
summary(glm2)
glm3 <- update(glm2, .~. -soil_depth)
anova(glm2, glm3)
summary(glm3)
model_final <- glm3write_model_table(model.result = model_final, file.name = "Model.result.tables_2024.xlsx")summary(model_final)
Call:
glm(formula = cbind(nb.flower, nb.stem) ~ management2 + soil_water +
HL_cover, family = binomial, data = data_2024)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.092329 0.654805 -0.141 0.88787
management2Kiefer (entb.) -0.649906 0.233826 -2.779 0.00545 **
management2Buche 0.180868 0.286321 0.632 0.52758
management2Fichte 0.005277 0.372174 0.014 0.98869
soil_water -0.043477 0.021553 -2.017 0.04367 *
HL_cover 0.014137 0.005412 2.612 0.00899 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 78.394 on 79 degrees of freedom
Residual deviance: 61.577 on 74 degrees of freedom
AIC: 196.7
Number of Fisher Scoring iterations: 4Testen der Modell-Voraussetzungen.
par(mfrow = c(2, 2))
plot(model_final)
par(mfrow = c(1, 1))
check_model(model_final)
Visualisierungen
Blüherfolg ~ Bodenfeuchte
glm_soil_water <- glm(cbind(nb.flower, nb.stem) ~ soil_water,
family = binomial,
data = data_2024)
new.y <- data.frame(
soil_water = seq(min(data_2024$soil_water), max(data_2024$soil_water), length.out = 100)
)
new.y$fit <- predict(glm_soil_water, newdata = new.y, type = "response")
ci <- predict(glm_soil_water, newdata = new.y, type = "link", se.fit = TRUE)
new.y$lower <- plogis(ci$fit - 1.96 * ci$se.fit)
new.y$upper <- plogis(ci$fit + 1.96 * ci$se.fit)
plot <- ggplot() +
geom_point(data = data_2024,
aes(x = soil_water, y = nb.flower / nb.stem),
colour = "black",
alpha = 0.25,
shape = 16,
size = 0.6) +
geom_line(data = new.y,
aes(x = soil_water, y = fit),
colour = "darkgreen",
linewidth = 1) +
geom_ribbon(data = new.y,
aes(x = soil_water, ymin = lower, ymax = upper),
fill = "green",
alpha = 0.35) +
labs(x = "Bodenfeuchte [%]",
y = "Anteil blühender Sprosse [%]") +
theme(panel.background = element_blank(),
panel.border = element_rect(colour = "black", fill = NA),
plot.margin = unit(c(0, 0.05, 0.2, 0), "cm"),
legend.position = "none")
plot
ggsave(filename = "Plots/regline_prop.flower~soil_water_2024.png", plot = plot, width = 7.4, height = 7.4, units = "cm")Blüherfolg ~ Deckung Krautschicht
glm_HL_cover <- glm(cbind(nb.flower, nb.stem) ~ HL_cover,
family = binomial,
data = data_2024)
new.y <- data.frame(
HL_cover = seq(min(data_2024$HL_cover), max(data_2024$HL_cover), length.out = 100)
)
new.y$fit <- predict(glm_HL_cover, newdata = new.y, type = "response")
ci <- predict(glm_HL_cover, newdata = new.y, type = "link", se.fit = TRUE)
new.y$lower <- plogis(ci$fit - 1.96 * ci$se.fit)
new.y$upper <- plogis(ci$fit + 1.96 * ci$se.fit)
plot <- ggplot() +
geom_point(data = data_2024,
aes(x = HL_cover, y = nb.flower / nb.stem),
colour = "black",
alpha = 0.25,
shape = 16,
size = 0.6) +
geom_line(data = new.y,
aes(x = HL_cover, y = fit),
colour = "darkgreen",
linewidth = 1) +
geom_ribbon(data = new.y,
aes(x = HL_cover, ymin = lower, ymax = upper),
fill = "green",
alpha = 0.35) +
labs(x = "Deckung Krautschicht [%]",
y = "Anteil blühender Sprosse [%]") +
theme(panel.background = element_blank(),
panel.border = element_rect(colour = "black", fill = NA),
plot.margin = unit(c(0, 0.05, 0.2, 0), "cm"),
legend.position = "none")
plot
ggsave(filename = "Plots/regline_prop.flower~HL_cover_2024.png", plot = plot, width = 7.4, height = 7.4, units = "cm")Modell zu Sprosshöhe (stem.height)
Erstellen und Simplifikation des Modells. Außerdem Exportieren der Ergebnisse in Excel-Datei.
lm1 <- lm(log(stem.height) ~ management2 + exposition + slope + soil_depth + soil_water + PAR + HL_cover + SL_cover + TL_cover + soil_cover + moss_cover + vh.max + vh.90, data = data_2024)
#plot(lm1)
summary(lm1)
#lm2 <- stepAIC(lm1)
#summary(lm2)
#edit 2024: left out stepAIC to manually test if there is really nothing significant to model
lm2 <- update(lm1, .~. -soil_water)
anova(lm1, lm2)
summary(lm2)
lm3 <- update(lm2, .~. -soil_depth)
anova(lm2, lm3)
summary(lm3)
lm4 <- update(lm3, .~. -management2)
anova(lm4, lm3)
summary(lm4)
lm5 <- update(lm4, .~. -slope)
anova(lm4, lm5)
summary(lm5)
lm6 <- update(lm5, .~. -PAR)
anova(lm5, lm6)
summary(lm6)
lm7 <- update(lm6, .~. -exposition)
anova(lm6, lm7)
summary(lm7)
lm8 <- update(lm7, .~. -HL_cover)
anova(lm7, lm8)
summary(lm8)
lm9 <- update(lm8, .~. -vh.max)
anova(lm8, lm9)
summary(lm9)
lm10 <- update(lm9, .~. -vh.90)
anova(lm9, lm10)
summary(lm10)
lm11 <- update(lm10, .~. -SL_cover)
anova(lm10, lm11)
summary(lm11)
lm12 <- update(lm11, .~. -moss_cover)
anova(lm11, lm12)
summary(lm12)
lm13 <- update(lm12, .~. -TL_cover)
anova(lm12, lm13)
summary(lm13)
model_final <- lm13write_model_table(model.result = model_final, file.name = "Model.result.tables_2024.xlsx")summary(model_final)
Call:
lm(formula = log(stem.height) ~ soil_cover, data = data_2024)
Residuals:
Min 1Q Median 3Q Max
-1.71033 -0.22416 0.09463 0.34052 0.71761
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 3.331358 0.072592 45.891 <2e-16 ***
soil_cover -0.003864 0.002354 -1.642 0.105
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.4596 on 78 degrees of freedom
Multiple R-squared: 0.0334, Adjusted R-squared: 0.02101
F-statistic: 2.695 on 1 and 78 DF, p-value: 0.1047Modell weist keine signifikante Korrelation mit Umweltvariablen auf!
Testen der Modell-Voraussetzungen.
par(mfrow = c(2, 2))
plot(model_final)
par(mfrow = c(1, 1))
check_model(model_final)
Modell zu Sprossanzahl pro m² Horstgröße (stem.per.sqm)
Erstellen und Simplifikation des Modells. Außerdem Exportieren der Ergebnisse in Excel-Datei.
lm1 <- lm(stem.per.sqm ~ management2 + exposition + slope + soil_depth + soil_water + PAR + HL_cover + SL_cover + TL_cover + soil_cover + moss_cover + vh.max + vh.90, data = data_2024)
summary(lm1)
lm2 <- stepAIC(lm1)
summary(lm2)
lm3 <- update(lm2, .~. -PAR)
anova(lm2, lm3)
summary(lm3)
lm4 <- update(lm3, .~. -HL_cover)
anova(lm3, lm4)
summary(lm4)
model_final <- lm4write_model_table(model.result = model_final, file.name = "Model.result.tables_2024.xlsx")summary(model_final)
Call:
lm(formula = stem.per.sqm ~ management2 + vh.90, data = data_2024)
Residuals:
Min 1Q Median 3Q Max
-154.567 -39.060 -2.531 23.248 254.825
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -7.3966 23.9074 -0.309 0.757886
management2Kiefer (entb.) 78.4965 24.3436 3.225 0.001870 **
management2Buche 14.8977 23.9651 0.622 0.536064
management2Fichte 11.6178 24.4336 0.475 0.635825
vh.90 3.7545 0.9233 4.066 0.000117 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 75.64 on 75 degrees of freedom
Multiple R-squared: 0.357, Adjusted R-squared: 0.3227
F-statistic: 10.41 on 4 and 75 DF, p-value: 9.227e-07Testen der Modell-Voraussetzungen.
par(mfrow = c(2, 2))
plot(model_final)
par(mfrow = c(1, 1))
check_model(model_final)
Visualisierungen
Sprossanzahl pro m² Horstgröße ~ Höhe 90% der Vegetation
lm.vh.90 <- lm(stem.per.sqm ~ vh.90, data = data_2024)
range.x <- range(data_2024$vh.90)
new.x <- data.frame(vh.90 = seq(from = range.x[1], to = range.x[2], 0.1))
new.y <- predict(lm.vh.90,
newdata = new.x,
se.fit = TRUE) %>%
as.data.frame() %>%
mutate(vh.90 = seq(from = range.x[1], to = range.x[2], 0.1),
lower = (fit - 1.96*se.fit),
point.estimate = (fit),
upper = (fit + 1.96*se.fit))
plot <- ggplot(data_2024, aes(x = vh.90, y = stem.per.sqm)) +
geom_point(colour = "black",
alpha = 0.25,
shape = 16,
size = 0.6) +
geom_line(aes(x = vh.90, y = point.estimate),
data = new.y,
colour = "darkgreen",
linewidth = 1) +
geom_ribbon(aes(x = vh.90, ymin = lower, ymax = upper),
data = new.y,
color = NA, fill = "green",
alpha = 0.35,
inherit.aes = FALSE) +
labs(x = "Höhe 90% der Vegetation [cm]",
y = "Sprossanzahl pro m² Horstgröße") +
theme(panel.background = element_blank(),
panel.border = element_rect(colour = "black", fill = NA),
legend.position = "none")
plot
ggsave(filename = "Plots/regline_stem.per.sqm~vh.90_2024.png", plot = plot, width = 7.4, height = 7.4, units = "cm")Modell zu Blattfläche (leaf.area)
Erstellen und Simplifikation des Modells. Außerdem Exportieren der Ergebnisse in Excel-Datei.
lm1 <- lm(leaf.area ~ management2 + exposition + slope + soil_depth + soil_water + PAR + HL_cover + SL_cover + TL_cover + soil_cover + moss_cover + vh.max + vh.90, data = data_2024)
summary(lm1)
lm2 <- stepAIC(lm1)
summary(lm2)
lm3 <- update(lm2, .~. -exposition)
anova(lm2, lm3)
summary(lm3)
model_final <- lm3write_model_table(model.result = model_final, file.name = "Model.result.tables_2024.xlsx")summary(model_final)
Call:
lm(formula = leaf.area ~ soil_cover, data = data_2024)
Residuals:
Min 1Q Median 3Q Max
-72.303 -20.092 2.297 21.706 60.853
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 109.0023 4.6775 23.304 <2e-16 ***
soil_cover -0.3678 0.1516 -2.425 0.0176 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 29.61 on 78 degrees of freedom
Multiple R-squared: 0.07011, Adjusted R-squared: 0.05819
F-statistic: 5.881 on 1 and 78 DF, p-value: 0.01762Testen der Modell-Voraussetzungen.
par(mfrow = c(2, 2))
plot(model_final)
par(mfrow = c(1, 1))
check_model(model_final)
Visualisierungen
Blattfläche ~ Deckung offener Boden
lm.soil_cover <- lm(leaf.area ~ soil_cover, data = data_2024)
range.x <- range(data_2024$soil_cover)
new.x <- data.frame(soil_cover = seq(from = range.x[1], to = range.x[2], 0.1))
new.y <- predict(lm.soil_cover,
newdata = new.x,
se.fit = TRUE) %>%
as.data.frame() %>%
mutate(soil_cover = seq(from = range.x[1], to = range.x[2], 0.1),
lower = (fit - 1.96*se.fit),
point.estimate = (fit),
upper = (fit + 1.96*se.fit))
plot <- ggplot(data_2024, aes(x = soil_cover, y = leaf.area)) +
geom_point(colour = "black",
alpha = 0.25,
shape = 16,
size = 0.6) +
geom_line(aes(x = soil_cover, y = point.estimate),
data = new.y,
colour = "darkgreen",
linewidth = 1) +
geom_ribbon(aes(x = soil_cover, ymin = lower, ymax = upper),
data = new.y,
color = NA, fill = "green",
alpha = 0.35,
inherit.aes = FALSE) +
labs(x = "Deckung offener Boden [%]",
y = "Blattfläche [cm²]") +
theme(panel.background = element_blank(),
panel.border = element_rect(colour = "black", fill = NA),
legend.position = "none")
plot
ggsave(filename = "Plots/regline_leaf.area~soil_cover_2024.png", plot = plot, width = 7.4, height = 7.4, units = "cm")Modell zu Deckung Krautschicht (HL_cover)
Erstellen und Simplifikation des Modells. Außerdem Exportieren der Ergebnisse in Excel-Datei.
lm1 <- lm(HL_cover ~ management2 + exposition + slope + soil_depth + soil_water + PAR + SL_cover + TL_cover + soil_cover + moss_cover + vh.max + vh.90, data = data_2024)
summary(lm1)
lm2 <- stepAIC(lm1)
summary(lm2)
lm3 <- update(lm2, .~. -vh.max)
anova(lm2, lm3)
summary(lm3)
model_final <- lm3write_model_table(model.result = model_final, file.name = "Model.result.tables_2024.xlsx")summary(model_final)
Call:
lm(formula = HL_cover ~ management2 + SL_cover + soil_cover +
moss_cover, data = data_2024)
Residuals:
Min 1Q Median 3Q Max
-32.193 -10.497 1.817 11.529 33.159
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 90.76683 9.30330 9.756 7.21e-15 ***
management2Kiefer (entb.) -12.61287 5.58755 -2.257 0.0270 *
management2Buche -9.47061 5.77094 -1.641 0.1051
management2Fichte -35.86095 6.39038 -5.612 3.40e-07 ***
SL_cover -0.20160 0.08978 -2.245 0.0278 *
soil_cover -0.56435 0.12036 -4.689 1.25e-05 ***
moss_cover -0.24705 0.09990 -2.473 0.0157 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 15.84 on 73 degrees of freedom
Multiple R-squared: 0.5696, Adjusted R-squared: 0.5343
F-statistic: 16.1 on 6 and 73 DF, p-value: 1.051e-11Testen der Modell-Voraussetzungen.
par(mfrow = c(2, 2))
plot(model_final)
par(mfrow = c(1, 1))
check_model(model_final)
Visualisierungen
Deckung Krautschicht ~ Deckung Strauchschicht
lm.SL_cover <- lm(HL_cover ~ SL_cover, data = data_2024)
range.x <- range(data_2024$SL_cover)
new.x <- data.frame(SL_cover = seq(from = range.x[1], to = range.x[2], 0.1))
new.y <- predict(lm.SL_cover,
newdata = new.x,
se.fit = TRUE) %>%
as.data.frame() %>%
mutate(SL_cover = seq(from = range.x[1], to = range.x[2], 0.1),
lower = (fit - 1.96*se.fit),
point.estimate = (fit),
upper = (fit + 1.96*se.fit))
plot <- ggplot(data_2024, aes(x = SL_cover, y = HL_cover)) +
geom_point(colour = "black",
alpha = 0.25,
shape = 16,
size = 0.6) +
geom_line(aes(x = SL_cover, y = point.estimate),
data = new.y,
colour = "darkgreen",
linewidth = 1) +
geom_ribbon(aes(x = SL_cover, ymin = lower, ymax = upper),
data = new.y,
color = NA, fill = "green",
alpha = 0.35,
inherit.aes = FALSE) +
labs(x = "Deckung Strauchschicht [%]",
y = "Deckung Krautschicht [%]") +
theme(panel.background = element_blank(),
panel.border = element_rect(colour = "black", fill = NA),
legend.position = "none")
plot
ggsave(filename = "Plots/regline_HL_cover~SL_cover_2024.png", plot = plot, width = 7.4, height = 7.4, units = "cm")Deckung Krautschicht ~ Deckung offener Boden
lm.soil_cover <- lm(HL_cover ~ soil_cover, data = data_2024)
range.x <- range(data_2024$soil_cover)
new.x <- data.frame(soil_cover = seq(from = range.x[1], to = range.x[2], 0.1))
new.y <- predict(lm.soil_cover,
newdata = new.x,
se.fit = TRUE) %>%
as.data.frame() %>%
mutate(soil_cover = seq(from = range.x[1], to = range.x[2], 0.1),
lower = (fit - 1.96*se.fit),
point.estimate = (fit),
upper = (fit + 1.96*se.fit))
plot <- ggplot(data_2024, aes(x = soil_cover, y = HL_cover)) +
geom_point(colour = "black",
alpha = 0.25,
shape = 16,
size = 0.6) +
geom_line(aes(x = soil_cover, y = point.estimate),
data = new.y,
colour = "darkgreen",
linewidth = 1) +
geom_ribbon(aes(x = soil_cover, ymin = lower, ymax = upper),
data = new.y,
color = NA, fill = "green",
alpha = 0.35,
inherit.aes = FALSE) +
labs(x = "Deckung offener Boden [%]",
y = "Deckung Krautschicht [%]") +
theme(panel.background = element_blank(),
panel.border = element_rect(colour = "black", fill = NA),
legend.position = "none")
plot
ggsave(filename = "Plots/regline_HL_cover~soil_cover_2024.png", plot = plot, width = 7.4, height = 7.4, units = "cm")Deckung Krautschicht ~ Deckung Moosschicht
lm.moss_cover <- lm(HL_cover ~ moss_cover, data = data_2024)
range.x <- range(data_2024$moss_cover)
new.x <- data.frame(moss_cover = seq(from = range.x[1], to = range.x[2], 0.1))
new.y <- predict(lm.moss_cover,
newdata = new.x,
se.fit = TRUE) %>%
as.data.frame() %>%
mutate(moss_cover = seq(from = range.x[1], to = range.x[2], 0.1),
lower = (fit - 1.96*se.fit),
point.estimate = (fit),
upper = (fit + 1.96*se.fit))
plot <- ggplot(data_2024, aes(x = moss_cover, y = HL_cover)) +
geom_point(colour = "black",
alpha = 0.25,
shape = 16,
size = 0.6) +
geom_line(aes(x = moss_cover, y = point.estimate),
data = new.y,
colour = "darkgreen",
linewidth = 1) +
geom_ribbon(aes(x = moss_cover, ymin = lower, ymax = upper),
data = new.y,
color = NA, fill = "green",
alpha = 0.35,
inherit.aes = FALSE) +
labs(x = "Deckung Moosschicht [%]",
y = "Deckung Krautschicht [%]") +
theme(panel.background = element_blank(),
panel.border = element_rect(colour = "black", fill = NA),
legend.position = "none")
plot
ggsave(filename = "Plots/regline_HL_cover~moss_cover_2024.png", plot = plot, width = 7.4, height = 7.4, units = "cm")Verschieben der Modell-Ergebnistabelle
Direktes Abspeichern nach jedem Modell im Ordner führte zu Problemen. Daher nachträgliches Verschieben in den vorgesehenen Ordner.
file.rename("Model.result.tables_2024.xlsx", "Result Tables/Model.result.tables_2024.xlsx")3. Ökologische Modellierung 2019 + 2024
Zusammenführen von Daten
data_2024$year <- 2024
data_2019$year <- 2019
common_columns <- intersect(names(data_2024), names(data_2019))
data_2024_common <- data_2024[, c(common_columns, "year")]
data_2019_common <- data_2019[, c(common_columns, "year")]
names(data_2024_common) <- common_columns
names(data_2019_common) <- common_columns
data_merged <- rbind(data_2024_common, data_2019_common)Korrelationsmatrix
subsetcorr <- subset(data_merged, select=c(nb.stem, stem.per.sqm, prop.flower, area.bunch, nb.flower, stem.height, exposition, slope, soil_water, PAR, HL_cover, SL_cover, soil_cover, moss_cover))
#make matrix
cormatrix <- cor(subsetcorr)
corpvalues <- cor.mtest(subsetcorr)
#matrix with numbers
png("Plots/corrmatrix_numbers_merged.png", width = 25, height = 25, units = "cm", res = 300)
corrplot(cormatrix, p.mat=corpvalues$p, type="upper", method="number")
dev.off()
#matrix with circles
png("Plots/corrmatrix_circle_merged.png", width = 25, height = 25, units = "cm", res = 300)
corrplot(cormatrix, p.mat=corpvalues$p, type="upper", method="circle")
dev.off()corrplot(cormatrix, p.mat=corpvalues$p, type="upper", method="number")
corrplot(cormatrix, p.mat=corpvalues$p, type="upper", method="circle")
Modell zu Sprosszahl pro Horst (nb.stem)
Erstellen und Simplifikation des Modells. Außerdem Exportieren der Ergebnisse in Excel-Datei.
glm1 <- glm.nb(nb.stem ~ management2 + exposition + slope + soil_depth + soil_water + PAR + HL_cover + SL_cover + soil_cover + moss_cover, data = data_merged)
summary(glm1)
glm2 <- stepAIC(glm1)
summary(glm2)
model_final <- glm2write_model_table(model.result = model_final, file.name = "Model.result.tables_2024.xlsx")summary(model_final)
Call:
glm.nb(formula = nb.stem ~ management2 + soil_depth + PAR + HL_cover +
soil_cover, data = data_merged, init.theta = 2.74491993,
link = log)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 2.294495 0.279576 8.207 2.27e-16 ***
management2Kiefer (entb.) -0.181741 0.144879 -1.254 0.209684
management2Buche -0.804246 0.209838 -3.833 0.000127 ***
management2Fichte -0.454714 0.217364 -2.092 0.036443 *
soil_depth -0.040801 0.010301 -3.961 7.47e-05 ***
PAR 0.009686 0.004938 1.961 0.049824 *
HL_cover 0.006566 0.003208 2.047 0.040694 *
soil_cover -0.006735 0.003126 -2.155 0.031187 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for Negative Binomial(2.7449) family taken to be 1)
Null deviance: 263.26 on 170 degrees of freedom
Residual deviance: 163.19 on 163 degrees of freedom
AIC: 922.67
Number of Fisher Scoring iterations: 1
Theta: 2.745
Std. Err.: 0.426
2 x log-likelihood: -904.672 Testen der Modell-Voraussetzungen.
par(mfrow = c(2, 2))
plot(model_final)
par(mfrow = c(1, 1))
check_model(model_final)
Visualisierungen
Sprosszahl pro Horst ~ Alle signifikanten Umweltvariablen
pplot_reg.lines.single.plots(model_final, data_merged, "merged")Modell zu Blüherfolg (cbind(nb.flower, nb.stem))
Erstellen und Simplifikation des Modells. Außerdem Exportieren der Ergebnisse in Excel-Datei.
glm1 <- glm(cbind(nb.flower, nb.stem) ~ management2 + exposition + slope + soil_depth + soil_water + PAR + HL_cover + SL_cover + soil_cover + moss_cover, data = data_merged, family = binomial)
summary(glm1)
glm2 <- stepAIC(glm1)
summary(glm2)
glm3 <- update(glm2, .~. -PAR)
anova(glm2, glm3)
summary(glm3)
glm4 <- update(glm3, .~. -soil_cover)
anova(glm3, glm4)
summary(glm4)
model_final <- glm4write_model_table(model.result = model_final, file.name = "Model.result.tables_2024.xlsx")summary(model_final)
Call:
glm(formula = cbind(nb.flower, nb.stem) ~ soil_depth + HL_cover,
family = binomial, data = data_merged)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.64135 0.21172 -3.029 0.002452 **
soil_depth -0.02747 0.01001 -2.745 0.006053 **
HL_cover 0.00793 0.00239 3.318 0.000906 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 238.68 on 170 degrees of freedom
Residual deviance: 219.47 on 168 degrees of freedom
AIC: 532.86
Number of Fisher Scoring iterations: 4Testen der Modell-Voraussetzungen.
par(mfrow = c(2, 2))
plot(model_final)
par(mfrow = c(1, 1))
check_model(model_final)
Visualisierungen
Blüherfolg ~ Bodentiefe
glm_soil_depth <- glm(cbind(nb.flower, nb.stem) ~ soil_depth,
family = binomial,
data = data_merged)
new.y <- data.frame(
soil_depth = seq(min(data_merged$soil_depth), max(data_merged$soil_depth), length.out = 100)
)
new.y$fit <- predict(glm_soil_depth, newdata = new.y, type = "response")
ci <- predict(glm_soil_depth, newdata = new.y, type = "link", se.fit = TRUE)
new.y$lower <- plogis(ci$fit - 1.96 * ci$se.fit)
new.y$upper <- plogis(ci$fit + 1.96 * ci$se.fit)
plot <- ggplot() +
geom_point(data = data_merged,
aes(x = soil_depth, y = nb.flower / nb.stem),
colour = "black",
alpha = 0.25,
shape = 16,
size = 0.6) +
geom_line(data = new.y,
aes(x = soil_depth, y = fit),
colour = "darkgreen",
linewidth = 1) +
geom_ribbon(data = new.y,
aes(x = soil_depth, ymin = lower, ymax = upper),
fill = "green",
alpha = 0.35) +
labs(x = "Bodentiefe [cm]",
y = "Anteil blühender Sprosse [%]") +
theme(panel.background = element_blank(),
panel.border = element_rect(colour = "black", fill = NA),
plot.margin = unit(c(0, 0.05, 0.2, 0), "cm"),
legend.position = "none")
plot
ggsave(filename = "Plots/regline_prop.flower~soil_depth_merged.png", plot = plot, width = 7.4, height = 7.4, units = "cm")Blüherfolg ~ Deckung Krautschicht
glm_HL_cover <- glm(cbind(nb.flower, nb.stem) ~ HL_cover,
family = binomial,
data = data_merged)
new.y <- data.frame(
HL_cover = seq(min(data_merged$HL_cover), max(data_merged$HL_cover), length.out = 100)
)
new.y$fit <- predict(glm_HL_cover, newdata = new.y, type = "response")
ci <- predict(glm_HL_cover, newdata = new.y, type = "link", se.fit = TRUE)
new.y$lower <- plogis(ci$fit - 1.96 * ci$se.fit)
new.y$upper <- plogis(ci$fit + 1.96 * ci$se.fit)
plot <- ggplot() +
geom_point(data = data_merged,
aes(x = HL_cover, y = nb.flower / nb.stem),
colour = "black",
alpha = 0.25,
shape = 16,
size = 0.6) +
geom_line(data = new.y,
aes(x = HL_cover, y = fit),
colour = "darkgreen",
linewidth = 1) +
geom_ribbon(data = new.y,
aes(x = HL_cover, ymin = lower, ymax = upper),
fill = "green",
alpha = 0.35) +
labs(x = "Deckung Krautschicht [%]",
y = "Anteil blühender Sprosse [%]") +
theme(panel.background = element_blank(),
panel.border = element_rect(colour = "black", fill = NA),
plot.margin = unit(c(0, 0.05, 0.2, 0), "cm"),
legend.position = "none")
plot
ggsave(filename = "Plots/regline_prop.flower~HL_cover_merged.png", plot = plot, width = 7.4, height = 7.4, units = "cm")Modell zu Sprossanzahl pro m² Horstgröße (stem.per.sqm)
Erstellen und Simplifikation des Modells. Außerdem Exportieren der Ergebnisse in Excel-Datei.
lm1 <- lm(stem.per.sqm ~ management2 + exposition + slope + soil_depth + soil_water + PAR + HL_cover + SL_cover + soil_cover + moss_cover, data = data_merged)
summary(lm1)
lm2 <- stepAIC(lm1)
summary(lm2)
lm3 <- update(lm2, .~. -slope)
anova(lm2, lm3)
summary(lm3)
model_final <- lm3write_model_table(model.result = model_final, file.name = "Model.result.tables_2024.xlsx")summary(model_final)
Call:
lm(formula = stem.per.sqm ~ management2 + soil_water, data = data_merged)
Residuals:
Min 1Q Median 3Q Max
-227.66 -85.20 -20.32 38.82 604.14
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 284.980 32.624 8.735 2.53e-15 ***
management2Kiefer (entb.) 26.447 28.928 0.914 0.3619
management2Buche 26.350 34.933 0.754 0.4517
management2Fichte -66.359 30.472 -2.178 0.0308 *
soil_water -6.529 1.334 -4.892 2.34e-06 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 143.9 on 166 degrees of freedom
Multiple R-squared: 0.1699, Adjusted R-squared: 0.1499
F-statistic: 8.491 on 4 and 166 DF, p-value: 2.942e-06Testen der Modell-Voraussetzungen.
par(mfrow = c(2, 2))
plot(model_final)
par(mfrow = c(1, 1))
check_model(model_final)
Visualisierungen
Sprossanzahl pro m² Horstgröße ~ Bodenfeuchte
lm.soil_water <- lm(stem.per.sqm ~ soil_water, data = data_merged)
range.x <- range(data_merged$soil_water)
new.x <- data.frame(soil_water = seq(from = range.x[1], to = range.x[2], 0.1))
new.y <- predict(lm.soil_water,
newdata = new.x,
se.fit = TRUE) %>%
as.data.frame() %>%
mutate(soil_water = seq(from = range.x[1], to = range.x[2], 0.1),
lower = (fit - 1.96*se.fit),
point.estimate = (fit),
upper = (fit + 1.96*se.fit))
plot <- ggplot(data_merged, aes(x = soil_water, y = stem.per.sqm)) +
geom_point(colour = "black",
alpha = 0.25,
shape = 16,
size = 0.6) +
geom_line(aes(x = soil_water, y = point.estimate),
data = new.y,
colour = "darkgreen",
linewidth = 1) +
geom_ribbon(aes(x = soil_water, ymin = lower, ymax = upper),
data = new.y,
color = NA, fill = "green",
alpha = 0.35,
inherit.aes = FALSE) +
labs(x = "Bodenfeuchte [%]",
y = "Sprossanzahl pro m² Horstgröße") +
theme(panel.background = element_blank(),
panel.border = element_rect(colour = "black", fill = NA),
legend.position = "none")
plot
ggsave(filename = "Plots/regline_stem.per.sqm~soil_water_merged.png", plot = plot, width = 7.4, height = 7.4, units = "cm")Modell zu Deckung Krautschicht (HL_cover)
Erstellen und Simplifikation des Modells. Außerdem Exportieren der Ergebnisse in Excel-Datei.
lm1 <- lm(HL_cover ~ management2 + exposition + slope + soil_depth + soil_water + PAR + SL_cover + soil_cover + moss_cover, data = data_merged)
summary(lm1)
lm2 <- stepAIC(lm1)
summary(lm2)
model_final <- lm2write_model_table(model.result = model_final, file.name = "Model.result.tables_2024.xlsx")summary(model_final)
Call:
lm(formula = HL_cover ~ slope + PAR + soil_cover + moss_cover,
data = data_merged)
Residuals:
Min 1Q Median 3Q Max
-48.896 -9.758 0.702 12.382 39.135
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 54.50930 4.84252 11.256 < 2e-16 ***
slope 0.42169 0.07799 5.407 2.20e-07 ***
PAR 0.22613 0.09013 2.509 0.0131 *
soil_cover -0.59437 0.07066 -8.411 1.78e-14 ***
moss_cover -0.37585 0.05829 -6.448 1.19e-09 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 17.02 on 166 degrees of freedom
Multiple R-squared: 0.4594, Adjusted R-squared: 0.4464
F-statistic: 35.26 on 4 and 166 DF, p-value: < 2.2e-16Testen der Modell-Voraussetzungen.
par(mfrow = c(2, 2))
plot(model_final)
par(mfrow = c(1, 1))
check_model(model_final)
Visualisierungen
Deckung Krautschicht ~ Deckung Strauchschicht
lm.slope <- lm(HL_cover ~ slope, data = data_merged)
range.x <- range(data_merged$slope)
new.x <- data.frame(slope = seq(from = range.x[1], to = range.x[2], 0.1))
new.y <- predict(lm.slope,
newdata = new.x,
se.fit = TRUE) %>%
as.data.frame() %>%
mutate(slope = seq(from = range.x[1], to = range.x[2], 0.1),
lower = (fit - 1.96*se.fit),
point.estimate = (fit),
upper = (fit + 1.96*se.fit))
plot <- ggplot(data_merged, aes(x = slope, y = HL_cover)) +
geom_point(colour = "black",
alpha = 0.25,
shape = 16,
size = 0.6) +
geom_line(aes(x = slope, y = point.estimate),
data = new.y,
colour = "darkgreen",
linewidth = 1) +
geom_ribbon(aes(x = slope, ymin = lower, ymax = upper),
data = new.y,
color = NA, fill = "green",
alpha = 0.35,
inherit.aes = FALSE) +
labs(x = "Deckung Strauchschicht [%]",
y = "Deckung Krautschicht [%]") +
theme(panel.background = element_blank(),
panel.border = element_rect(colour = "black", fill = NA),
legend.position = "none")
plot
ggsave(filename = "Plots/regline_HL_cover~slope_merged.png", plot = plot, width = 7.4, height = 7.4, units = "cm")Deckung Krautschicht ~ Photosynthetisch aktive Strahlung
lm.PAR <- lm(HL_cover ~ PAR, data = data_merged)
range.x <- range(data_merged$PAR)
new.x <- data.frame(PAR = seq(from = range.x[1], to = range.x[2], 0.1))
new.y <- predict(lm.PAR,
newdata = new.x,
se.fit = TRUE) %>%
as.data.frame() %>%
mutate(PAR = seq(from = range.x[1], to = range.x[2], 0.1),
lower = (fit - 1.96*se.fit),
point.estimate = (fit),
upper = (fit + 1.96*se.fit))
plot <- ggplot(data_merged, aes(x = PAR, y = HL_cover)) +
geom_point(colour = "black",
alpha = 0.25,
shape = 16,
size = 0.6) +
geom_line(aes(x = PAR, y = point.estimate),
data = new.y,
colour = "darkgreen",
linewidth = 1) +
geom_ribbon(aes(x = PAR, ymin = lower, ymax = upper),
data = new.y,
color = NA, fill = "green",
alpha = 0.35,
inherit.aes = FALSE) +
labs(x = "Photosynthetisch aktive Strahlung [%]",
y = "Deckung Krautschicht [%]") +
theme(panel.background = element_blank(),
panel.border = element_rect(colour = "black", fill = NA),
legend.position = "none")
plot
ggsave(filename = "Plots/regline_HL_cover~PAR_merged.png", plot = plot, width = 7.4, height = 7.4, units = "cm")Deckung Krautschicht ~ Deckung offener Boden
lm.soil_cover <- lm(HL_cover ~ soil_cover, data = data_merged)
range.x <- range(data_merged$soil_cover)
new.x <- data.frame(soil_cover = seq(from = range.x[1], to = range.x[2], 0.1))
new.y <- predict(lm.soil_cover,
newdata = new.x,
se.fit = TRUE) %>%
as.data.frame() %>%
mutate(soil_cover = seq(from = range.x[1], to = range.x[2], 0.1),
lower = (fit - 1.96*se.fit),
point.estimate = (fit),
upper = (fit + 1.96*se.fit))
plot <- ggplot(data_merged, aes(x = soil_cover, y = HL_cover)) +
geom_point(colour = "black",
alpha = 0.25,
shape = 16,
size = 0.6) +
geom_line(aes(x = soil_cover, y = point.estimate),
data = new.y,
colour = "darkgreen",
linewidth = 1) +
geom_ribbon(aes(x = soil_cover, ymin = lower, ymax = upper),
data = new.y,
color = NA, fill = "green",
alpha = 0.35,
inherit.aes = FALSE) +
labs(x = "Deckung offener Boden [%]",
y = "Deckung Krautschicht [%]") +
theme(panel.background = element_blank(),
panel.border = element_rect(colour = "black", fill = NA),
legend.position = "none")
plot
ggsave(filename = "Plots/regline_HL_cover~soil_cover_merged.png", plot = plot, width = 7.4, height = 7.4, units = "cm")Deckung Krautschicht ~ Deckung Moosschicht
lm.moss_cover <- lm(HL_cover ~ moss_cover, data = data_merged)
range.x <- range(data_merged$moss_cover)
new.x <- data.frame(moss_cover = seq(from = range.x[1], to = range.x[2], 0.1))
new.y <- predict(lm.moss_cover,
newdata = new.x,
se.fit = TRUE) %>%
as.data.frame() %>%
mutate(moss_cover = seq(from = range.x[1], to = range.x[2], 0.1),
lower = (fit - 1.96*se.fit),
point.estimate = (fit),
upper = (fit + 1.96*se.fit))
plot <- ggplot(data_merged, aes(x = moss_cover, y = HL_cover)) +
geom_point(colour = "black",
alpha = 0.25,
shape = 16,
size = 0.6) +
geom_line(aes(x = moss_cover, y = point.estimate),
data = new.y,
colour = "darkgreen",
linewidth = 1) +
geom_ribbon(aes(x = moss_cover, ymin = lower, ymax = upper),
data = new.y,
color = NA, fill = "green",
alpha = 0.35,
inherit.aes = FALSE) +
labs(x = "Deckung Moosschicht [%]",
y = "Deckung Krautschicht [%]") +
theme(panel.background = element_blank(),
panel.border = element_rect(colour = "black", fill = NA),
legend.position = "none")
plot
ggsave(filename = "Plots/regline_HL_cover~moss_cover_merged.png", plot = plot, width = 7.4, height = 7.4, units = "cm")Verschieben der Modell-Ergebnistabelle
Direktes Abspeichern nach jedem Modell im Ordner führte zu Problemen. Daher nachträgliches Verschieben in den vorgesehenen Ordner.
file.rename("Model.result.tables_merged.xlsx", "Result Tables/Model.result.tables_merged.xlsx")