Code
library(tidyverse)
library(readxl)
library(MuMIn)Что обуславливает возникновение неблагоприятных условий?
К неблагоприятным (для мелких млекопитающих) отнесли дни, когда снег покрывал менее половины поверхности почвы, был слежавшимся, иногда со льдом в толще или на поверхности. С какими параметрами ассоциированы такие дни?
Шаг 1: Загрузка пакетов
Шаг 2: Кастомные функции
read_files - функция, которая считывает исходные файлы по обычной схеме, но также отмечает переход температур через 0 градусов (если произошел переход, то +1 в столбце)
Code
read_files <- function(x){
read_excel(x) |>
select(Local_time, Temp = T, Sn_description = `E'`, sss) |>
mutate(Local_time = as_datetime(Local_time, format = "%d.%m.%Y %R")) |>
mutate(Year = year(Local_time),
Month = as.factor(month(Local_time)),
Day = mday(Local_time)) |>
filter(is.na(Temp) == F) |>
mutate(
more_then_zero = Temp > 0, # логический столбец, параметр больше 0? (TRUE | FALSE)
chng1 = cumsum(more_then_zero != lag(more_then_zero, def = first(more_then_zero)))
) |>
group_by(Year, Month, Day) |>
mutate(max_chng = max(chng1),
min_chng = min(chng1),
diff = max_chng-min_chng) |>
summarise(Sn = mean(sss, na.rm = T),
Sn_description = as.factor(str_flatten(Sn_description, ", ", na.rm = T)),
Tavg = mean(Temp, na.rm = T),
diff = max(diff))
}
daily_temp <- function(x){
read_excel(x) |>
select(Local_time, Temp = T, Sn_description = `E'`, sss) |>
mutate(Local_time = as_datetime(Local_time, format = "%d.%m.%Y %R")) |>
mutate(Year = year(Local_time),
Month = as.factor(month(Local_time)),
Day = mday(Local_time)) |>
filter(is.na(Temp) == F)
}Шаг 3: Импорт необходимых файлов и проведение расчетов
На выходе получаем файл, где для каждого дня есть данные о глубине снежного покрова (Sn), качестве снежного покрова (Sn_description, елси 0, то условия благоприятные, если 1, то неблагоприятные), среднесуточной температуре (Tavg) и отметка о том, пересекали ли температуры 0 градусов в течение суток.
Code
paths <- list.files("../initial_data/climate/2005_2023", pattern = "[.]xls$", full.names = TRUE) # Просканировали все файлы в директории
total_2005_2023 <- paths |>
map_df(read_files)
s <- levels(as.factor(total_2005_2023$Sn_description))[c(4, 5,6)] # Отбор уровней фактора, когда снег = лед, или снег не покрывает всю поверхность почвы
total_2005_2023 <- total_2005_2023 |>
mutate(Sn_description = ifelse(Sn_description %in% s,1,0)) |>
mutate(Sn = case_when(is.na(Sn) ~ 0, .default = Sn))
total_2005_2023Шаг 4: Построение модели
Зависимая переменная бинарная (1 - неблагоприятные условия, 0 - благоприятные).
Предикторы - глубина снежного покрова (Sn), среднесуточная температура (Tavg), количество переходов температур через 0 градусов в течение суток (diff). Выведем все коэффициенты модели и построим доверительные интервалы (для периода с 2005 по 2023 гг)
Code
# Logreg ------------------------------------------------------------------
for_log_reg <- total_2005_2023 |>
ungroup() |>
select(Sn_description, Sn, Tavg, diff)
model <- glm(Sn_description~., for_log_reg, family = "binomial")
summary(model)
Call:
glm(formula = Sn_description ~ ., family = "binomial", data = for_log_reg)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -3.019350 0.096926 -31.151 < 2e-16 ***
Sn -0.092481 0.009137 -10.122 < 2e-16 ***
Tavg -0.043577 0.006070 -7.179 7.02e-13 ***
diff 1.260204 0.076101 16.560 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 2016.7 on 6704 degrees of freedom
Residual deviance: 1487.5 on 6701 degrees of freedom
AIC: 1495.5
Number of Fisher Scoring iterations: 9Code
confint(model) 2.5 % 97.5 %
(Intercept) -3.21384562 -2.83363921
Sn -0.11189942 -0.07589912
Tavg -0.05528946 -0.03144222
diff 1.11146347 1.41003728Теперь построим график количества неблагоприятных дней
Code
a <- total_2005_2023 |>
filter(Year > 2007,
Sn_description == 1) |>
filter(diff>1) |>
group_by(Year, Month) |>
summarise(N = n()) |>
mutate(cond = case_when(Year > 2012 ~ "good", .default = "bad"))
variability <- ggplot(a, aes(Year, N))+
geom_col() +
theme_minimal(base_size = 18)+
labs(y = "Дней",
x = "Год")+
scale_y_continuous(breaks = c(2,5,7))
print(variability)
Построим такую же модель, но для данных до 2005 года. Количество переходов температур в течение суток через 0 градусов посчитать не получится (данные изначально суточного разрешения). Выведем все коэффициенты и построим доверительные интервалы.
Code
tavg_1961_2004 <- read_excel("../initial_data/climate/1961_2005/23776_TTTR.xlsx") |>
select(
Year = год,
Month = месяц,
Day = день,
Tavg = Тср)
snow_quality_1961_2004 <- read.csv2("../initial_data/climate/1961_2005/23776_snow.csv", fileEncoding = "windows-1251") |>
select(
Station = Станция,
Year = Год,
Month = Месяц,
Day = День,
Sn = Высота_снежного_покрова,
Sn_description = Снежный_покров_.степень_покрытия,
Q1,
Q2,
Q3
)|>
filter(Q1 == 0) |> # по описанию данных 0 - данные о снежном покрове верные
filter(Sn_description != 99) |>
mutate(Sn_description = ifelse(Sn_description <= 5,1,0)) |>
left_join(tavg_1961_2004, by = c("Year", "Month", "Day")) |>
select(Year, Month, Day, Tavg, Sn, Sn_description) |>
filter(Year %in% c(1976:1994)) |>
filter(Month %in% c(4,5,10,11)) |>
select(Sn_description, Sn, Tavg) |>
mutate(Sn_description = as.factor(Sn_description)) |>
filter(is.na(Tavg)==F)
model_XX <- glm(Sn_description~., snow_quality_1961_2004, family = "binomial")
summary(model_XX)
Call:
glm(formula = Sn_description ~ ., family = "binomial", data = snow_quality_1961_2004)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.76001 0.20496 -3.708 0.000209 ***
Sn -0.30244 0.04863 -6.219 5.01e-10 ***
Tavg 0.31908 0.04177 7.639 2.19e-14 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 654.86 on 1737 degrees of freedom
Residual deviance: 329.99 on 1735 degrees of freedom
AIC: 335.99
Number of Fisher Scoring iterations: 10Code
confint(model_XX) 2.5 % 97.5 %
(Intercept) -1.1664476 -0.3612355
Sn -0.4068428 -0.2158333
Tavg 0.2415705 0.4056723Как видим, чем меньше глубина снега, тем выше шанс того, что условия будут неблагоприятными. Чем больше переходов температур через 0 градусов в течение суток, тем шанс возникновения неблагоприятных условий тоже выше. Влияние среднесуточных температур нелинейно: для данных до 2005 года чем выше среднесуточная температура, тем выше шанс того, что условия будут неблагоприятными. Для периода после 2005 года влияние обратное.