Code
paste0("You are seeing data from survey round ", unique(data_quarto_use$survey_round))[1] "You are seeing data from survey round 1"Overview Pre-test Data 1
Test Link
Current survey version: https://s.onfly.pl/Forestde/
Data Version
Basic checks
Code
print("Look at only NA columns")[1] "Look at only NA columns"Code
data_complete %>%
select(where(~ all(is.na(.)))) %>%
names() [1] "JEZYK" "M1_PANEL_1" "M2_PANEL"
[4] "M6_PANEL" "M7_PANEL" "TEST"
[7] "ZRODLO" "N10_5_other" "N24_2_none"
[10] "QCombined_podszyt" "QCombined_zreb" "QCombined_pow"
[13] "Q36_2_4_other" "Q36_3_4_other" "Q36_4_4_other"
[16] "Q36_5_4_other" "EIG4_6_other" Code
summary(data_complete$duration_min) Min. 1st Qu. Median Mean 3rd Qu. Max.
6.679 13.392 18.213 26.576 22.932 737.275 Code
ggplot(data = data_complete) +
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_pre$M3v, "labels")18-24 25-34 35-44 45-54 55-64 65+
1 2 3 4 5 6 Code
ggplot(data = data_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_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 == 2
table(data_pre$N1, useNA = "ifany") # 167 drop outs, 161 people surveyed
1 2 <NA>
157 167 5 Code
sum(table(data_pre$N1, useNA = "ifany")) # 333 respondents overall[1] 329Code
# 2. NCONTROL (before DCE) =! "Stimme teilweise zu"
table(data_pre$NCONTROL, useNA = "ifany") # no further drop outs
2 <NA>
157 172 Code
# 3. Survey version on drought: A10_99 =! "trifft eher zu" (4)
attr(data_pre$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_pre$A10_99, "labels")1 strongly disagree 2 partly disagree 3 neutral 4 partly agree
1 2 3 4
5 strongly agree
5 Code
table(data_pre$A10_99, useNA = "ifany") # no further drop outs (80 (of 161) people in drought version of the survey)
4 <NA>
78 251 Code
# check how many people answer the first drought question (asked in both survey versions)
attr(data_pre$A1, "label")[1] "A1. Have the forest areas that you regularly visit changed in the last 10 years?"Code
attr(data_pre$A1, "labels") No Yes, the condition has improved
1 2
Yes, the condition has deteriorated I can't say
3 4 Code
table(data_pre$A1, useNA = "ifany") # all 161 people
1 2 3 4 <NA>
41 21 61 34 172 Code
sum(table(data_pre$A1))[1] 157Code
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_pre$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_pre$A4_1, "labels") 1 Major Problem 2 Rather a problem 3 Neutral
1 2 3
4 Rather an opportunity 5 Major Opportunity
4 5 Code
table(data_pre$A4_1, useNA = "ifany") # 80 people left in survey version on drought
1 2 3 4 5 <NA>
41 17 17 1 2 251 Code
sum(table(data_pre$A4_1))[1] 78Code
# check how many answered the question after the last control question
attr(data_pre$A11, "label")[1] "A11. Are you concerned about the long-term availability of healthy forests in Germany?"Code
attr(data_pre$A11, "labels")Yes, very strongly Yes, a little I can’t judge No, hardly
1 2 3 4
No, not at all
5 Code
table(data_pre$A11, useNA = "ifany")
1 2 3 4 5 <NA>
26 33 9 7 3 251 Code
sum(table(data_pre$A11)) # still all 80 in survey after last control question[1] 78Code
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_complete$N2, useNA = "ifany") # 126 Waldbesuch einziges Ziel, 35 respondents Waldbesuch mit anderen Aktivitäten
1 2
122 35 Code
block_N3_to_N10 <- c("N3_1", "N3_2", "N4_1", "N4_2", "N4_9_none", "N5",
paste0("N6_", 1:17), "N6_17_other", "N7_1", "N8_1", "N8_2", "N9", "N10", "N10_5_other")
block_N11_to_N22 <- c("N11", "N11_7_other", "N12_1_1", "N13_1", "N13_2", "N14_1", "N14_2", "N14_9_none", "N15",
paste0("N16_", 1:17), "N16_17_other", "N17_1", "N18_1", "N18_2", "N19_1", "N19_2", "N19_9_none",
"N20_1", "N21", "N22", "N22_5_other")
# Create Flags for Response in Each Block (TRUE/FALSE)
data_complete <- data_complete %>%
mutate(
block1_answered = rowSums(select(., all_of(block_N3_to_N10)) != "" & !is.na(select(., all_of(block_N3_to_N10)))) > 0,
block2_answered = rowSums(select(., all_of(block_N11_to_N22)) != "" & !is.na(select(., all_of(block_N11_to_N22)))) > 0
)
# Check Consistency with N2 Routing
table(data_complete$N2, data_complete$block1_answered, data_complete$block2_answered, useNA = "ifany"), , = FALSE
FALSE TRUE
1 0 122
2 0 0
, , = TRUE
FALSE TRUE
1 0 0
2 35 0Code
count(data_complete, N2, block1_answered, block2_answered)# A tibble: 2 × 4
N2 block1_answered block2_answered n
<dbl> <lgl> <lgl> <int>
1 1 TRUE FALSE 122
2 2 FALSE TRUE 35Code
# delete helping vars
data_complete$block_N3_to_N10 <- NULL
data_complete$block_N11_to_N22 <- NULL
data_complete$block1_answered <- NULL
data_complete$block2_answered <- NULL
print("2. Single, smaller conditions:")[1] "2. Single, smaller conditions:"Code
print("2.1: N23: Only shown if N10 or N22 == 3")[1] "2.1: N23: Only shown if N10 or N22 == 3"Code
table(data_complete$N10, useNA = "ifany") # N10 == 3
1 2 3 4 <NA>
44 14 5 59 35 Code
table(data_complete$N22, useNA = "ifany") # N22 == 3
1 2 3 4 5 <NA>
9 3 6 16 1 122 Code
table(data_complete$N23[!(data_complete$N10 == 3 | data_complete$N22 == 3)], useNA = "ifany") # only NA
<NA>
146 Code
count(data_complete, N2, N23, N10, N22)# A tibble: 10 × 5
N2 N23 N10 N22 n
<dbl> <dbl> <dbl> <dbl> <int>
1 1 1 3 NA 4
2 1 2 3 NA 1
3 1 NA 1 NA 44
4 1 NA 2 NA 14
5 1 NA 4 NA 59
6 2 1 NA 3 6
7 2 NA NA 1 9
8 2 NA NA 2 3
9 2 NA NA 4 16
10 2 NA NA 5 1Code
print("2.2: N23a: Only if N23 == 2 (Nein)")[1] "2.2: N23a: Only if N23 == 2 (Nein)"Code
table(data_complete$N23a_1[data_complete$N23 != 2], useNA = "ifany")
<NA>
156 Code
count(data_complete, N2, N23, N23a_1)# A tibble: 5 × 4
N2 N23 N23a_1 n
<dbl> <dbl> <dbl> <int>
1 1 1 NA 4
2 1 2 4 1
3 1 NA NA 117
4 2 1 NA 6
5 2 NA NA 29Code
print("2.3: N24: Only if N10 or N22 == 4 (car)")[1] "2.3: N24: Only if N10 or N22 == 4 (car)"Code
table(data_complete$N24_1[!(data_complete$N10 == 4 | data_complete$N22 == 4)], useNA = "ifany")
<NA>
82 Code
count(data_complete, N2, N24_1, N10, N22)# A tibble: 15 × 5
N2 N24_1 N10 N22 n
<dbl> <dbl> <dbl> <dbl> <int>
1 1 1 4 NA 37
2 1 2 4 NA 14
3 1 3 4 NA 7
4 1 4 4 NA 1
5 1 NA 1 NA 44
6 1 NA 2 NA 14
7 1 NA 3 NA 5
8 2 1 NA 4 11
9 2 2 NA 4 2
10 2 3 NA 4 2
11 2 4 NA 4 1
12 2 NA NA 1 9
13 2 NA NA 2 3
14 2 NA NA 3 6
15 2 NA NA 5 1Code
print("2.4: N25: only if N10 or N22 == 4")[1] "2.4: N25: only if N10 or N22 == 4"Code
table(data_complete$N25[!(data_complete$N10 == 4 | data_complete$N22 == 4)], useNA = "ifany")
<NA>
82 Code
count(data_complete, N2, N25, N25_1_other, N10, N22)# A tibble: 36 × 6
N2 N25 N25_1_other N10 N22 n
<dbl> <dbl> <chr> <dbl> <dbl> <int>
1 1 1 1 4 NA 1
2 1 1 10 4 NA 1
3 1 1 11 4 NA 1
4 1 1 12 4 NA 1
5 1 1 13 4 NA 1
6 1 1 4,2 4 NA 1
7 1 1 4;5 4 NA 1
8 1 1 5 4 NA 3
9 1 1 5,0 4 NA 1
10 1 1 5,3 4 NA 1
# ℹ 26 more rowsCode
print("2.5: N26: Only if N25 == 2 (Nein)")[1] "2.5: N26: Only if N25 == 2 (Nein)"Code
table(data_complete$N26[data_complete$N25 != 2], useNA = "ifany")
<NA>
127 Code
count(data_complete, N2, N26, N25)# A tibble: 12 × 4
N2 N26 N25 n
<dbl> <dbl> <dbl> <int>
1 1 2 2 2
2 1 3 2 7
3 1 4 2 12
4 1 5 2 2
5 1 6 2 1
6 1 NA 1 35
7 1 NA NA 63
8 2 2 2 1
9 2 3 2 4
10 2 5 2 1
11 2 NA 1 10
12 2 NA NA 19Code
print("2.6: N29–N32 if N28 == 1")[1] "2.6: N29–N32 if N28 == 1"Code
attr(data_pre$N28, "label")[1] "N28. Was your overall trip that included the forest visit longer or shorter than 1 day?"Code
attr(data_pre$N28, "labels")Longer than one day 1 day or shorter
1 2 Code
table(data_complete$N28, useNA = "ifany")
1 2
24 133 Code
# Show who answered when they shouldn't have
table(data_complete$N29_1[data_complete$N28 != 1], useNA = "ifany")
<NA>
133 Code
table(data_complete$N30[data_complete$N28 != 1], useNA = "ifany")
<NA>
133 Code
table(data_complete$N31[data_complete$N28 != 1], useNA = "ifany")
<NA>
133 Code
table(data_complete$N32_1[data_complete$N28 != 1], useNA = "ifany")
<NA>
133 Code
table(data_complete$N32_2_none[data_complete$N28 != 1], useNA = "ifany")
<NA>
133 Code
print("Check C: DCE-PART")[1] "Check C: DCE-PART"Code
print("Q26 Routing Logic")[1] "Q26 Routing Logic"Code
# Q26=1 -> trigger Q26a
# Q26=2 -> trigger Q26b
# Q26=3 -> trigger Q26c
# Count inconsistencies
data_complete %>%
filter((Q26 == 1 & is.na(Q26a)) |
(Q26 == 2 & is.na(Q26b)) |
(Q26 == 3 & is.na(Q26c))) %>%
nrow()[1] 0Code
count(data_complete, Q26, Q26a, Q26b, Q26c) # A tibble: 11 × 5
Q26 Q26a Q26b Q26c n
<dbl> <dbl> <dbl> <dbl> <int>
1 1 1 NA NA 5
2 1 2 NA NA 4
3 1 3 NA NA 6
4 2 NA 1 NA 6
5 2 NA 2 NA 3
6 2 NA 3 NA 12
7 3 NA NA 1 7
8 3 NA NA 2 15
9 3 NA NA 3 40
10 3 NA NA 4 9
11 3 NA NA 5 50Code
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 consistency
data_complete %>%
filter((Q28 %in% c(1,2) & is.na(Q29)) |
(Q28 == 3 & is.na(Q29a))) %>%
nrow()[1] 0Code
# test for false positives
data_complete %>%
filter((Q28 %in% c(1,2) & !is.na(Q29a)) |
(Q28 == 3 & !is.na(Q29))) %>%
nrow()[1] 0Code
count(data_complete, Q28, Q29, Q29a) # A tibble: 8 × 4
Q28 Q29 Q29a n
<dbl> <dbl> <dbl> <int>
1 1 1 NA 2
2 1 3 NA 1
3 2 1 NA 8
4 2 2 NA 20
5 2 3 NA 35
6 3 NA 1 11
7 3 NA 2 15
8 3 NA 3 65Code
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 answered
data_complete %>%
filter(is.na(Q31_1)) %>%
nrow()[1] 0Code
data_complete %>%
filter(Q31_1 < 5 & is.na(Q31a_1)) %>%
nrow()[1] 0Code
count(data_complete, Q31_1, Q31a_1)# A tibble: 12 × 3
Q31_1 Q31a_1 n
<dbl> <chr> <int>
1 2 "Wald wächst nicht in gerader Linie" 1
2 3 "ca. 80 % des Bestandes sieht schlecht aus" 1
3 3 "ich kenne mich damit nicht aus, aber das ist der Eindruck den i… 1
4 3 "viele sturmschäden" 1
5 4 "abholzung" 1
6 4 "braun" 1
7 5 "" 8
8 6 "" 11
9 7 "" 38
10 8 "" 47
11 9 "" 23
12 10 "" 24Code
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 interesting
test <- data_complete %>%
filter(Q28 == 3 & Q30 %in% c(1, 2)) # 42 obs
rm(test)
print("Other considered attributes: Q36 show if Q35_1a==2")[1] "Other considered attributes: Q36 show if Q35_1a==2"Code
attr(data_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_complete$Q35_1a, useNA = "ifany") # does not indicate anything
0
157 Code
attr(data_pre$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_complete$Q35_97, useNA = "ifany") # 62 respondents overlooked none of the attributes
0 1
97 60 Code
#View(data_complete[, c("Q35_1a", "Q35_1", "Q35_2", "Q35_3", "Q35_4", "Q35_5", "Q35_6", "Q35_7", "Q35_97")])
data_complete %>%
filter(Q35_1 == 1 |
Q35_2 == 1 |
Q35_3 == 1 |
Q35_4 == 1 |
Q35_5 == 1 |
Q35_6 == 1 |
Q35_7 == 1) %>%
nrow() # 99[1] 97Code
table(data_complete$Q35_1a, data_complete$Q35_97, useNA = "ifany")
0 1
0 97 60Code
count(data_complete, Q35_97, Q36_1, Q36_2, Q36_3, Q36_4, Q36_5, Q36_6, Q36_7)# A tibble: 60 × 9
Q35_97 Q36_1 Q36_2 Q36_3 Q36_4 Q36_5 Q36_6 Q36_7 n
<dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <int>
1 0 1 1 1 1 1 1 NA 2
2 0 1 1 1 1 1 NA NA 4
3 0 1 1 1 NA NA 1 NA 1
4 0 1 1 NA NA 1 1 NA 1
5 0 1 1 NA NA NA 1 1 1
6 0 1 1 NA NA NA NA NA 1
7 0 1 2 NA NA NA NA NA 1
8 0 1 NA 1 1 1 NA 4 1
9 0 1 NA 1 1 NA 1 NA 1
10 0 1 NA 3 NA NA 4 NA 1
# ℹ 50 more rowsCode
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 questions
attr(data_complete$Q38, "labels") NULLCode
table(data_complete$Q38)
1 2
63 94 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_complete %>% filter(Q38 == 2 & !is.na(Q39)) %>% nrow()[1] 0Code
data_complete %>% filter(Q38 == 2 & is.na(Q39)) %>% nrow()[1] 94Code
# 2. Q40a–Q42a only if Q39 == 2
data_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] 0Code
# 3. Q41b–Q43b only if Q39 in 3:6
data_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] 0Code
# 4. Optional: Check if expected data is missing when it should be present
# 4.1 Q39 == 2 --> Check for missing responses in forest
data_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(.))
)
)# A tibble: 1 × 7
Q40a_1 Q40a_2 Q41a_1 Q41a_2 Q41a_9_none Q41ax Q42a_1
<int> <int> <int> <int> <int> <int> <int>
1 0 0 1 1 13 0 0Code
# 4.2 Q39 ∈ 3:6 --> Check for missing responses in forest
data_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() # no missing values[1] 0Code
print("Berries & Mushrooms")[1] "Berries & Mushrooms"Code
# 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_complete$Q45)
1 2
37 120 Code
data_complete %>%
filter(Q45 != 1 & Q46_1 != "" & !is.na(Q46_1)) %>%
nrow()[1] 0Code
#View(data_complete[,c("Q45", "Q46_1")])
table(data_complete$Q47)
1 2
24 133 Code
# Should NOT be answered if Q47 != 1
data_complete %>%
filter(Q47 != 1 & Q48_1 != "" & !is.na(Q48_1)) %>%
nrow()[1] 0Code
# Should be answered if Q47 == 1
data_complete %>%
filter(Q47 == 1 & Q48_1 == "") %>%
nrow()[1] 0Code
#View(data_complete[,c("Q47", "Q48_1")])Code
print("Check D: DEMOGRAPHY")[1] "Check D: DEMOGRAPHY"Code
print("M11: M10 >= 2017")[1] "M11: M10 >= 2017"Code
data_complete %>%
filter((M10_1 < 2017) & !is.na(M11)) %>%
nrow()[1] 0Code
#View(data_complete[,c("M10_1","M10_2_none","M11")])Code
print("Check E: Ownership")[1] "Check E: Ownership"Code
print("Trigger several questions if EIG1==1")[1] "Trigger several questions if EIG1==1"Code
data_complete %>%
filter((EIG1 !=1) & (!is.na(EIG2_1) & !is.na(EIG2_2_none))) %>%
nrow()[1] 0Code
data_complete %>%
filter((EIG1 ==1) & (is.na(EIG2_1) & is.na(EIG2_2_none))) %>%
nrow()[1] 0Code
data_complete %>%
filter((EIG1 !=1) & (!is.na(EIG3_1) & !is.na(EIG3_2_none))) %>%
nrow()[1] 0Code
#View(data_complete[,c("EIG1","EIG2_2_none","EIG2_1","EIG3_2_none","EIG3_1")])Code
print("Check F: Drought") # 80 respondents in this survey split (81 on other survey split)[1] "Check F: Drought"Code
print("Trigger several questions only if A2 AND A3 != “No, not at all”")[1] "Trigger several questions only if A2 AND A3 != “No, not at all”"Code
table(data_complete$A2, useNA = "ifany")
1 2 3 4 5
26 51 36 28 16 Code
table(data_complete$A3, useNA = "ifany")
1 2 3 4 5
23 54 44 25 11 Code
table(data_complete$A2, data_complete$A3, useNA = "ifany")
1 2 3 4 5
1 17 9 0 0 0
2 6 34 9 2 0
3 0 8 26 2 0
4 0 2 7 16 3
5 0 1 2 5 8Code
# data_complete %>% filter(A2 == 5 & A3 == 5) %>%
# select(A2, A3, A5, A6, A8) %>%
# View()
#%>% nrow() # 8 times the case that questions are not triggered
sum(!(data_complete$A2 == 5 & data_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?[1] 149Code
data_complete %>%
filter((A2 ==5 & A3==5) & (!is.na(A5) & !is.na(A6) & !is.na(A8))) %>%
nrow()[1] 0Code
data_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] 72Code
table(data_complete$A5, useNA = "ifany") # answered by 74, NA:87
1 2 3 4 5 6 <NA>
4 15 40 1 3 9 85 Code
table(data_complete$A6, useNA = "ifany") # answered by 74, NA:87
2 3 4 5 6 <NA>
7 43 4 2 16 85 Code
table(data_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>
11 55 91 Code
table(data_complete$A9_1, useNA = "ifany")
149
abwechslung beim laufen
1
abwechsung
1
bunker wald
1
Ich habe meine Waldbesuche auf andere Gebiete verlegt, weil ich neue Wege und Landschaften entdecken wollte und mir Abwechslung wichtig ist. Außerdem waren manche der früher besuchten Wälder stärker frequentiert, was die Ruhe und Erholung beeinträchtigt hat. I
1
Ich will mehr andere Wälde sehen.
1
mehr kennen lernen
1
nein
1
Neugier
1 Code
# View(data_complete[,c("A2","A3","A5","A6","A8")])
print("Trigger A7 only if A6==1|2|4|5")[1] "Trigger A7 only if A6==1|2|4|5"Code
data_complete %>%
filter((A6 == 3 | A6 == 6) & A7_1 != "") %>%
nrow()[1] 0Code
print("Trigger A9 only if A8==1")[1] "Trigger A9 only if A8==1"Code
data_complete %>%
filter(A8 == 2 & A9_1 != "") %>%
nrow()[1] 0Spatial Distribution of Forest Visits
Code
pal <- colorNumeric(
palette = "viridis",
domain = log1p(data_complete_sf$num_visits)
)
leaflet(data_complete_sf) %>%
addProviderTiles("OpenStreetMap") %>%
addCircleMarkers(
radius = 5,
fillColor = ~pal(log1p(num_visits)), # log1p here
fillOpacity = 0.7,
color = "#444444",
weight = 0.5,
popup = ~paste("Forest Visits:", num_visits)
) %>%
addLegend(
"bottomright",
pal = pal,
values = ~log1p(num_visits), # log1p here
title = "Forest Visits (log scale)",
opacity = 0.7,
labFormat = labelFormat(transform = function(x) round(expm1(x))) # back-transform legend labels
)Code
print("Check distance to nearby village location")[1] "Check distance to nearby village location"Code
summary(data_complete_sf_check$dist_check/1000) Min. 1st Qu. Median Mean 3rd Qu. Max.
0.02049 2.02034 4.99903 20.26536 9.52247 289.89298 Code
ggplot(data = data_complete_sf_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_complete_sf$geometry3 <- st_sfc(
Map(function(lon, lat) st_point(c(lon, lat)),
data_complete_sf$N8_2, # longitude
data_complete_sf$N8_1), # latitude
crs = 4326
)
leaflet(data_complete_sf) %>%
addTiles() %>% # Adds OpenStreetMap base layer
addMarkers(data = st_set_geometry(data_complete_sf, "geometry3"))Code
data_complete_sf <- data_complete_sf %>%
mutate(
geometry_proj1 = st_transform(geometry, 3857),
geometry_proj3 = st_transform(geometry3, 3857)
)
data_complete_sf$dist_por_for <- st_distance(
data_complete_sf$geometry_proj1,
data_complete_sf$geometry_proj3,
by_element = TRUE
)
print("Distance place of residence to forest")[1] "Distance place of residence to forest"Code
summary(data_complete_sf$dist_por_for/1000) Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
0.0194 3.3612 8.5461 49.0327 22.4418 1172.8614 35 Socio-demographic Overview
Code
ggplot(data=data_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_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_complete$M3v)
1 2 3 4 5 6
1 11 16 20 48 61 Code
print("Gender")[1] "Gender"Code
table(data_complete$M1)
1 2
88 69 Code
attr(data_pre$M2, "labels")Rural area or settlement with less than 500 inhabitants
1
Rural area or settlement 500 - 2.999 inhabitants
2
A small town with 3.000 – 19.999 inhabitants
3
A town with 20.000 - 100.000 inhabitants
4
More than 100.000 inhabitants
5 Code
ggplot(data = data_complete) +
geom_bar(aes(x = M2), position = "dodge") +
labs(title="Settlement type")
Code
attr(data_pre$M4, "labels") Primary school Secondary school Vocational training
1 2 3
Undergraduate studies Postgraduate studies
4 5 Code
ggplot(data = data_complete) +
geom_bar(aes(x = M4), position = "dodge") +
labs(title="Education")
Code
ggplot(data=data_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_complete) +
geom_bar(aes(x = as.factor(M9)), position = "dodge") +
labs(title = "Net Income", x = "", y = "Count") +
theme_minimal() +
scale_fill_brewer(palette = "Set3")
Code
ggplot(data=data_complete) +
geom_boxplot(aes(x="", y=M5_1)) +
labs(title="HH Size")
Code
ggplot(data=data_complete) +
geom_bar(aes(x=M5_1)) +
labs(title="HH Size")
Code
ggplot(data=data_complete) +
geom_bar(aes(x=M7_1)) +
labs(title="Kids in HH")
Code
ggplot(data=data_complete) +
geom_bar(aes(x=M10_1)) +
labs(title="Living in Region since")
Forest Visit Stats
Code
attr(data_pre$M2a, "labels") Baden-Württemberg Bayern Berlin
1 2 3
Brandenburg Bremen Hamburg
4 5 6
Hessen Mecklenburg-Vorpommern Niedersachsen
7 8 9
Nordrhein-Westfalen Rheinland-Pfalz Saarland
10 11 12
Sachsen Sachsen-Anhalt Schleswig-Holstein
13 14 15
Thüringen
16 Code
ggplot(data = data_complete) +
geom_boxplot(aes(x = factor(M2a), y = N7_1, 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)) 
Spatial Forest Characteristics
Code
pal <- colorFactor(
palette = c("green", "yellow", "red"),
domain = c(1, 2, 3)
)
leaflet(data_complete_sf) %>%
addProviderTiles("OpenStreetMap") %>%
addCircleMarkers(
radius = 5,
fillColor = ~pal(Q30),
fillOpacity = 0.7,
color = "#444444",
weight = 0.5,
popup = ~paste("Deadwood Indicator:", Q30)
) %>%
addLegend(
"bottomright",
pal = pal,
values = ~Q30,
title = "Deadwood Indicator",
opacity = 0.7,
labFormat = labelFormat(
transform = identity,
digits = 0
)
)Code
forest_labels <- c("1" = "Coniferous", "2" = "Broadleaved", "3" = "Mixed")
pal_q26 <- colorFactor(
palette = c("darkgreen", "saddlebrown", "orange"),
domain = c(1, 2, 3)
)
leaflet(data_complete_sf) %>%
addProviderTiles("OpenStreetMap") %>%
addCircleMarkers(
radius = 5,
fillColor = ~pal_q26(Q26),
fillOpacity = 0.7,
color = "#444444",
weight = 0.5,
popup = ~paste("Forest Type:", forest_labels[as.character(Q26)])
) %>%
addLegend(
"bottomright",
pal = pal_q26,
values = ~Q26,
title = "Forest Type",
opacity = 0.7,
labels = c("1" = "Coniferous", "2" = "Broadleaved", "3" = "Mixed")
)Code
attr(data_pre$N5, "labels") Less than 30 min 31-60 min 1-2 hours 2-3 hours
1 2 3 4
3-5 hours 5-8 hours more than 8 hours
5 6 7 Code
ggplot(data=data_complete) +
geom_bar(aes(x = N5), position = "dodge") +
labs(title = "Time Spend in Forest", x = "", y = "Count") +
theme_minimal() 
Code
attr(data_pre$N9, "labels") Less than 1 km 1-4 km 5-9 km 10-19 km
1 2 3 4
20-39 km 40-69 km 70-99 km 100-150 km
5 6 7 8
More than 150 km
9 Code
breaks_km <- c(-Inf, 1, 4, 9, 19, 39, 69, 99, 150, Inf)
labels_n9 <- c(
"Less than 1 km", "1-4 km", "5-9 km", "10-19 km",
"20-39 km", "40-69 km", "70-99 km", "100-150 km", "More than 150 km"
)
data_complete$for_dis_labelled <- factor(data_complete$for_dis, levels = 1:9, labels = labels_n9)
ggplot(data = data_complete) +
geom_bar(aes(x = for_dis_labelled), position = "dodge", fill = "steelblue") +
labs(title = "Distance to Forest", x = "Distance Category", y = "Count") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
Code
data_complete_sf$dist_bin <- cut(
data_complete_sf$dist_por_for/ 1000, # convert to km
breaks = breaks_km,
labels = labels_n9,
right = TRUE, # interval is (a, b]
include.lowest = TRUE
)
true_counts <- data_complete_sf %>%
mutate(bin = factor(for_dis, levels = 1:9, labels = labels_n9)) %>%
count(bin, source = "Survey")
# For computed distances (assuming already labeled with same strings)
comp_counts <- data_complete_sf %>%
filter(!is.na(dist_bin)) %>%
mutate(bin = factor(dist_bin, levels = labels_n9)) %>%
count(bin, source = "Computed")
# Combine and align factor levels
combined_counts <- bind_rows(true_counts, comp_counts) %>%
mutate(bin = factor(bin, levels = labels_n9))
ggplot(combined_counts, aes(x = bin, y = n, fill = source)) +
geom_bar(stat = "identity",
position = "identity",
alpha = 0.4,
data = subset(combined_counts, source == "Survey")) +
geom_bar(stat = "identity",
position = "identity",
alpha = 1, width = 0.5,
data = subset(combined_counts, source == "Computed")) +
scale_fill_manual(values = c("Survey" = "grey70", "Computed" = "indianred")) +
theme_minimal() +
labs(title = "Distance to Forest: Reported vs Computed",
x = "Distance Category",
y = "Count",
fill = "") +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
Choice Experiment
Code
attr(data_pre$Q26, "labels")Coniferous composed mostly of coniferous tree e.g. Pine, Spruce, Fir
1
Broadleaved composed mostly of broadleved tree e.g. Oak, Beech, Ash
2
Mix of both composed of coniferous and broadleved tree species
3 Code
ggplot(data=data_complete) +
geom_bar(aes(x=Q26), fill="forestgreen") +
labs(title="Forest Type")
Code
attr(data_pre$Q28, "labels")Recently planted Height around 8 meters Growing Height around 18 meters
1 2
Mature Height around 24 meters or more
3 Code
ggplot(data=data_complete) +
geom_bar(aes(x=Q28), fill="grey30") +
labs(title="Tree Height")
Code
attr(data_pre$Q29, "labels")Single-aged composed of trees are of the same age and similar size
1
Two-aged composed of trees that are of two age and size classes
2
Multi-aged composed of trees of varying age and size classes
3 Code
ggplot(data=data_complete) +
geom_bar(aes(x=Q29), fill="darkorange") +
labs(title="Age Composition")
Code
attr(data_pre$Q30, "labels") Keine Keine Bäume, die dem natürlichen Verfall überlassen werden
1
Wenige Wenige Bäume, die dem natürlichen Verfall überlassen werden; im Durchschnitt finden Sie alle 50 m einen solchen Baum
2
Medium Einige Bäume, die dem natürlichen Verfall überlassen werden; im Durchschnitt finden Sie alle 25 m einen solchen Baum
3 Code
ggplot(data=data_complete) +
geom_bar(aes(x=Q30), fill="cyan") +
labs(title="Deadwood")
Code
ggplot(data=data_complete) +
geom_bar(aes(x=as.factor(Q31_1))) +
labs(title="How well is forest displayed?")
Code
print("check randomization of blocks")[1] "check randomization of blocks"Code
table(data_complete$CJ)
1 2 3
54 52 51 Choices
Code
print("First choice")[1] "First choice"Code
table(data_complete$CJ_1_1)
1 2 3
21 31 105 Code
table(data_complete$CJ_2_1)
1 2 3
24 26 107 Code
table(data_complete$CJ_3_1)
1 2 3
18 17 122 Code
table(data_complete$CJ_4_1)
1 2 3
23 18 116 Code
table(data_complete$CJ_5_1)
1 2 3
10 29 118 Code
table(data_complete$CJ_6_1)
1 2 3
17 24 116 Code
table(data_complete$CJ_7_1)
1 2 3
17 24 116 Code
table(data_complete$CJ_8_1)
1 2 3
43 11 103 Code
table(data_complete$CJ_9_1)
1 2 3
27 22 108 Code
table(data_complete$CJ_10_1)
1 2 3
21 25 111 Code
table(data_complete$CJ_11_1)
1 2 3
10 34 113 Code
table(data_complete$CJ_12_1)
1 2 3
22 14 121 Code
print("Second choice")[1] "Second choice"Code
table(data_complete$CJ_1_2)
1 2 3
60 59 38 Code
table(data_complete$CJ_2_2)
1 2 3
62 56 39 Code
table(data_complete$CJ_3_2)
1 2 3
59 73 25 Code
table(data_complete$CJ_4_2)
1 2 3
77 47 33 Code
table(data_complete$CJ_5_2)
1 2 3
23 104 30 Code
table(data_complete$CJ_6_2)
1 2 3
53 72 32 Code
table(data_complete$CJ_7_2)
1 2 3
54 70 33 Code
table(data_complete$CJ_8_2)
1 2 3
78 36 43 Code
table(data_complete$CJ_9_2)
1 2 3
62 56 39 Code
table(data_complete$CJ_10_2)
1 2 3
77 42 38 Code
table(data_complete$CJ_11_2)
1 2 3
48 74 35 Code
table(data_complete$CJ_12_2)
1 2 3
66 61 30 Check design with 2017 version
Code
print("Distance levels 2017")[1] "Distance levels 2017"Code
table(database_2017$dist_1)
0.5 1 2 3 4 6 8 9 10 12 15 18 25 30 32 40
1276 857 634 1083 621 1519 467 190 490 1502 204 481 684 203 434 348
45 50 60 65 90 125
116 210 100 195 194 216 Code
print("Distance levels 2025")[1] "Distance levels 2025"Code
table(database_2025$alt1.dist_trans)
1 2 4 5 6 10 12 15 20 25 30 40 50 60 80 100 120
229 141 87 176 93 254 67 90 224 36 93 153 102 98 23 8 10 Code
prop.table(table(database_2017$dist_1))
0.5 1 2 3 4 6
0.106121091 0.071274118 0.052727878 0.090069860 0.051646707 0.126330672
8 9 10 12 15 18
0.038838989 0.015801730 0.040751830 0.124916833 0.016966068 0.040003327
25 30 32 40 45 50
0.056886228 0.016882901 0.036094478 0.028942116 0.009647372 0.017465070
60 65 90 125
0.008316700 0.016217565 0.016134398 0.017964072 Code
prop.table(table(database_2025$alt1.dist_trans))
1 2 4 5 6 10
0.121549894 0.074840764 0.046178344 0.093418259 0.049363057 0.134819533
12 15 20 25 30 40
0.035562633 0.047770701 0.118895966 0.019108280 0.049363057 0.081210191
50 60 80 100 120
0.054140127 0.052016985 0.012208068 0.004246285 0.005307856 Code
print("SQ Distances in 2017")[1] "SQ Distances in 2017"Code
table(database_2017$dist_3)
0.5 1 5 10 20 40 70 100 150
1440 2952 1812 1632 1380 600 372 444 1392 Code
print("SQ Distances in 2025")[1] "SQ Distances in 2025"Code
table(database_2025$alt3.dist)
0.5 1 5 10 20 40 70 100 150
204 456 312 468 216 108 36 36 48 Code
print("Check broadleaf 1")[1] "Check broadleaf 1"Code
prop.table(table(database_2025$broad_1))
0 1
0.7776008 0.2223992 Code
prop.table(table(database_2017$broad_1))
0 1
0.7762808 0.2237192 Code
print("Check broadleaf 2")[1] "Check broadleaf 2"Code
prop.table(table(database_2025$broad_2))
0 1
0.7494692 0.2505308 Code
prop.table(table(database_2017$broad_2))
0 1
0.7525782 0.2474218 Code
print("Check mix 1")[1] "Check mix 1"Code
prop.table(table(database_2025$mix_1))
0 1
0.3906582 0.6093418 Code
prop.table(table(database_2017$mix_1))
0 1
0.3897206 0.6102794 Code
print("Check mix 2")[1] "Check mix 2"Code
prop.table(table(database_2025$mix_2))
0 1
0.4437367 0.5562633 Code
prop.table(table(database_2017$mix_2))
0 1
0.4455256 0.5544744 Code
prop.table(table(database_2025$infra2_1))
0 1
0.7478769 0.2521231 Code
prop.table(table(database_2017$infra2_1))
0 1
0.7532435 0.2467565 Code
prop.table(table(database_2025$infra3_1))
0 1
0.7489384 0.2510616 Code
prop.table(table(database_2017$infra3_1))
0 1
0.7480872 0.2519128 Code
prop.table(table(database_2025$infra4_1))
0 1
0.7515924 0.2484076 Code
prop.table(table(database_2017$infra4_1))
0 1
0.7493347 0.2506653 Code
prop.table(table(database_2025$infra2_2))
0 1
0.7521231 0.2478769 Code
prop.table(table(database_2017$infra2_2))
0 1
0.7467565 0.2532435 Code
prop.table(table(database_2025$infra3_2))
0 1
0.7489384 0.2510616 Code
prop.table(table(database_2017$infra3_2))
0 1
0.7480872 0.2519128 Code
prop.table(table(database_2025$infra4_2))
0 1
0.75 0.25 Code
prop.table(table(database_2017$infra4_2))
0 1
0.75 0.25 Clogits
Code
stargazer(clg_2025_out, clg_2017_out,
type = "html",
column.labels = c("2025", "2017"),
title = "Comparison of Models (2025 vs 2017)")| Estimate | s.e. | t.rat.(0) | Rob.s.e. | Rob.t.rat.(0) | |
| type_broad | -0.223 | 0.125 | -1.782 | 0.223 | -1.003 |
| type_mix | 0.224 | 0.130 | 1.730 | 0.200 | 1.119 |
| species | 0.096 | 0.053 | 1.791 | 0.088 | 1.086 |
| height | 0.031 | 0.010 | 3.083 | 0.016 | 1.955 |
| age_two | 0.084 | 0.116 | 0.725 | 0.206 | 0.407 |
| age_multi | 0.563 | 0.108 | 5.228 | 0.196 | 2.873 |
| deadwood_med | 0.104 | 0.098 | 1.057 | 0.156 | 0.667 |
| deadwood_high | 0.319 | 0.097 | 3.290 | 0.143 | 2.226 |
| infrastr2 | -0.092 | 0.116 | -0.793 | 0.165 | -0.555 |
| infrastr3 | 0.058 | 0.090 | 0.640 | 0.164 | 0.352 |
| infrastr4 | 0.443 | 0.099 | 4.455 | 0.140 | 3.167 |
| dist | -0.014 | 0.002 | -8.133 | 0.004 | -3.163 |
| ASC_sq | 1.866 | 0.214 | 8.700 | 0.373 | 5.007 |
| Estimate | s.e. | t.rat.(0) | Rob.s.e. | Rob.t.rat.(0) | |
| type_broad | 0.017 | 0.045 | 0.367 | 0.067 | 0.247 |
| type_mix | 0.091 | 0.047 | 1.928 | 0.071 | 1.285 |
| species | 0.164 | 0.020 | 8.299 | 0.033 | 4.987 |
| height | 0.039 | 0.004 | 10.516 | 0.005 | 7.307 |
| age_two | 0.009 | 0.044 | 0.212 | 0.069 | 0.134 |
| age_multi | 0.349 | 0.040 | 8.740 | 0.065 | 5.338 |
| deadwood_med | 0.045 | 0.035 | 1.260 | 0.051 | 0.882 |
| deadwood_high | 0.211 | 0.035 | 5.952 | 0.050 | 4.260 |
| infrastr2 | -0.086 | 0.044 | -1.947 | 0.059 | -1.465 |
| infrastr3 | 0.315 | 0.033 | 9.421 | 0.061 | 5.185 |
| infrastr4 | 0.548 | 0.036 | 15.362 | 0.060 | 9.108 |
| dist | -0.014 | 0.0004 | -31.011 | 0.001 | -11.745 |
| ASC_sq | 1.132 | 0.025 | 44.840 | 0.048 | 23.661 |
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;">')Code
cat('<div style="flex: 1;">', paste(html_2025, collapse = "\n"), '</div>')| Estimate | s.e. | t.rat.(0) | Rob.s.e. | Rob.t.rat.(0) | |
| type_broad | -0.223 | 0.125 | -1.782 | 0.223 | -1.003 |
| type_mix | 0.224 | 0.130 | 1.730 | 0.200 | 1.119 |
| species | 0.096 | 0.053 | 1.791 | 0.088 | 1.086 |
| height | 0.031 | 0.010 | 3.083 | 0.016 | 1.955 |
| age_two | 0.084 | 0.116 | 0.725 | 0.206 | 0.407 |
| age_multi | 0.563 | 0.108 | 5.228 | 0.196 | 2.873 |
| deadwood_med | 0.104 | 0.098 | 1.057 | 0.156 | 0.667 |
| deadwood_high | 0.319 | 0.097 | 3.290 | 0.143 | 2.226 |
| infrastr2 | -0.092 | 0.116 | -0.793 | 0.165 | -0.555 |
| infrastr3 | 0.058 | 0.090 | 0.640 | 0.164 | 0.352 |
| infrastr4 | 0.443 | 0.099 | 4.455 | 0.140 | 3.167 |
| dist | -0.014 | 0.002 | -8.133 | 0.004 | -3.163 |
| ASC_sq | 1.866 | 0.214 | 8.700 | 0.373 | 5.007 |
Code
cat('<div style="flex: 1;">', paste(html_2017, collapse = "\n"), '</div>')| Estimate | s.e. | t.rat.(0) | Rob.s.e. | Rob.t.rat.(0) | |
| type_broad | 0.017 | 0.045 | 0.367 | 0.067 | 0.247 |
| type_mix | 0.091 | 0.047 | 1.928 | 0.071 | 1.285 |
| species | 0.164 | 0.020 | 8.299 | 0.033 | 4.987 |
| height | 0.039 | 0.004 | 10.516 | 0.005 | 7.307 |
| age_two | 0.009 | 0.044 | 0.212 | 0.069 | 0.134 |
| age_multi | 0.349 | 0.040 | 8.740 | 0.065 | 5.338 |
| deadwood_med | 0.045 | 0.035 | 1.260 | 0.051 | 0.882 |
| deadwood_high | 0.211 | 0.035 | 5.952 | 0.050 | 4.260 |
| infrastr2 | -0.086 | 0.044 | -1.947 | 0.059 | -1.465 |
| infrastr3 | 0.315 | 0.033 | 9.421 | 0.061 | 5.185 |
| infrastr4 | 0.548 | 0.036 | 15.362 | 0.060 | 9.108 |
| dist | -0.014 | 0.0004 | -31.011 | 0.001 | -11.745 |
| ASC_sq | 1.132 | 0.025 | 44.840 | 0.048 | 23.661 |
Code
cat('</div>')Code
ggplot(filter(ci_data, Group == "wtp"), aes(x = Coefficient, y = Mean, color = Model)) +
geom_point(size = 3, position = position_dodge(width = 0.6)) +
geom_errorbar(aes(ymin = Lower_CI, ymax = Upper_CI),
width = 0.2, position = position_dodge(width = 0.6)) +
geom_hline(yintercept = 0, linetype = "dashed") +
labs(
title = "WTT Estimates with 95% Confidence Intervals",
x = "Coefficient",
y = "WTT Estimate (km)",
color = "Model"
) +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
Code
ci_data_com <- ci_data %>% filter(!grepl("^ASC", Coefficient))
ggplot(filter(ci_data_com, Group == "wtp_first_vs_second"), aes(x = Coefficient, y = Mean, color = Model)) +
geom_rect(data = stripe_data[1:6, ],
aes(xmin = xmin, xmax = xmax, ymin = ymin, ymax = ymax),
fill = "gray96", inherit.aes = FALSE) +
geom_point(size = 3, position = position_dodge(width = 0.6)) +
geom_errorbar(aes(ymin = Lower_CI, ymax = Upper_CI),
width = 0.2, position = position_dodge(width = 0.6)) +
geom_hline(yintercept = 0, linetype = "dashed") +
labs(
title = "WTT Estimates with 95% Confidence Intervals First vs Second Choice",
x = "Coefficient",
y = "WTT Estimate (km)",
color = "Model"
) +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
Mixed Logits
Code
stargazer(mxl_2025_out,
type = "html",
title = "Mixed Logit Model Output (2025) WTP Space",
column.labels = "2025")| Estimate | s.e. | t.rat.(0) | Rob.s.e. | Rob.t.rat.(0) | |
| beta_type_broad | 2.501 | 4.054 | 0.617 | 4.377 | 0.571 |
| beta_type_mix | 8.313 | 4.045 | 2.055 | 3.817 | 2.178 |
| beta_age_multi | 10.627 | 3.577 | 2.971 | 4.269 | 2.489 |
| beta_age_two | 5.522 | 3.623 | 1.524 | 4.278 | 1.291 |
| beta_deadwood_high | 6.543 | 2.835 | 2.308 | 2.983 | 2.193 |
| beta_deadwood_med | 6.120 | 2.867 | 2.135 | 2.575 | 2.377 |
| beta_height | 0.878 | 0.283 | 3.097 | 0.358 | 2.453 |
| beta_infrastr2 | 2.958 | 3.620 | 0.817 | 3.940 | 0.751 |
| beta_infrastr3 | 7.405 | 3.533 | 2.096 | 4.577 | 1.618 |
| beta_infrastr4 | 15.726 | 3.683 | 4.270 | 5.065 | 3.105 |
| beta_species | 4.462 | 1.837 | 2.429 | 2.065 | 2.161 |
| beta_ASC_sq | 58.061 | 8.846 | 6.563 | 14.215 | 4.085 |
| beta_dist | 0.063 | 0.008 | 7.786 | 0.013 | 4.989 |
| sig_type_broad | 11.621 | 2.961 | 3.925 | 2.269 | 5.121 |
| sig_type_mix | -5.761 | 3.496 | -1.648 | 3.001 | -1.920 |
| sig_age_multi | 0.614 | 3.206 | 0.192 | 1.572 | 0.391 |
| sig_age_two | -4.233 | 4.850 | -0.873 | 4.786 | -0.885 |
| sig_deadwood_high | 2.159 | 4.249 | 0.508 | 2.317 | 0.932 |
| sig_deadwood_med | -5.871 | 4.337 | -1.354 | 4.718 | -1.244 |
| sig_height | -0.202 | 0.198 | -1.016 | 0.125 | -1.609 |
| sig_infrastr2 | -5.826 | 4.128 | -1.411 | 3.229 | -1.805 |
| sig_infrastr3 | 13.479 | 3.637 | 3.706 | 5.831 | 2.312 |
| sig_infrastr4 | 12.243 | 4.451 | 2.751 | 5.944 | 2.060 |
| sig_species | 6.860 | 1.661 | 4.129 | 1.953 | 3.513 |
| sig_ASC_sq | 39.144 | 5.456 | 7.175 | 9.888 | 3.959 |
| sig_dist | 0.041 | 0.006 | 6.872 | 0.009 | 4.813 |
Code
stargazer(mxl_2017_out,
type = "html",
title = "Mixed Logit Model Output (2017) WTP Space",
column.labels = "2017")| Estimate | s.e. | t.rat.(0) | Rob.s.e. | Rob.t.rat.(0) | |
| beta_type_broad | 1.517 | 0.858 | 1.769 | 1.132 | 1.341 |
| beta_type_mix | 2.314 | 0.842 | 2.748 | 0.995 | 2.327 |
| beta_age_multi | 2.098 | 0.740 | 2.835 | 0.798 | 2.630 |
| beta_age_two | -1.199 | 0.807 | -1.485 | 1.027 | -1.167 |
| beta_deadwood_high | 1.589 | 0.572 | 2.778 | 0.555 | 2.861 |
| beta_deadwood_med | 1.422 | 0.578 | 2.460 | 0.544 | 2.613 |
| beta_height | 0.948 | 0.091 | 10.452 | 0.153 | 6.205 |
| beta_infrastr2 | 0.790 | 0.827 | 0.955 | 0.914 | 0.864 |
| beta_infrastr3 | 6.982 | 0.794 | 8.795 | 1.233 | 5.664 |
| beta_infrastr4 | 9.393 | 0.855 | 10.982 | 1.507 | 6.232 |
| beta_species | 2.398 | 0.396 | 6.060 | 0.500 | 4.797 |
| beta_ASC_sq | 17.237 | 1.167 | 14.771 | 2.379 | 7.246 |
| beta_dist | 0.100 | 0.005 | 19.455 | 0.011 | 9.495 |
| sig_type_broad | 1.028 | 2.082 | 0.494 | 3.422 | 0.301 |
| sig_type_mix | -0.235 | 1.207 | -0.194 | 1.188 | -0.197 |
| sig_age_multi | 2.055 | 1.590 | 1.292 | 2.260 | 0.909 |
| sig_age_two | 0.244 | 1.379 | 0.177 | 1.647 | 0.148 |
| sig_deadwood_high | 1.126 | 1.114 | 1.011 | 1.065 | 1.057 |
| sig_deadwood_med | -0.137 | 1.553 | -0.088 | 1.950 | -0.070 |
| sig_height | 0.846 | 0.090 | 9.409 | 0.159 | 5.315 |
| sig_infrastr2 | 4.128 | 1.420 | 2.907 | 1.504 | 2.745 |
| sig_infrastr3 | 5.118 | 1.133 | 4.517 | 1.565 | 3.270 |
| sig_infrastr4 | 5.847 | 1.202 | 4.864 | 1.736 | 3.368 |
| sig_species | 1.581 | 0.959 | 1.648 | 1.797 | 0.880 |
| sig_ASC_sq | 23.723 | 1.357 | 17.476 | 3.462 | 6.853 |
| sig_dist | -0.082 | 0.005 | -17.193 | 0.010 | -8.482 |
Code
ci_data_mod <- ci_data %>% filter(Model == "2017" | Model == "2025")
ggplot(filter(ci_data_mod, Group == "mxl"), aes(x = Coefficient, y = Mean, color = Model)) +
geom_rect(data = stripe_data[1:13, ],
aes(xmin = xmin, xmax = xmax, ymin = ymin, ymax = ymax),
fill = "gray96", inherit.aes = FALSE) +
geom_point(size = 3, position = position_dodge(width = 0.6)) +
geom_errorbar(aes(ymin = Lower_CI, ymax = Upper_CI),
width = 0.2, position = position_dodge(width = 0.6)) +
geom_hline(yintercept = 0, linetype = "dashed") +
labs(
title = "WTT Estimates with 95% Confidence Intervals",
x = "Coefficient",
y = "WTT Estimate (km)",
color = "Model"
) +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
Code
ci_data_17 <- ci_data %>% filter(Model == "2017 rep") %>%
filter(grepl("^beta_", Coefficient)) # %>% filter(!grepl("^beta_ASC", Coefficient))
ggplot(filter(ci_data_17, Group == "mxl"), aes(x = Coefficient, y = Mean, color = Model), color="black") +
geom_rect(data = stripe_data[1:8, ],
aes(xmin = xmin, xmax = xmax, ymin = ymin, ymax = ymax),
fill = "gray96", inherit.aes = FALSE) +
geom_errorbar(aes(ymin = Lower_CI, ymax = Upper_CI),
width = 0.2, position = position_dodge(width = 0.6), color="black") +
geom_point(size = 2, position = position_dodge(width = 0.6), color="red") +
geom_hline(yintercept = 0, linetype = "dashed") +
labs(
title = "WTT Estimates with 95% Confidence Intervals (2017)",
x = "Coefficient",
y = "WTT Estimate (km)",
color = "Model"
) +
theme_minimal() +
coord_flip() 
Drought Questions
Code
attr(data_pre$A1, "label")[1] "A1. Have the forest areas that you regularly visit changed in the last 10 years?"Code
attr(data_pre$A1, "labels") No Yes, the condition has improved
1 2
Yes, the condition has deteriorated I can't say
3 4 Code
table(data_complete$A1)
1 2 3 4
41 21 61 34 Code
attr(data_pre$A2, "label")[1] "A2. Did you notice signs of drought damage during your last visit to the forest?"Code
attr(data_pre$A2, "labels")Yes, clearly visible Yes, a little I can't judge
1 2 3
No, hardly No, not at all
4 5 Code
table(data_complete$A2)
1 2 3 4 5
26 51 36 28 16 Code
attr(data_pre$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_pre$A3, "labels")Yes, very strongly Yes, a little I can't judge No, hardly
1 2 3 4
No, not at all
5 Code
table(data_complete$A3)
1 2 3 4 5
23 54 44 25 11 Code
attr(data_pre$A6, "label")[1] "A6. Have you changed the frequency of your forest visits due to drought damage?"Code
attr(data_pre$A6, "labels") I go to the forest much less often I go to the forest somewhat less often
1 2
My frequency of visits has not changed I go to the forest somewhat more often
3 4
I go to the forest much more often Does not apply
5 6 Code
table(data_complete$A6)
2 3 4 5 6
7 43 4 2 16 Code
attr(data_pre$A8, "label")[1] "A8. Have you changed the forest areas you usually visit because of drought damage?"Code
attr(data_pre$A8, "labels")Yes, I now visit other forest areas for recreational activities.
1
No, I have not changed the forest area I use for recreation.
2 Code
table(data_complete$A8)
1 2
11 55 Code
attr(data_pre$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_pre$A13, "labels") 0???% – 10???% 11???% – 20???% 21???% – 30???% 31???% – 40???% 41 % - 50 %
1 2 3 4 5
51 % - 60 % 61 % - 70% More than 70 % I can’t tell
6 7 8 9 Code
ggplot(data=data_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_complete_sf) %>%
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
)Code
leaflet(data_complete_sf) %>%
addProviderTiles("OpenStreetMap") %>%
addCircleMarkers(
radius = 5,
fillColor = ~pal(A3),
fillOpacity = 0.7,
color = "#444444",
weight = 0.5,
) %>%
addLegend(
"bottomright",
pal = pal,
values = ~A3,
title = "Drought Index (general)",
opacity = 0.7
)Forest Owners
Code
table(data_complete$EIG1)
1 2
8 149 Code
table(data_complete$EIG2_1)
1 2 4 10 50 12000
2 1 1 1 1 1 Code
table(data_complete$EIG3_1)
5 30 100 10000
1 2 1 1 Check trees and leaves comparison
Code
check_pair_responses <- function(data, type = c("LEAVES", "TREES"), n_pairs = 5) {
type <- match.arg(type)
# Inner function for one comparison
check_one <- function(i) {
prefix <- paste0("PAIR", type, "_", i)
show_col <- paste0("SHOWPAIR", type, "_", i)
data %>%
select(all_of(c(paste0(prefix, "_1"), paste0(prefix, "_2"), show_col))) %>%
rename(left = 1, right = 2, answer = 3) %>%
mutate(comparison = i) %>%
separate(left, into = c("left_id", NA), sep = ",", remove = FALSE) %>%
separate(right, into = c("right_id", NA), sep = ",", remove = FALSE) %>%
mutate(
answer = as.character(answer),
valid_answer = answer %in% c(left_id, right_id, "99", NA)
) %>%
select(comparison, answer, left_id, right_id, valid_answer)
}
# Apply to all comparisons
purrr::map_dfr(1:n_pairs, check_one)
}
leaf_check <- check_pair_responses(data_complete, type = "LEAVES", n_pairs = 5)
tree_check <- check_pair_responses(data_complete, type = "TREES", n_pairs = 5)
leaf_check %>% filter(!valid_answer)# A tibble: 0 × 5
# ℹ 5 variables: comparison <int>, answer <chr>, left_id <chr>, right_id <chr>,
# valid_answer <lgl>Code
tree_check %>% filter(!valid_answer)# A tibble: 0 × 5
# ℹ 5 variables: comparison <int>, answer <chr>, left_id <chr>, right_id <chr>,
# valid_answer <lgl>Code
head(leaf_check)# A tibble: 6 × 5
comparison answer left_id right_id valid_answer
<int> <chr> <chr> <chr> <lgl>
1 1 99 "67" "87" TRUE
2 1 <NA> "" "" TRUE
3 1 13 "13" "1" TRUE
4 1 81 "76" "81" TRUE
5 1 <NA> "" "" TRUE
6 1 <NA> "" "" TRUE Code
head(tree_check)# A tibble: 6 × 5
comparison answer left_id right_id valid_answer
<int> <chr> <chr> <chr> <lgl>
1 1 99 "67" "87" TRUE
2 1 <NA> "" "" TRUE
3 1 1 "13" "1" TRUE
4 1 76 "76" "81" TRUE
5 1 <NA> "" "" TRUE
6 1 <NA> "" "" TRUE