Pulgarin-2026-ips-typographus-finland-damage-model/codeGithub.R at main · alexpulga/Pulgarin-2026-ips-typographus-finland-damage-model · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
############################################################
### Code for the paper:
### Pulgarín Díaz, J.A., Melin, M., Mehtätalo, L., Polade, S., Aalto, J., Peltola, H., Tikkanen, O.-P.
### Stand, landscape and climatic attributes contributing to the probability of Ips typographus damage in Finland.
### https://doi.org/10.1016/j.foreco.2025.123436
############################################################

# download data from https://zenodo.org/uploads/19150352
# Data for the paper:
# Pulgarín Díaz, J.A., Melin, M., Mehtätalo, L., Polade, S., Aalto, J., Peltola, H., Tikkanen, O.-P.
# Stand, landscape and climatic attributes contributing to the probability of Ips typographus damage in Finland.
# https://doi.org/10.1016/j.foreco.2025.123436
#
# See also the README.txt file for details!
#
# Some practical tipps:
# - Check the version of the packages at the end of this file.
#
# @author: Pulgarín Díaz, J.A.


# Startign point

#############################################################
# Model fitting
#############################################################

#Load the dataset as "data" from: https://zenodo.org/uploads/19150352

library("glmmTMB")
model <- glmmTMB(  damage ~   meandiameter + dWpy + dWppy+ dBpy + dBppy +  dcc  + otherVol + volSprMosa +
                       prPre + prCur + hotPre + hotCur + GDDSPre + GDDSCur +
                       (1  | year) +(1  | province) , family= binomial, data= data)

#############################################################
# Model fit assessment----
#############################################################

# Overdispersion
sum((residuals(model, type = "pearson"))^2) /df.residual(model)

library("DHARMa")
simulatedRes  <-  simulateResiduals(fittedModel = model)

# Time dependency
res = recalculateResiduals(simulatedRes, group = ( data$year)) # como numerico o no dan lo mismo en la prueba
testTemporalAutocorrelation(res, time = unique(( data$year) ), plot=TRUE  )

# Time dependency
library("gstat")
library("sp")
dat = data.frame(res = residuals(simulatedRes), x = data$x, y = data$y)
coordinates(dat) = ~x+y
vario = variogram(res~1, data = dat, cutoff = 2000)
plot(vario)

# Variance inflation factor
library("performance")
check_collinearity(model)

#coefficients of determination following Nakagawa et al., (2017)
r2(model)

#  Area under the receiver operating characteristic curve
library("pROC")
observed <- data$damage
pred_probs <- predict(model, type = "response")
roc_obj <- roc(observed, pred_probs)
auc(roc_obj)
plot(roc_obj, print.auc = TRUE)

#############################################################
# End
# @author: Pulgarín Díaz, J.A.

#############################################################
# Metadata about the R session at run time:

#R version 4.4.3 (2025-02-28)
#Platform: x86_64-pc-linux-gnu
#Running under: Rocky Linux 8.10 (Green Obsidian)

#Matrix products: default
#BLAS/LAPACK: /usr/lib64/libopenblasp-r0.3.15.so;  LAPACK version 3.9.0

#locale:
#   [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C               LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8     LC_MONETARY=en_US.UTF-8
#[6] LC_MESSAGES=en_US.UTF-8    LC_PAPER=en_US.UTF-8       LC_NAME=C                  LC_ADDRESS=C               LC_TELEPHONE=C
#[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C

#attached base packages:
#   [1] stats     graphics  grDevices utils     datasets  methods   base