Min. 1st Qu. Median Mean 3rd Qu. Max.
9.067 14.227 19.479 27.234 29.035 1475.502
Code
data_2025_complete_filt <- data_2025_complete %>%filter(duration_min <=3*60)ggplot(data = data_2025_complete_filt) +geom_histogram(aes(x = duration_min), bins =30, fill ="steelblue", color ="white") +labs(x ="Duration (minutes)", y ="Count", title ="Survey Completion Time Distribution") +theme_minimal()
Code
attr(data_2025_complete$M3v, "labels")
NULL
Code
ggplot(data = data_2025_complete) +geom_boxplot(aes(x =as.factor(M3v), y = duration_min, fill =as.factor(M3v)), outlier.shape =NA) +coord_cartesian(ylim =c(0, 45)) +labs(x ="Age Group", y ="Duration (minutes)", fill="Age Group", title ="Survey Duration by Age Group") +theme_minimal()
Code
ggplot(data = data_2025_complete) +geom_boxplot(aes(x =as.factor(WERSJA_A_B), y = duration_min, fill =as.factor(WERSJA_A_B)), outlier.shape =NA) +coord_cartesian(ylim =c(0, 45)) +labs(x ="Survey version", y ="Duration (minutes)", fill="Survey version", title ="Survey Duration by Version") +theme_minimal()
Check Conditions
Check Conditions
Code
print("Check A: Drop outs")
[1] "Check A: Drop outs"
Code
# Drop out possibilities:# 1. N1 == 2table(data_2025_complete$N1, useNA ="ifany") # 167 drop outs, 161 people surveyed
# 2. NCONTROL (before DCE) =! "Stimme teilweise zu"table(data_2025_complete$NCONTROL, useNA ="ifany") # no further drop outs
2
2263
Code
# 3. Survey version on drought: A10_99 =! "trifft eher zu" (4)attr(data_2025_complete$A10_99, "label")
[1] "A10_99. To what extent does perceptible drought damage influence which forest areas you visit? Please rate each statement on a scale from 1 = “strongly disagree” to 5 = “strongly agree”. (... please select “partly agree” for this statement.)"
Code
attr(data_2025_complete$A10_99, "labels")
NULL
Code
table(data_2025_complete$A10_99, useNA ="ifany") # no further drop outs (80 (of 161) people in drought version of the survey)
4 <NA>
1106 1157
Code
# check how many people answer the first drought question (asked in both survey versions) attr(data_2025_complete$A1, "label")
[1] "A1. Have the forest areas that you regularly visit changed in the last 10 years?"
Code
attr(data_2025_complete$A1, "labels")
NULL
Code
table(data_2025_complete$A1, useNA ="ifany") # all 161 people
1 2 3 4
475 567 781 440
Code
sum(table(data_2025_complete$A1))
[1] 2263
Code
print("After survey version split 80 people left in survey version on drought")
[1] "After survey version split 80 people left in survey version on drought"
Code
attr(data_2025_complete$A4_1, "label")
[1] "A4_1. How would you rate the impact of drought damage in forests in relation to the following aspects? Please rate each statement on a scale from 1 = “Major problem” to 5 = “Major opportunity”. (… for nature and the environment)"
Code
attr(data_2025_complete$A4_1, "labels")
NULL
Code
table(data_2025_complete$A4_1, useNA ="ifany") # 80 people left in survey version on drought
1 2 3 4 5 <NA>
503 288 139 108 68 1157
Code
sum(table(data_2025_complete$A4_1))
[1] 1106
Code
# check how many answered the question after the last control questionattr(data_2025_complete$A11, "label")
[1] "A11. Are you concerned about the long-term availability of healthy forests in Germany?"
Code
attr(data_2025_complete$A11, "labels")
NULL
Code
table(data_2025_complete$A11, useNA ="ifany")
1 2 3 4 5 <NA>
302 580 128 71 25 1157
Code
sum(table(data_2025_complete$A11)) # still all 80 in survey after last control question
[1] 1106
Code
print("Check B: PRE-DCE PART")
[1] "Check B: PRE-DCE PART"
Code
print("1. N2-based branching into either N3-N10 (N2 == 1) or N11-N22 (N2==2)")
[1] "1. N2-based branching into either N3-N10 (N2 == 1) or N11-N22 (N2==2)"
Code
table(data_2025_complete$N2, useNA ="ifany") # 126 Waldbesuch einziges Ziel, 35 respondents Waldbesuch mit anderen Aktivitäten
Simple feature collection with 11 features and 5 fields
Geometry type: MULTIPOINT
Dimension: XY
Bounding box: xmin: -115.2727 ymin: 33.43144 xmax: 33.1082 ymax: 64.10578
Geodetic CRS: WGS 84
# A tibble: 11 × 6
Q26 Q26a Q26b Q26c n geometry
* <dbl> <dbl> <dbl> <dbl> <int> <MULTIPOINT [°]>
1 1 1 NA NA 184 ((24.28938 54.04921), (14.26094 52.15201), (14…
2 1 2 NA NA 88 ((6.328122 33.43144), (6.328119 33.43144), (6.…
3 1 3 NA NA 86 ((14.49389 52.27651), (14.60507 51.96806), (13…
4 2 NA 1 NA 126 ((33.1082 40.97438), (14.69295 51.09008), (13.…
5 2 NA 2 NA 70 ((32.91229 46.0284), (14.29066 51.66917), (14.…
6 2 NA 3 NA 188 ((9.519691 55.50398), (13.66483 51.03029), (12…
7 3 NA NA 1 118 ((14.63419 51.68781), (13.1917 48.97575), (8.2…
8 3 NA NA 2 187 ((-1.858765 43.34529), (14.04597 51.92539), (1…
9 3 NA NA 3 417 ((8.255586 55.65127), (14.65548 51.93426), (13…
10 3 NA NA 4 204 ((14.39528 52.39488), (14.42663 51.57859), (14…
11 3 NA NA 5 595 ((-4.148713 33.64863), (8.173572 55.85895), (2…
Code
print("Check Conditional Routing for Q29 and Q29a")
[1] "Check Conditional Routing for Q29 and Q29a"
Code
# If Q28 ∈ [1,2], the respondent should answer Q29# If Q28 == 3, they should answer Q29a# Check routing consistencydata_2025_complete %>%filter((Q28 %in%c(1,2) &is.na(Q29)) | (Q28 ==3&is.na(Q29a))) %>%nrow()
[1] 0
Code
# test for false positivesdata_2025_complete %>%filter((Q28 %in%c(1,2) &!is.na(Q29a)) | (Q28 ==3&!is.na(Q29))) %>%nrow()
print("Check if Summary Description and Q31 Were Displayed Correctly")
[1] "Check if Summary Description and Q31 Were Displayed Correctly"
Code
# Everyone who saw the summary image (status quo mock-up) should have answered Q31# If Q31 < 5 --> Q31a should be answereddata_2025_complete %>%filter(is.na(Q31_1)) %>%nrow()
Simple feature collection with 84 features and 3 fields
Geometry type: GEOMETRY
Dimension: XY
Bounding box: xmin: -115.2727 ymin: 33.43144 xmax: 33.1082 ymax: 64.10578
Geodetic CRS: WGS 84
# A tibble: 84 × 4
Q31_1 Q31a_1 n geometry
* <dbl> <chr> <int> <GEOMETRY [°]>
1 1 "" 3 MULTIPOINT ((10.58988 51…
2 1 "abgestorbene bäume" 1 POINT (14.21837 52.16511)
3 1 "es sind mehr Baumarten" 1 POINT (8.576292 48.76298)
4 1 "vom Borkenkäfer aufgefressen" 1 POINT (10.61691 51.72486)
5 2 "" 1 POINT (8.383942 50.78814)
6 2 "Bäume waren abgestorben" 1 POINT (11.13062 51.48121)
7 2 "Es war kurz nach dem Sturm und Hagel … 1 POINT (8.934975 48.61725)
8 2 "da gibt es ZU viele verschiedene Arte… 1 POINT (8.418358 48.98416)
9 2 "viele vetrocknete bäume" 1 POINT (14.21974 52.16634)
10 2 "zu viele nadelbäume" 1 POINT (9.127526 48.76064)
# ℹ 74 more rows
Code
print("Check for Internal Consistency") # can be extended
[1] "Check for Internal Consistency"
Code
# E.g. for a claimed high forest height (Q28 = 3) but low naturalness in Q30 might be interestingtest <- data_2025_complete %>%filter(Q28 ==3& Q30 %in%c(1, 2)) # 42 obsrm(test)print("Other considered attributes: Q36 show if Q35_1a==2")
[1] "Other considered attributes: Q36 show if Q35_1a==2"
Code
attr(data_2025_complete$Q35_1a, "label")
[1] "Q35. Were there any characteristics that you did not take into account when making your choices? (Yes, I did not consider)"
Code
table(data_2025_complete$Q35_1a, useNA ="ifany") # does not indicate anything
0
2263
Code
attr(data_2025_complete$Q35_97, "label")
[1] "Q35. Were there any characteristics that you did not take into account when making your choices? (No, I considered all forest characteristics)"
Code
table(data_2025_complete$Q35_97, useNA ="ifany") # 62 respondents overlooked none of the attributes
Simple feature collection with 417 features and 9 fields
Geometry type: GEOMETRY
Dimension: XY
Bounding box: xmin: -115.2727 ymin: 33.43144 xmax: 33.1082 ymax: 64.10578
Geodetic CRS: WGS 84
# A tibble: 417 × 10
Q35_97 Q36_1 Q36_2 Q36_3 Q36_4 Q36_5 Q36_6 Q36_7 n
* <dbl> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <int>
1 0 I did not notice it in the … I di… I di… I di… I di… I di… I di… 1
2 0 I did not notice it in the … I di… I di… I di… <NA> I di… I di… 2
3 0 I did not notice it in the … I di… I di… I di… <NA> <NA> <NA> 1
4 0 I did not notice it in the … I di… I di… <NA> I di… <NA> <NA> 1
5 0 I did not notice it in the … I di… I di… <NA> <NA> I di… I di… 1
6 0 I did not notice it in the … I di… I di… <NA> <NA> I di… <NA> 1
7 0 I did not notice it in the … I di… I di… <NA> <NA> <NA> <NA> 4
8 0 I did not notice it in the … I di… I fo… I fo… It i… <NA> <NA> 1
9 0 I did not notice it in the … I di… I fo… <NA> <NA> It i… I fo… 1
10 0 I did not notice it in the … I di… I fo… <NA> <NA> <NA> <NA> 1
# ℹ 407 more rows
# ℹ 1 more variable: geometry <GEOMETRY [°]>
Code
print("Forest visit: Q38 and Q39 Routing Logic")
[1] "Forest visit: Q38 and Q39 Routing Logic"
Code
# If Q38 = 1 (yes, visited other forests) --> Show Q39# If Q38 = 2, skip Q39 and all follow-up questionsattr(data_2025_complete$Q38, "labels")
NULL
Code
table(data_2025_complete$Q38)
1 2
1120 1143
Code
#Q39: If Q38 == 1 (visited other forests)#Q40a–Q42a: If Q39 == 2 (visited 2nd forest)#Q41b–Q43b: If Q39 %in% 3:6 (visited 3–5 forests)#Q41c–Q43c: If Q39 %in% 4:6 (visited 4–6 forests)# 1. Q39 should be NA if Q38 == 2 (didn't visit other forests)data_2025_complete %>%filter(Q38 ==2&!is.na(Q39)) %>%nrow()
# 2. Q40a–Q42a only if Q39 == 2data_2025_complete %>%filter(Q39 !=2& (!is.na(Q40a_1) |!is.na(Q40a_2) |!is.na(Q41a_1) |!is.na(Q41a_2) |!is.na(Q41a_9_none) |!is.na(Q41ax) |!is.na(Q42a_1))) %>%nrow()
[1] 312
Code
# 3. Q41b–Q43b only if Q39 in 3:6data_2025_complete %>%filter(!(Q39 %in%3:6) & (!is.na(Q41b_1) |!is.na(Q41b_2) |!is.na(Q42b_1) |!is.na(Q42b_2) |!is.na(Q42b_9_none) |!is.na(Q42bx) |!is.na(Q43b_1) |!is.na(Q41c_1) |!is.na(Q41c_2) |!is.na(Q42c_1) |!is.na(Q42c_2) |!is.na(Q42c_9_none) |!is.na(Q42cx) |!is.na(Q43c_1))) %>%nrow()
[1] 754
Code
# 4. Optional: Check if expected data is missing when it should be present# 4.1 Q39 == 2 --> Check for missing responses in forest data_2025_complete %>%filter(Q39 ==2& (is.na(Q40a_1) |is.na(Q40a_2) |is.na(Q41a_1) |is.na(Q41a_2) |is.na(Q41a_9_none) |is.na(Q41ax) |is.na(Q42a_1) )) %>%summarise(across(c(Q40a_1, Q40a_2, Q41a_1, Q41a_2, Q41a_9_none, Q41ax, Q42a_1),~sum(is.na(.)) ) )
Simple feature collection with 1 feature and 7 fields (with 1 geometry empty)
Geometry type: GEOMETRYCOLLECTION
Dimension: XY
Bounding box: xmin: NA ymin: NA xmax: NA ymax: NA
Geodetic CRS: WGS 84
# A tibble: 1 × 8
Q40a_1 Q40a_2 Q41a_1 Q41a_2 Q41a_9_none Q41ax Q42a_1 geometry
<int> <int> <int> <int> <int> <int> <int> <GEOMETRYCOLLECTION [°]>
1 0 0 0 0 0 0 0 GEOMETRYCOLLECTION EMPTY
# Q46 Should be answered (non-NA) only if Q45 == 1 (i.e., mushrooms were collected)# Q48 Should be answered (non-NA) only if Q47 == 1 (i.e., berries were collected)table(data_2025_complete$Q45)
# data_2025_complete %>% filter(A2 == 5 & A3 == 5) %>% # select(A2, A3, A5, A6, A8) %>%# View()#%>% nrow() # 8 times the case that questions are not triggeredsum(!(data_2025_complete$A2 ==5& data_2025_complete$A3 ==5), na.rm =TRUE) # --> but 8 is counting both survey versions. QUestion: how many got not triggered only in survey version for drought?
data_2025_complete %>%filter((A2 !=5| A5!=5) & (!is.na(A5) |!is.na(A6) |!is.na(A8))) %>%nrow() # 74 responses, so for 6 respondents questions not triggered
[1] 1073
Code
table(data_2025_complete$A5, useNA ="ifany") # answered by 74, NA:87
1 2 3 4 5 6 <NA>
27 249 505 76 60 157 1189
Code
table(data_2025_complete$A6, useNA ="ifany") # answered by 74, NA:87
1 2 3 4 5 6 <NA>
9 71 672 102 33 187 1189
Code
table(data_2025_complete$A8, useNA ="ifany") # answered by 68, NA: 93 --> lower respondents rate and higher NA than there should be (mandatory question and no further drop-outs), also the question did not appear in test version of the survey when I tried again (same for A9 but this is conditional on A9 therefore not surprising), check again
1 2 <NA>
193 881 1189
Code
table(data_2025_complete$A9_1, useNA ="ifany")
2142
-
1
Ich suche Abwechslung und neue Impressionen. Andere Gebiete bieten verschiedene Landschaften und Pflanzen.
1
ruhigere oder weniger frequentierte Orte sucht
1
Vielfalt der Natur Erklärung: Ich suche nach neuen Erlebnissen und möchte verschiedene Landschaften und Ökosysteme entdecken. Andere Gebiete bieten möglicherweise einzigartige Flora und Fauna, die ich erleben möchte.
1
"Ich habe meine Waldbesuche auf andere Gebiete verlegt, weil die Wälder in meiner Region stärker von Dürreschäden betroffen sind."
1
alles okay und gut!
1
andre eindrücke
1
asd a sda d
1
auch mal andere zu sehen.-
1
Aufgrund der besseren Infrastruktur und der vielfältigeren Wandermöglichkeiten in anderen Gebieten.
1
Aufgrund von Bauarbeiten im ursprünglichen Gebiet habe ich meine Waldbesuche umgeleitet.
1
Baum, okay.
1
Because the air quality is better in other places.
2
Consider the air quality and environment of the forest.
1
Damit ich auch andere Gebieten besuchen kann, auch um neue Pilz und Beeren Flecken zu finden. Auch wenn Dürreschäden sind, so gehe ich in denselben Wald, weil ich ihn von klein auf kannte.
1
Das weiß ich nicht
1
Das wird es nicht.
1
Der frühere Wald war zu stark frequentiert, also suchte ich einen ruhigeren Ort für Entspannung.
1
Die Änderungen der Häufigkeit, mit der ich die Wälder in Deutschland besuche, hängen hauptsächlich von der Jahreszeit und der Arbeitsrhythmus ab. Im Frühling und Sommer wirkt der Wald wie verzaubert
1
die natürliche Umwelt
1
Die Umgebung ist nicht mehr geeignet.
1
Die ursprüngliche regionale Ökologie wird überlastet oder geschädigt
1
Die Wälder in anderen Gebieten sehen besser aus.
1
Discovered other places with a better environment.
2
dort sind mehr Tiere
1
Eine andere Landschaft sehen
1
Eine bessere Umgebung anstreben
1
Einfach so
1
Enim non sed aliquip quod quod vitae ea minim magna
1
es gibt noch viel zu entdecken
1
es ist immer wichtig und sio
1
Es war ein Besuch von Verwandtschaft und ein Spaziergang in der näheren Umgebung dieser Familie
1
Facere ut perferendis voluptate aperiam
1
für balanziert
1
Gegend hat mich interessiert
1
good
2
hab ich so nicht nachgedacht
1
hatte lust drauf
1
i couldn‘t cover the distance like befor
1
I don't want to worsen the condition of the forest by further intervention
1
ich ahbe kein ehanung
1
Ich besuche die Wälder, wo ich in der Nähe bin (Urlaub oder andere Aktivitäten)
1
Ich besuche jetzt lieber Wälder in der nähe von Flüssen und Bächen.
1
Ich bleibe lieber in der Nähe
1
ich denke das es am entspannendesten sit
1
Ich habe meine Waldbesuche auf andere Gebiete verlegt, weil die ursprünglichen Gebiete zu überfüllt waren und ich eine ruhigere und entspanntere Erfahrung suche.
1
Ich habe meine Waldbesuche auf andere Gebiete verlegt, weil diese ruhiger, naturnäher oder abwechslungsreicher sind als mein bisheriger Wald. Weitere mögliche Gründe könnten sein: Bessere Wege oder Infrastruktur (z.???B. Parkplätze, Sitzbänke, ausgeschilderte Routen) Schönere Landschaft oder mehr Vielfalt in Flora und Fauna Weniger Menschen – für mehr Erholung und Stille
1
Ich habe meine Waldbesuche verlegt, um bessere Zugänglichkeit, schönere Landschaften und weniger Touristen zu finden – alles, was das Erlebnis verbessert.
1
Ich halte ausschau nach interessanteren wäldern
1
Ich möchte einen neuen Ort besuchen
1
Ich suche immer mal wieder Abwechslung
1
ich wechsel häufig die Waldgebiete um neues kennenzulernen
1
Ich wollte neue natürliche Umgebungen erkunden.
1
If it refers to the metaphorical expression "forest journey", I can understand it as referring to a certain plan, goal, or change in one's life path. The reason for changing the "forest journey" might be to seek new experiences and challenges, or because the original path has encountered unforeseen obstacles or a situation where it cannot continue to progress. Just like traveling in a forest, sometimes you may find that the path no longer suits you, or you may desire to explore different landscapes and unknown places.
1
In der Regel wegen der Luftqualitätsprobleme und der Umwelt
1
interesse
1
intuitiv
1
ist näher an meinem Wohort
1
JA
1
just to explore
1
keine
1
keine ahnung
1
Köttinger wALD
1
Lieber den weinberg.
3
mal was anders
1
mal was neues
1
Man wechselt die Waldreise, um verschiedene Waldtypen kennenzulernen – von dunkeln Nadelwäldern bis zu sonnigen Laubwäldern. Jeder hat seine eigene Flora, Fauna und Atmosphäre, und so entdeckt man die Vielfalt der Natur auf neue Weise.
1
mehr sehen
1
Mir gefällt die umgebung
1
möchte auch mal woanders gehen und neu Pilzgebiete finden
1
Nähe
1
Neue Wälder kennenzulernen und suche nach neuen Alternativen
1
neue Wege erkunden
1
none
1
play in fragmented time
1
Ruhigere oder weniger überfüllte Orte
1
Schauen Sie sich verschiedene Bäume an
1
Schöner anzusehen
1
schonung
1
schwer zu sagen
1
sehr gut
2
sind schöner
1
sind und ich habe gerade gesehen habe
1
um etwas neues zu endeken
2
Um neue Gebiete kennenzulernen
1
um Neues kennen zulernen
1
um was anderes zu sehen als sonst
1
um was neues zu sehen
1
Umwelt schlechter, neuer Ort.
1
umwelt zuliebe
1
umzug
1
Vielseitigkeit, um andere kennenzulernen.
1
Waldbesuche werden aufgrund einer Reihe von Faktoren häufig in andere Gebiete verlagert, darunter Änderungen der Waldzugänglichkeit, der Umweltbedingungen, Naturschutzbemühungen und persönliche Vorlieben für neue Erfahrungen.
1
war zu überfüllt mit BesucherInnen, sodass ich keine Ruhe im Wald genießen konnte. Deshalb wechselte ich den Ort, um die Natur in Frieden zu erleben.
1
was neues zur sehen.
1
Weil das ehemalige Waldgebiet durch übermäige Tourismus - Entwicklung beeinträchtigt wurde.
1
Weil die Hitze das nicht erlaubte
1
Weil die neuen Gebiete näher an meinem Zuhause sind.
1
Weil die Wälder dort in einem besseren Zustand sind
1
Weil es besser ist
2
weil es schöner dort ist
1
weil es schöner gebiete gibt
1
weil ich abwechslung wollte
1
weil ich den Wald nicht noch weiter zerstören wollte.
1
Weil ich gerne Abwechslung hätte
1
weil ich mochte andere walde sehen
1
Weil ich zu einem anderen Wald wechselte, um die natürliche Ökologie des aktuellen Waldgebiets zu schützen, sodass die natürliche Ökosystem-Umgebung erhalten bleibt.
1
Weil jeder Wald einzigartig ist: andere Bäume, Tiere, Landschaften. Neue Pfade bringen Frische, Entdeckungen und eine andere Verbindung zur Natur – ein Wechsel stärkt die Lebensfreude.
1
weil sie schön sind
1
weil sie zu weit weg waren
1
werden kann ich
1
Wollte was Neues sehen
1
Zu voll, Ruhe gesucht.
1
Code
# View(data_2025_complete[,c("A2","A3","A5","A6","A8")])print("Trigger A7 only if A6==1|2|4|5")
# print("Check distance to nearby village location")# summary(data_2025_complete_check$dist_check/1000)# # print("Nearby village more than 20km away?")# table(data_2025_complete_check$dist_check >= 20000)# # print("Nearby village more than 100km away?")# table(data_2025_complete_check$dist_check >= 100000)# # ggplot(data = data_2025_complete_check) +# geom_boxplot(aes(x = factor(M2a), y = dist_check / 1000), na.rm = TRUE) +# labs(x = "M2a", y = "Distance (km)")
Locations of the respondents
Code
data_2025_complete_filt <- data_2025_complete %>%mutate(lat = N8_1,lon = N8_2,# reason for invalidity (first match wins)invalid_reason =case_when(is.na(lat) &is.na(lon) ~NA_character_, # both missing = not "invalid", just missingxor(is.na(lat), is.na(lon)) ~"only_one_missing",!is.finite(lat) |!is.finite(lon) ~"non_finite", lat ==0& lon ==0~"null_island_0_0", lat <-90| lat >90~"lat_out_of_range", lon <-180| lon >180~"lon_out_of_range",TRUE~"" ),valid_combo = invalid_reason =="" ) %>%filter(valid_combo =="TRUE")data_2025_complete_filt$geometry3 <-st_sfc(Map(function(lon, lat) st_point(c(lon, lat)), data_2025_complete_filt$N8_2, # longitude data_2025_complete_filt$N8_1), # latitudecrs =4326)leaflet(data_2025_complete_filt) %>%addTiles() %>%# Adds OpenStreetMap base layeraddMarkers(data =st_set_geometry(data_2025_complete_filt, "geometry3"))
Code
data_2025_complete <- data_2025_complete %>%mutate(geometry_proj1 =st_transform(geometry, 3857), geometry_proj4 =st_transform(geometry_depart, 3857) )data_2025_complete$dist_por_for <-st_distance( data_2025_complete$geometry_proj1, data_2025_complete$geometry_proj4,by_element =TRUE)print("Distance departure point to forest")
[1] "Distance departure point to forest"
Code
summary(data_2025_complete$dist_por_for/1000)
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.053 3.314 13.666 97.268 69.507 14067.930
Socio-demographic Overview
Code
ggplot(data=data_2025_complete) +geom_bar(aes(x = M2a, fill = M2a), position ="dodge") +labs(title ="Federal States", x ="", y ="Count") +theme_minimal() +scale_fill_brewer(palette ="Set3")
Code
ggplot(data=data_2025_complete) +geom_bar(aes(x = M3_1, fill = M3_1), position ="dodge") +labs(title ="Age Distribution", x ="", y ="Count") +theme_minimal() +scale_fill_brewer(palette ="Set3")
Code
table(data_2025_complete$M3v)
1 2 3 4 5 6
82 456 627 358 338 402
Code
print("Gender")
[1] "Gender"
Code
table(data_2025_complete$M1)
1 2 3
1181 1078 4
Code
attr(data_2025_complete$M2, "labels")
NULL
Code
ggplot(data = data_2025_complete) +geom_bar(aes(x = M2), position ="dodge") +labs(title="Settlement type")
Code
attr(data_2025_complete$M4, "labels")
NULL
Code
ggplot(data = data_2025_complete) +geom_bar(aes(x = M4), position ="dodge") +labs(title="Education")
Code
ggplot(data=data_2025_complete) +geom_bar(aes(x =as.factor(M4a)), position ="dodge") +labs(title ="Life Satistfaction", x ="", y ="Count") +theme_minimal() +scale_fill_brewer(palette ="Set3")
Code
ggplot(data=data_2025_complete) +geom_bar(aes(x =as.factor(M9)), position ="dodge") +labs(title ="Net Income", x ="", y ="Count") +theme_minimal() +scale_fill_brewer(palette ="Set3")
ggplot(data=data_2025_complete) +geom_bar(aes(x=M7_1)) +labs(title="Kids in HH")
Code
ggplot(data=data_2025_complete) +geom_bar(aes(x=M10_1)) +labs(title="Living in Region since")
Flag weird entries
Code
table(data_2025_complete$red_flag_2000)
0 1
2197 66
Code
table(data_2025_complete$red_flag)
0 1
2251 12
Code
table(data_2025_complete$speeder_flag)
0
2263
Code
table(data_2025_complete$flag_out_germany)
0 1
2244 19
Compare 2017 and 2025
Code
combined <-rbind(data.frame(age =2017- data_2017_complete$M3, year =2017),data.frame(age =2025- data_2025_complete$M3_1, year =2025))# --- Option A: BINS (recommended for comparison) ---combined <- combined %>%mutate(age_group =cut( age,breaks =c(-Inf, 29, 39, 49, 60, Inf),labels =c("18–29", "30–39", "40–49", "50–60", "60+"),right =TRUE ))# Count per (year, age_group) and convert to proportion within each yearbin_props <- combined %>%filter(!is.na(age_group)) %>%count(year, age_group, name ="n") %>%group_by(year) %>%mutate(prop = n /sum(n)) %>%ungroup()ggplot(bin_props, aes(x = age_group, y = prop, fill =factor(year))) +geom_col(position ="dodge") +labs(title ="Age Distribution Comparison",x ="Age Group", y ="Proportion", fill ="Sample") +theme_minimal() +scale_fill_manual(values =c("2017"="orange", "2025"="forestgreen"))
Code
#### Compare distancedf <-bind_rows(data.frame(year ="2017", dist =factor(database_2017$dist_3)),data.frame(year ="2025", dist =factor(database_2025$dist_3)))# Ensure consistent ordering of categoriesdf$dist <-factor(df$dist, levels =sort(unique(df$dist)))freq <- df %>%count(year, dist) %>%group_by(year) %>%mutate(pct = n /sum(n))ggplot(freq, aes(x = dist, y = pct, fill = year)) +geom_col(position ="dodge") +scale_y_continuous(labels = scales::percent) +labs(x ="Distance category (in km)", y ="Relative share",title ="Comparison of forest travel distances") +scale_fill_manual(values =c("2017"="orange", "2025"="forestgreen")) +theme_minimal()
Code
#### Federal state distribution of respondents 2017 vs 2025 ####count_share <-function(df, year_label) { df %>% sf::st_drop_geometry() %>%filter(!is.na(NAME_2)) %>%filter(!is.na(NAME_1)) %>%count(NAME_1, name ="n") %>%mutate(year = year_label) %>%group_by(year) %>%mutate(share = n /sum(n)) %>%ungroup()}dist17 <-count_share(data_2017_complete, "2017")dist25 <-count_share(data_2025_complete, "2025")# Combine; ensure all counties exist for both years (fill with zeros)all_counties <-union(dist17$NAME_1, dist25$NAME_1)dist <-bind_rows(dist17, dist25) %>%complete(NAME_1 = all_counties, year =c("2017","2025"),fill =list(n =0, share =0)) %>%mutate(NAME_2 =fct_reorder(NAME_1, share, .fun = sum, .desc =TRUE))n2017 <-sum(dist17$n); n2025 <-sum(dist25$n)ggplot(dist, aes(x = NAME_1, y = share, fill = year)) +geom_col(position =position_dodge(width =0.8), width =0.7, alpha =0.85) +scale_y_continuous(labels =percent_format(accuracy =1)) +scale_fill_manual(values =c("2017"="orange", "2025"="forestgreen")) +labs(title ="Federal state distribution of respondents (relative shares)",x ="", y ="Share of sample", fill ="Year") +theme_minimal(base_size =12) +theme(legend.position ="top", panel.grid.major.y =element_blank()) +theme(axis.text.x =element_text(angle =45, hjust =1))
Forest Visit Stats
Code
attr(data_2025_complete$M2a, "labels")
NULL
Code
ggplot(data = data_2025_complete) +geom_boxplot(aes(x =factor(M2a), y = N7, fill =factor(M2a)), position ="dodge", outlier.shape =NA) +labs(title ="Forest Visits", x ="M2a", y ="N7_1") +theme_minimal() +scale_y_continuous(limits =c(0, 80))
stargazer(clg_2025_out, clg_2017_out, type ="html", column.labels =c("2025", "2017"),title ="Comparison of Models (2025 vs 2017)")
Comparison of Models (2025 vs 2017)
Estimate
s.e.
t.rat.(0)
Rob.s.e.
Rob.t.rat.(0)
type_broad
-0.543
1.935
-0.281
3.247
-0.167
type_mix
1.033
2.129
0.485
3.284
0.314
species
14.729
0.981
15.012
1.755
8.391
height
2.439
0.168
14.557
0.253
9.655
age_two
-8.542
1.822
-4.687
3.212
-2.660
age_multi
8.609
1.797
4.790
3.105
2.772
deadwood_med
3.550
1.535
2.313
2.400
1.479
deadwood_high
11.887
1.581
7.521
2.405
4.943
infrastr2
3.513
1.852
1.897
2.779
1.264
infrastr3
18.339
1.596
11.487
2.962
6.191
infrastr4
40.556
1.845
21.977
3.558
11.399
dist
-0.013
0.0003
-39.317
0.001
-17.262
ASC_sq
125.566
4.405
28.506
8.484
14.800
Comparison of Models (2025 vs 2017)
Estimate
s.e.
t.rat.(0)
Rob.s.e.
Rob.t.rat.(0)
type_broad
1.181
3.279
0.360
4.882
0.242
type_mix
6.663
3.451
1.931
5.171
1.288
species
11.942
1.488
8.025
2.629
4.543
height
2.814
0.280
10.055
0.441
6.387
age_two
0.718
3.187
0.225
5.044
0.142
age_multi
25.370
2.978
8.518
5.213
4.867
deadwood_med
3.253
2.575
1.263
3.676
0.885
deadwood_high
15.399
2.619
5.879
3.822
4.029
infrastr2
-6.279
3.230
-1.944
4.333
-1.449
infrastr3
22.950
2.493
9.207
4.708
4.875
infrastr4
39.910
2.767
14.422
5.264
7.581
dist
0.014
0.0004
30.968
0.001
11.721
ASC_sq
82.378
2.778
29.656
7.485
11.006
Code
html_2025 <-capture.output(stargazer(clg_2025_out, type ="html", title ="2025", column.labels ="2025"))html_2017 <-capture.output(stargazer(clg_2017_out, type ="html", title ="2017", column.labels ="2017"))cat('<div style="display: flex; gap: 20px;">')
[1] "A1. Have the forest areas that you regularly visit changed in the last 10 years?"
Code
attr(data_2025_complete$A1, "labels")
NULL
Code
table(data_2025_complete$A1)
1 2 3 4
475 567 781 440
Code
ggplot(data = data_2025_complete) +geom_bar(aes(x =as.factor(A1), y =after_stat(count /sum(count))),position ="dodge") +scale_y_continuous(labels = scales::percent) +labs(title ="Has forest changed in last 10 years \n(1=No, 2=Improved, 3=Deteriorated, 4=Can't tell)?",x ="",y ="Percent" ) +theme_minimal()
Code
attr(data_2025_complete$A2, "label")
[1] "A2. Did you notice signs of drought damage during your last visit to the forest?"
Code
attr(data_2025_complete$A2, "labels")
NULL
Code
table(data_2025_complete$A2)
1 2 3 4 5
387 964 436 378 98
Code
ggplot(data = data_2025_complete) +geom_bar(aes(x =as.factor(A2), y =after_stat(count /sum(count))),position ="dodge") +scale_y_continuous(labels = scales::percent) +labs(title ="Visible drought damages last forest visit (yes (1) to no (5))?",x ="",y ="Percent" ) +theme_minimal()
Code
attr(data_2025_complete$A3, "label")
[1] "A3. Do you have the impression, regardless of your last visit to the forest, that the forests in your region are increasingly affected by drought damage?"
Code
attr(data_2025_complete$A3, "labels")
NULL
Code
table(data_2025_complete$A3)
1 2 3 4 5
358 997 456 362 90
Code
ggplot(data = data_2025_complete) +geom_bar(aes(x =as.factor(A3), y =after_stat(count /sum(count))),position ="dodge") +scale_y_continuous(labels = scales::percent) +labs(title ="Perception of drought damages in region (yes (1) to no (5))?",x ="",y ="Percent" ) +theme_minimal()
Code
attr(data_2025_complete$A6, "label")
[1] "A6. Have you changed the frequency of your forest visits due to drought damage?"
Code
attr(data_2025_complete$A6, "labels")
NULL
Code
table(data_2025_complete$A6)
1 2 3 4 5 6
9 71 672 102 33 187
Code
ggplot(data = data_2025_complete %>%filter(!is.na(A6) & A6 !=6)) +geom_bar(aes(x =as.factor(A6), y =after_stat(count /sum(count))),position ="dodge") +scale_y_continuous(labels = scales::percent) +labs(title ="Changed forest visiting frequency because of drought \n(1=Less often, 3=Unchanged, 5=More often)?",x ="",y ="Percent" ) +theme_minimal()
Code
attr(data_2025_complete$A8, "label")
[1] "A8. Have you changed the forest areas you usually visit because of drought damage?"
Code
attr(data_2025_complete$A8, "labels")
NULL
Code
table(data_2025_complete$A8)
1 2
193 881
Code
ggplot(data = data_2025_complete %>%filter(!is.na(A8))) +geom_bar(aes(x =as.factor(A8), y =after_stat(count /sum(count))),position ="dodge") +scale_y_continuous(labels = scales::percent) +labs(title ="Changed forest ares because of drought (1=Yes, 2=No)?",x ="",y ="Percent" ) +theme_minimal()
Code
attr(data_2025_complete$A13, "label")
[1] "A13. How many trees in German forests do you estimate are currently affected by natural events such as drought, storms and bark beetle infestations?"
Code
attr(data_2025_complete$A13, "labels")
NULL
Code
ggplot(data=data_2025_complete) +geom_bar(aes(x =as.factor(A13)), position ="dodge") +labs(title ="How much trees are sick?", x ="", y ="Count") +theme_minimal()
Code
pal <-colorBin(palette =c("red", "yellow", "green"),bins =c(0.5, 2.5, 3.5, 5.5), # edges: (0.5–2.5], (2.5–3.5], (3.5–5.5]domain =c(1, 2, 3, 4, 5),right =TRUE# ensures e.g. 2.0 falls into red, 3.0 into yellow, 5.0 into green)leaflet(data_2025_complete) %>%addProviderTiles("OpenStreetMap") %>%addCircleMarkers(radius =5,fillColor =~pal(A2),fillOpacity =0.7,color ="#444444",weight =0.5, ) %>%addLegend("bottomright",pal = pal,values =~A2,title ="Drought Index (last visit)",opacity =0.7 )
library(tidyterra)points_sf <- data_2025_complete # assuming this is an sf objectpoints_sf <-st_transform(points_sf, crs =crs(forest_loss))points_in_raster_extent <- points_sf[st_within(points_sf, st_as_sfc(st_bbox(forest_loss)), sparse =FALSE), ]ggplot() +geom_spatraster(data = forest_loss_100) +scale_fill_viridis_c(na.value ="transparent", name ="Forest Loss", trans="sqrt") +geom_sf(data = points_in_raster_extent, aes(color = A3), size =1) +scale_color_viridis_c(name ="Drought Index", direction =-1, option ="A") +theme_minimal() +labs(title ="Forest Loss and Drought Index in Germany (100x100m resolution)") +coord_sf()
Code
ggplot() +geom_spatraster(data = forest_loss_1km) +scale_fill_viridis_c(na.value ="transparent", name ="Forest Loss", trans="sqrt") +geom_sf(data = points_in_raster_extent, aes(color = A3), size =1) +scale_color_viridis_c(name ="Drought Index", direction =-1, option ="A") +theme_minimal() +labs(title ="Forest Loss and Drought Index in Germany (1x1km resolution)") +coord_sf()
Code
ggplot(data=data_2025_complete) +geom_point(aes(x=forest_loss_visit_1km, y=dead_share_of_forest_5000, col=for_dis)) +labs(x="Dead forest share within 1km×1km cell of forest visit (%)", y="Share of dead forest within 5km radius around place of residence (%)")
Code
ggplot(data=data_2025_complete) +geom_point(aes(x=forest_loss_visit_1km, y=dead_share_of_forest_10000, col=for_dis)) +labs(x="Dead forest share within 1km×1km cell of forest visit (%)", y="Share of dead forest within 10km radius around place of residence (%)")
Code
ggplot(data=data_2025_complete) +geom_point(aes(x=forest_loss_visit_1km, y=dead_share_of_forest_15000, col=for_dis)) +labs(x="Dead forest share within 1km×1km cell of forest visit (%)", y="Share of dead forest within 15km radius around place of residence (%)")
Code
ggplot(data=data_2025_complete) +geom_point(aes(x=dead_share_of_forest_15000, y=dead_share_of_forest_5000, col=for_dis)) +labs(x="Share of dead forest within 15km radius around place of residence (%)", y="Share of dead forest within 5km radius around place of residence (%)")
Code
tab3 <- data_2025_complete %>%count(fl_lv_third, fl_por_5) %>%group_by(fl_lv_third) %>%mutate(share = n /sum(n))ggplot(tab3, aes(x = fl_por_5, y = fl_lv_third, fill = share)) +geom_tile() +geom_text(aes(label = scales::percent(share, accuracy =1)), color ="black", size =3) +scale_fill_viridis_c(name ="Share within residence\nforest-loss level", labels = scales::percent) +labs(x ="Forest loss (place of residence 5km radius)",y ="Forest loss (last visited forest 1kmx1km cell)",title ="Relative distribution of forest-loss categories" ) +theme_minimal(base_size =12)
Code
tab <- data_2025_complete %>%count(fl_lv_third, fl_por_10) %>%group_by(fl_lv_third) %>%mutate(share = n /sum(n))ggplot(tab, aes(x = fl_por_10, y = fl_lv_third, fill = share)) +geom_tile() +geom_text(aes(label = scales::percent(share, accuracy =1)), color ="black", size =3) +scale_fill_viridis_c(name ="Share within residence\nforest-loss level", labels = scales::percent) +labs(x ="Forest loss (place of residence 10km radius)",y ="Forest loss (last visited forest 1kmx1km cell)",title ="Relative distribution of forest-loss categories" ) +theme_minimal(base_size =12)
Code
tab2 <- data_2025_complete %>%count(fl_lv_third, fl_por_15) %>%group_by(fl_lv_third) %>%mutate(share = n /sum(n))ggplot(tab2, aes(x = fl_por_15, y = fl_lv_third, fill = share)) +geom_tile() +geom_text(aes(label = scales::percent(share, accuracy =1)), color ="black", size =3) +scale_fill_viridis_c(name ="Share within residence\nforest-loss level", labels = scales::percent) +labs(x ="Forest loss (place of residence 15km radius)",y ="Forest loss (last visited forest 1kmx1km cell)",title ="Relative distribution of forest-loss categories" ) +theme_minimal(base_size =12)
Code
tab_17 <- data_2017_complete %>%filter(!is.na(fl_lv_third), !is.na(fl_por_10)) %>%count(fl_lv_third, fl_por_10) %>%group_by(fl_lv_third) %>%mutate(share = n /sum(n))ggplot(tab_17, aes(x = fl_por_10, y = fl_lv_third, fill = share)) +geom_tile() +geom_text(aes(label = scales::percent(share, accuracy =1)), color ="black", size =3) +scale_fill_viridis_c(name ="Share within residence\nforest-loss level", labels = scales::percent) +labs(x ="Forest loss (place of residence 10km radius)",y ="Forest loss (last visited forest 1kmx1km cell)",title ="Relative distribution of forest-loss categories in 2017" ) +theme_minimal(base_size =12)
Code
tab_fil <- data_2025_complete %>%filter(for_dis <=4) %>%count(fl_lv_third, fl_por_5) %>%group_by(fl_lv_third) %>%mutate(share = n /sum(n))ggplot(tab_fil, aes(x = fl_por_5, y = fl_lv_third, fill = share)) +geom_tile() +geom_text(aes(label = scales::percent(share, accuracy =1)), color ="black", size =3) +scale_fill_viridis_c(name ="Share within residence\nforest-loss level", labels = scales::percent) +labs(x ="Forest loss (place of residence 5km radius)",y ="Forest loss (last visited forest 1kmx1km cell)",title ="Relative distribution of forest-loss categories (travel distance below 20km)" ) +theme_minimal(base_size =12)
ggplot(data_2025_complete_ger) +geom_sf(data = germany_fed, fill =NA, color ="grey50") +geom_sf(aes(geometry = geometry_depart, color =as.factor(fl_por_10)), size =1) +labs(color="Groups") +scale_color_manual(values =c("blue", "gray70", "red"),labels =c("low", "medium", "high"))
Code
tab31 <- database_2025 %>%count(fl_lv_third, alt3.deadwood) %>%group_by(fl_lv_third) %>%mutate(share = n /sum(n))ggplot(tab31, aes(x = alt3.deadwood, y = fl_lv_third, fill = share)) +geom_tile() +geom_text(aes(label = scales::percent(share, accuracy =1)), color ="black", size =3) +scale_fill_viridis_c(name ="Share", labels = scales::percent) +labs(x ="Deadwood in forest (reported)",y ="Forest loss (last visited forest 1kmx1km cell)",title ="Relative distribution of forest-loss categories" ) +theme_minimal(base_size =12)
Code
tab31 <- database_2025 %>%count(fl_lv_third, alt3.deadwood) %>%mutate(share = n /sum(n)) # no group_by before thisggplot(tab31, aes(x = alt3.deadwood, y = fl_lv_third, fill = share)) +geom_tile() +geom_text(aes(label = scales::percent(share, accuracy =1)), color ="black", size =3) +scale_fill_viridis_c(name ="Share", labels = scales::percent) +labs(x ="Deadwood in forest (reported)",y ="Forest loss (last visited forest 1kmx1km cell)",title ="Relative distribution of forest-loss categories (global 100%)" ) +theme_minimal(base_size =12)
Code
# ggplot() +# geom_violin(data = final_2017_merged,# aes(x = "2017", y = forest_loss1, fill = "2017"),# alpha = 0.6, color = "darkorange3") +# geom_violin(data = final_2025_merged,# aes(x = "2025", y = forest_loss, fill = "2025"),# alpha = 0.6, color = "darkgreen") +# scale_fill_manual(name = "Sample year",# values = c("2017" = "orange", "2025" = "forestgreen")) +# theme_minimal(base_size = 14) +# labs(title = "Distribution of (Hypothetical) Forest Loss (2017 vs 2025)",# x = "", y = "Forest Loss within 1km×1km Cell")
Code
points_vect <-vect(points_sf)# Step 2: Extract raster values at point locationsextracted <- terra::extract(forest_loss_1km, points_vect)# Step 3: Combine extracted values with original sf# `extracted` has ID + raster value, match back by rowpoints_sf$raster_value <- extracted[[2]] # assuming it's the 2nd columnggplot(points_sf, aes(x = A2, y = raster_value)) +geom_jitter(width =0.2, height =0, alpha =0.5, color ="darkgreen") +geom_smooth(method ="lm", se =TRUE, color ="red", linewidth =1) +labs(title ="Relationship Between Drought Index (A2) and Forest Loss",x ="Drought Index (A2: 1 = high drought, 5 = no drought)",y ="Forest Loss within 1km^2 Raster Cell (%)" ) +theme_minimal()
Code
ggsave("figures/forest_loss_perceived_drought.png", width=9, height=5, dpi="print")ggplot(points_sf, aes(x = A2, y = dead_share_of_forest_10000, col=fl_por_10)) +geom_jitter(width =0.2, height =0, alpha =0.5) +geom_smooth(method ="lm", se =TRUE, color ="red", linewidth =1) +labs(title ="Relationship Between Drought Index (A2) and Forest Loss around Place of Residence",x ="Drought Index (A2: 1 = high drought, 5 = no drought)",y ="Forest Loss within 10km radius (%)" ) +theme_minimal()
Code
ggplot(points_sf, aes(x = A2, y = share_forest_10000, col=fl_por_10)) +geom_jitter(width =0.2, height =0, alpha =0.5) +geom_smooth(method ="lm", se =TRUE, color ="red", linewidth =1) +labs(title ="Relationship Between Drought Index (A2) and Forest Share around Place of Residence",x ="Drought Index (A2: 1 = high drought, 5 = no drought)",y ="Forest Share within 10km radius (%)" ) +theme_minimal()
Code
points_sf$drought_cat <-cut( points_sf$A2,breaks =c(0, 2, 3, 5),labels =c("High", "Medium", "Low"),right =TRUE,include.lowest =TRUE)ggplot(points_sf, aes(x = drought_cat, y = raster_value)) +geom_boxplot(fill ="lightblue", outlier.shape =NA) +geom_jitter(width =0.2, alpha =0.4, color ="darkgreen") +labs(title ="Forest Loss by Drought Severity Category",x ="Drought Severity (A2 grouped)",y ="Forest Loss within 1 km² Raster Cell (%)" ) +theme_minimal()
Code
points_sf <- points_sf %>%mutate(drought_agg = A2 + A3)points_sf$drought_agg_cat <-cut( points_sf$drought_agg,breaks =c(0, 4, 7, 10),labels =c("High", "Medium", "Low"),right =TRUE,include.lowest =TRUE)points_sf <- points_sf %>%mutate(drought_agg = A2 + A3)ggplot(points_sf, aes(x = drought_agg_cat, y = raster_value)) +geom_boxplot(fill ="lightblue", outlier.shape =NA) +geom_jitter(width =0.2, alpha =0.4, color ="darkgreen") +labs(title ="Forest Loss by Drought Severity Category",x ="Drought Severity (aggregated A2 and A3)",y ="Forest Loss within 1 km² Raster Cell (%)" ) +theme_minimal()
Code
ggplot(points_sf, aes(x = drought_agg, y = raster_value)) +geom_jitter(width =0.2, height =0, alpha =0.5, color ="darkgreen") +geom_smooth(method ="lm", se =TRUE, color ="red", linewidth =1) +labs(title ="Relationship Between Drought Index and Forest Loss",x ="Drought Index (2 = high drought, 10 = no drought)",y ="Forest Loss within 1km^2 Raster Cell (%)" ) +theme_minimal()
Check spatial correlation in drought index
Code
data_sf <- data_2025_complete# Extract coordinates from the now-active geometrycoords <- sf::st_coordinates(data_sf)# Variable of interesty_raw <- data_sf$A2# Keep only rows with finite y and valid coordsok <-is.finite(y_raw) &is.finite(coords[,1]) &is.finite(coords[,2])coords_ok <- coords[ok, , drop =FALSE]y <- y_raw[ok]# =============================================# 1) Build distance-band neighbors (great-circle)# =============================================# We work in lon/lat with great-circle distances using longlat=TRUE.# Define distance breaks in *kilometers* (tune to your scale).breaks_km <-seq(0, 120, by =2) # 0-50, 50-100, ..., up to 1000 km# Helper to compute Moran's I with permutation p-value for one bandband_moran <-function(d1, d2) {# Build neighbors for distances in [d1, d2) nb <-try(dnearneigh(coords_ok, d1 = d1, d2 = d2, longlat =TRUE), silent =TRUE)if (inherits(nb, "try-error")) return(NULL)if (length(nb) ==0) return(NULL) # no neighbors at this band# Some bands may produce isolates; allow zero.policy lw <-nb2listw(nb, style ="W", zero.policy =TRUE)# Permutation test for Moran's I (more robust than asymptotic) mc <-moran.mc(y, lw, nsim =999, zero.policy =TRUE, alternative ="two.sided")# Mean neighbor distance (km) in this band for a nice x-position dlist <-nbdists(nb, coords_ok, longlat =TRUE) mean_d <-mean(unlist(dlist))tibble(d_from = d1,d_to = d2,d_mid =ifelse(is.finite(mean_d), mean_d, (d1 + d2)/2),moranI =as.numeric(mc$statistic),p_val = mc$p.value,n_links =sum(lengths(nb)) # total neighbor links used )}# Loop over bands and bind results (skip empty bands)res_list <-lapply(seq_len(length(breaks_km) -1), function(i) {band_moran(breaks_km[i], breaks_km[i +1])})cor_df <-bind_rows(res_list) %>%filter(!is.na(moranI)) %>%mutate(sig = p_val <0.05)# =============================# 2) Plot the distance correlogram# =============================ggplot(cor_df, aes(x = d_mid, y = moranI)) +geom_hline(yintercept =0, linetype ="dashed", linewidth =0.4, color ="grey40") +geom_line(linewidth =0.6) +geom_point(aes(shape = sig), size =2) +scale_shape_manual(values =c(`TRUE`=16, `FALSE`=1),labels =c(`TRUE`="p < 0.05", `FALSE`="ns")) +labs(title ="Distance correlogram for Perceived Drought (A2)",subtitle ="Moran's I across distance bands; filled points indicate p < 0.05",x ="Distance (km, band mid-point)",y ="Moran's I" ) +theme_minimal(base_size =12) +theme(legend.title =element_blank())
Code
#### do plot for dead forest share #####data_sf <- sf::st_set_geometry(data_2025_complete, "geometry_depart")# Extract coordinates from the now-active geometrycoords <- sf::st_coordinates(data_sf)# Variable of interesty_raw <- data_sf$forest_loss_visit_1km# Keep only rows with finite y and valid coordsok <-is.finite(y_raw) &is.finite(coords[,1]) &is.finite(coords[,2])coords_ok <- coords[ok, , drop =FALSE]y <- y_raw[ok]# =============================================# 1) Build distance-band neighbors (great-circle)# =============================================# We work in lon/lat with great-circle distances using longlat=TRUE.# Define distance breaks in *kilometers* (tune to your scale).breaks_km <-seq(0, 120, by =5) # 0-50, 50-100, ..., up to 1000 km# Helper to compute Moran's I with permutation p-value for one bandband_moran <-function(d1, d2) {# Build neighbors for distances in [d1, d2) nb <-try(dnearneigh(coords_ok, d1 = d1, d2 = d2, longlat =TRUE), silent =TRUE)if (inherits(nb, "try-error")) return(NULL)if (length(nb) ==0) return(NULL) # no neighbors at this band# Some bands may produce isolates; allow zero.policy lw <-nb2listw(nb, style ="W", zero.policy =TRUE)# Permutation test for Moran's I (more robust than asymptotic) mc <-moran.mc(y, lw, nsim =999, zero.policy =TRUE, alternative ="two.sided")# Mean neighbor distance (km) in this band for a nice x-position dlist <-nbdists(nb, coords_ok, longlat =TRUE) mean_d <-mean(unlist(dlist))tibble(d_from = d1,d_to = d2,d_mid =ifelse(is.finite(mean_d), mean_d, (d1 + d2)/2),moranI =as.numeric(mc$statistic),p_val = mc$p.value,n_links =sum(lengths(nb)) # total neighbor links used )}# Loop over bands and bind results (skip empty bands)res_list <-lapply(seq_len(length(breaks_km) -1), function(i) {band_moran(breaks_km[i], breaks_km[i +1])})cor_df <-bind_rows(res_list) %>%filter(!is.na(moranI)) %>%mutate(sig = p_val <0.05)# =============================# 2) Plot the distance correlogram# =============================ggplot(cor_df, aes(x = d_mid, y = moranI)) +geom_hline(yintercept =0, linetype ="dashed", linewidth =0.4, color ="grey40") +geom_line(linewidth =0.6) +geom_point(aes(shape = sig), size =2) +scale_shape_manual(values =c(`TRUE`=16, `FALSE`=1),labels =c(`TRUE`="p < 0.05", `FALSE`="ns")) +labs(title ="Distance correlogram for dead forest share in forest \nusing place of residence",subtitle ="Moran's I across distance bands; filled points indicate p < 0.05",x ="Distance (km, band mid-point)",y ="Moran's I" ) +theme_minimal(base_size =12) +theme(legend.title =element_blank())
Code
#### do plot for forest type #####data_sf <- data_2025_complete# Extract coordinates from the now-active geometrycoords <- sf::st_coordinates(data_sf)# Variable of interesty_raw <- data_sf$Q26# Keep only rows with finite y and valid coordsok <-is.finite(y_raw) &is.finite(coords[,1]) &is.finite(coords[,2])coords_ok <- coords[ok, , drop =FALSE]y <- y_raw[ok]# =============================================# 1) Build distance-band neighbors (great-circle)# =============================================# We work in lon/lat with great-circle distances using longlat=TRUE.# Define distance breaks in *kilometers* (tune to your scale).breaks_km <-seq(0, 150, by =5) # 0-50, 50-100, ..., up to 1000 km# Helper to compute Moran's I with permutation p-value for one bandband_moran <-function(d1, d2) {# Build neighbors for distances in [d1, d2) nb <-try(dnearneigh(coords_ok, d1 = d1, d2 = d2, longlat =TRUE), silent =TRUE)if (inherits(nb, "try-error")) return(NULL)if (length(nb) ==0) return(NULL) # no neighbors at this band# Some bands may produce isolates; allow zero.policy lw <-nb2listw(nb, style ="W", zero.policy =TRUE)# Permutation test for Moran's I (more robust than asymptotic) mc <-moran.mc(y, lw, nsim =999, zero.policy =TRUE, alternative ="two.sided")# Mean neighbor distance (km) in this band for a nice x-position dlist <-nbdists(nb, coords_ok, longlat =TRUE) mean_d <-mean(unlist(dlist))tibble(d_from = d1,d_to = d2,d_mid =ifelse(is.finite(mean_d), mean_d, (d1 + d2)/2),moranI =as.numeric(mc$statistic),p_val = mc$p.value,n_links =sum(lengths(nb)) # total neighbor links used )}# Loop over bands and bind results (skip empty bands)res_list <-lapply(seq_len(length(breaks_km) -1), function(i) {band_moran(breaks_km[i], breaks_km[i +1])})cor_df <-bind_rows(res_list) %>%filter(!is.na(moranI)) %>%mutate(sig = p_val <0.05)# =============================# 2) Plot the distance correlogram# =============================ggplot(cor_df, aes(x = d_mid, y = moranI)) +geom_hline(yintercept =0, linetype ="dashed", linewidth =0.4, color ="grey40") +geom_line(linewidth =0.6) +geom_point(aes(shape = sig), size =2) +scale_shape_manual(values =c(`TRUE`=16, `FALSE`=1),labels =c(`TRUE`="p < 0.05", `FALSE`="ns")) +labs(title ="Distance correlogram for type of forest",subtitle ="Moran's I across distance bands; filled points indicate p < 0.05",x ="Distance (km, band mid-point)",y ="Moran's I" ) +theme_minimal(base_size =12) +theme(legend.title =element_blank())
points_vect <-vect(points_sf)# Step 2: Extract raster values at point locationsextracted <- terra::extract(land_cover, points_vect)# Step 3: Combine extracted values with original sf# `extracted` has ID + raster value, match back by rowpoints_sf$land_cover <- extracted[[2]] # assuming it's the 2nd columntable(points_sf$land_cover, dnn="Land cover of marked forest location")
Land cover of marked forest location
NoData Artificial Land Open Soil
0 182 1
High Seasonal Vegetation High Perennial Vegetation Low Seasonal Vegetation
982 675 116
Low Perennial Vegetation Water
246 17
Code
# Check land cover for place of residence/departure pointpoints_custom_geom <-st_set_geometry(points_sf, points_sf$geometry_proj4)# Convert to SpatVectorpoints_vect2 <-vect(points_custom_geom)extracted_custom <- terra::extract(land_cover, points_vect2)points_sf$land_cover_depart <- extracted_custom[[2]] # assuming it's the 2nd columntable(points_sf$land_cover_depart, dnn="Land cover of departure point")
Land cover of departure point
NoData Artificial Land Open Soil
0 1313 2
High Seasonal Vegetation High Perennial Vegetation Low Seasonal Vegetation
337 155 164
Low Perennial Vegetation Water
256 29
Code
check_points <- points_sf %>%filter(land_cover =="Artificial Land")check_points <-st_transform(check_points, crs =3035) buffers <- terra::buffer(vect(check_points), width =100) # use terra::buffer for SpatVector# 2. Extract raster values from land_cover in these buffers# `cells = TRUE` gets the cell IDs; `touches = TRUE` includes intersecting cellsbuffer_extract <- terra::extract(land_cover, buffers, touches =TRUE)# 3. Join point IDs from check_points to the extracted buffer values# The 'ID' column in buffer_extract refers to row numbers of check_pointsbuffer_extract <-left_join(buffer_extract, check_points %>%mutate(ID =row_number()), by ="ID")# 4. For each point, check if at least one neighboring value is NOT "Artificial Land"buffer_summary <- buffer_extract %>%group_by(ID) %>%summarise(any_natural =any(category !="Artificial Land", na.rm =TRUE))# 5. Join back to flagged check_pointscheck_points <- check_points %>%mutate(ID =row_number()) %>%left_join(buffer_summary, by ="ID")table(check_points$any_natural, dnn="Natural areas in surrounding for weird patterns")
Natural areas in surrounding for weird patterns
FALSE TRUE
28 154
Code
check_points <-st_transform(check_points, crs =4326) # Ensure CRS is set for leafletcheck_points <- check_points %>%filter(any_natural =="FALSE")leaflet(check_points) %>%addProviderTiles("OpenStreetMap") %>%addCircleMarkers(radius =5)
Matching
Code
source("scripts/matching.R")love.plot( m.out_base,stats ="mean.diffs",abs =TRUE,var.order ="unadjusted",line =TRUE,thresholds =c(m = .1), # optional threshold line at |SMD| = 0.1title ="Base Matching all equal")
Code
love.plot( m.out_exact,stats ="mean.diffs",abs =TRUE,var.order ="unadjusted",line =TRUE,thresholds =c(m = .1), # optional threshold line at |SMD| = 0.1title ="Exact Matching Bundesland")
Code
love.plot( m.out_mah,stats ="mean.diffs",abs =TRUE,var.order ="unadjusted",line =TRUE,thresholds =c(m = .1), # optional threshold line at |SMD| = 0.1title ="Mahalanobis Matching Lat Lon")
b <-bal.tab(m.out_exact_county, un =TRUE)# drop all rows whose name starts with "NAME_2"b$Balance <- b$Balance[!grepl("^NAME_2", rownames(b$Balance)), ]love.plot( b,stats ="mean.diffs",abs =TRUE,var.order ="unadjusted",line =TRUE,thresholds =c(m = .1), # optional threshold line at |SMD| = 0.1title ="Exact Matching Landkreis")
