### META-REGRESSION MODELS RELATED TO THE PUBLICATION: # NOBEL, A., LIZIN, S., BROUWER, R., BRUNS, S. B., STERN, D.I., AND MALINA, R. (2020). Are biodiversity losses valued # differently when they are caused by human activities? A meta-analysis of the non-use valuation literature. Env. Res. Lett. # # First created: April 1, 2019 # Last updated: March 18, 2020 # # This R code estimates meta-regression models by means of weighted least squares with standard errors clustered by primary publication. # The meta-regression models correspond to the related article as follows: # A. Model I (X=1 if at least one human threat) -> Table 2 (first column) # B. Model II (X=1 if human threats only) -> Table 2 (second column) # C. Model I (X=1 if at least one human threat) (OUTLIERS EXCLUDED) -> Table 3 (first column) # D. Model II (X=1 if only human threats) (OUTLIERS EXCLUDED) -> Table 3 (second column) # E. Model III with control for publication bias (X=1 if at least one human threat) (OUTLIERS EXCLUDED) -> Table 3 (third column) # The data consists of 159 estimates from 62 studies of the mean willingness-to-pay for biodiversity conservation # Load required packages library(openxlsx) require(prediction) library(plotrix) require(plm) require(lmtest) # Load meta-data data1 <- read.xlsx("Meta_data.xlsx",sheet = 1) pmodel1 <- pdata.frame(data1, index = c("StudyID", "ObsID")) # Set the first index to DataID to cluster by data collection event attach(pmodel1) ## Estimate a log-linear specification (dependent variable = ln(WTP)) with Weighted Least Squares # A. Model I (X=1 if at least one human threat) baseline_model2 <- lm(Ln_WTP ~ `Species.abundance`+`Species.richness`+Marine+Wetland+Grassland_Shrubland+Birds+Mammals+Other+`WTP.Measure.is.a.gain.(Compensating.Surplus)`+`Is.the.outcome.uncertain?`+Human_cause+ `Face-to-face` + `Payment.interval.(annual)` + `Payment.vehicle.(voluntary)`+ `Unit.(per.household)`+Protest.responses.are.removed+CE+CE_Mixed.logit+ OE+`CV-Payment.Card`+`CV_Non-parametric`+Study.year+`Non-use.value.isolated.by.allocation.of.TEV/NUV.to.existence.and.bequest.value.components?`+`Non-use.value.isolated.by.excluding.users.based.on.past.or.intended.recreational.use?`+ `North-America`+ ln_GDP, weights = Rel_weights, data=pmodel1) # Calculate clustered standard errors baseline_model2_pm <- plm(as.formula(baseline_model2), data = pmodel1, index = c("StudyID", "ObsID"), model = "pooling") G <- length(unique(pmodel1$Author)) N <- nobs(baseline_model2) dfa <- (G/(G - 1)) * (N - 1)/baseline_model2_pm$df.residual pmodel1 <- pdata.frame(data1, index = c("StudyID", "ObsID")) c_vcov <- dfa * vcovHC(baseline_model2_pm, type = "HC1", cluster = "group", drop.index=F, adjust = T) coeftest(baseline_model2, vcov = c_vcov) summary(baseline_model2)$adj.r.squared summary(baseline_model2)$fstatistic pf(summary(baseline_model2)$fstatistic[1],summary(baseline_model2)$fstatistic[2],summary(baseline_model2)$fstatistic[3],lower.tail=FALSE) # Calculate F value baseline_model2_without <- lm(Ln_WTP ~ Species.abundance+`Species.richness`+Marine+Wetland+Grassland_Shrubland+Birds+Mammals+Other+`WTP.Measure.is.a.gain.(Compensating.Surplus)`+`Is.the.outcome.uncertain?`+ `Face-to-face` + `Payment.interval.(annual)` + `Payment.vehicle.(voluntary)`+ `Unit.(per.household)`+Protest.responses.are.removed+CE+CE_Mixed.logit+ OE+`CV-Payment.Card`+`CV_Non-parametric`+Study.year+`Non-use.value.isolated.by.allocation.of.TEV/NUV.to.existence.and.bequest.value.components?`+`Non-use.value.isolated.by.excluding.users.based.on.past.or.intended.recreational.use?`+ `North-America`+ ln_GDP, weights = Rel_weights, data=pmodel1) anova(baseline_model2_without,baseline_model2) # Predict WTP with this model input_values <- read.xlsx("Meta_data.xlsx",sheet =1) fitted_model <- lm(as.formula(baseline_model2, env = parent.frame()), weights=Rel_weights,data=pmodel1) data_prediction <- find_data(fitted_model) data_prediction$`Payment.interval.(annual)`=1 data_prediction$`Unit.(per.household)`=0 # Set human threat variable to 1 during prediction data_prediction$Human_cause = 1 predicted <- prediction(fitted_model,data_prediction,allow.new.levels = TRUE) # Calculate predicted sample WTP as percentage of GDP per capita (US$ PPP 2017) WTP_Percentage_of_GDP_Mean_Human_Cause = mean(exp(predicted$fitted)/data_prediction$`GDP.(GDP.per.capita.in.current.(2017).US$.PPP.in.the.data.collection.year)`) WTP_Percentage_of_GDP_95LB_Human_Cause = (WTP_Percentage_of_GDP_Mean_Human_Cause-1.96*std.error(exp(predicted$fitted)/data_prediction$`GDP.(GDP.per.capita.in.current.(2017).US$.PPP.in.the.data.collection.year)`)) WTP_Percentage_of_GDP_95UB_Human_Cause = (WTP_Percentage_of_GDP_Mean_Human_Cause+1.96*std.error(exp(predicted$fitted)/data_prediction$`GDP.(GDP.per.capita.in.current.(2017).US$.PPP.in.the.data.collection.year)`)) WTP_Percentage_of_GDP_Mean_Human_Cause WTP_Percentage_of_GDP_95LB_Human_Cause WTP_Percentage_of_GDP_95UB_Human_Cause # Set human threat variable to 0 during prediction data_prediction$Human_cause = 0 predicted <- prediction(fitted_model,data_prediction,allow.new.levels = TRUE) # Calculate predicted sample WTP as percentage of GDP per capita (US$ PPP 2017) WTP_Percentage_of_GDP_Mean_No_Human_Cause = mean(exp(predicted$fitted)/data_prediction$`GDP.(GDP.per.capita.in.current.(2017).US$.PPP.in.the.data.collection.year)`) WTP_Percentage_of_GDP_95LB_No_Human_Cause = (WTP_Percentage_of_GDP_Mean_No_Human_Cause-1.96*std.error(exp(predicted$fitted)/data_prediction$`GDP.(GDP.per.capita.in.current.(2017).US$.PPP.in.the.data.collection.year)`)) WTP_Percentage_of_GDP_95UB_No_Human_Cause = (WTP_Percentage_of_GDP_Mean_No_Human_Cause+1.96*std.error(exp(predicted$fitted)/data_prediction$`GDP.(GDP.per.capita.in.current.(2017).US$.PPP.in.the.data.collection.year)`)) WTP_Percentage_of_GDP_Mean_No_Human_Cause WTP_Percentage_of_GDP_95LB_No_Human_Cause WTP_Percentage_of_GDP_95UB_No_Human_Cause # B. Model II (X=1 if human threats only) baseline_model2 <- lm(Ln_WTP ~ `Species.abundance`+`Species.richness`+Marine+Wetland+Grassland_Shrubland+Birds+Mammals+Other+`WTP.Measure.is.a.gain.(Compensating.Surplus)`+`Is.the.outcome.uncertain?`+Human_cause_alt+ `Face-to-face` + `Payment.interval.(annual)` + `Payment.vehicle.(voluntary)`+ `Unit.(per.household)`+Protest.responses.are.removed+CE+CE_Mixed.logit+ OE+`CV-Payment.Card`+`CV_Non-parametric`+Study.year+`Non-use.value.isolated.by.allocation.of.TEV/NUV.to.existence.and.bequest.value.components?`+`Non-use.value.isolated.by.excluding.users.based.on.past.or.intended.recreational.use?`+ `North-America`+ ln_GDP, weights = Rel_weights, data=pmodel1) # Calculate clustered standard errors baseline_model2_pm <- plm(as.formula(baseline_model2), data = pmodel1, index = c("StudyID", "ObsID"), model = "pooling") G <- length(unique(pmodel1$Author)) N <- nobs(baseline_model2) dfa <- (G/(G - 1)) * (N - 1)/baseline_model2_pm$df.residual pmodel1 <- pdata.frame(data1, index = c("StudyID", "ObsID")) c_vcov <- dfa * vcovHC(baseline_model2_pm, type = "HC1", cluster = "group", drop.index=F, adjust = T) coeftest(baseline_model2, vcov = c_vcov) summary(baseline_model2)$adj.r.squared summary(baseline_model2)$fstatistic pf(summary(baseline_model2)$fstatistic[1],summary(baseline_model2)$fstatistic[2],summary(baseline_model2)$fstatistic[3],lower.tail=FALSE) # Calculate F value baseline_model2_without <- lm(Ln_WTP ~ Species.abundance+`Species.richness`+Marine+Wetland+Grassland_Shrubland+Birds+Mammals+Other+`WTP.Measure.is.a.gain.(Compensating.Surplus)`+`Is.the.outcome.uncertain?`+ `Face-to-face` + `Payment.interval.(annual)` + `Payment.vehicle.(voluntary)`+ `Unit.(per.household)`+Protest.responses.are.removed+CE+CE_Mixed.logit+ OE+`CV-Payment.Card`+`CV_Non-parametric`+Study.year+`Non-use.value.isolated.by.allocation.of.TEV/NUV.to.existence.and.bequest.value.components?`+`Non-use.value.isolated.by.excluding.users.based.on.past.or.intended.recreational.use?`+ `North-America`+ ln_GDP, weights = Rel_weights, data=pmodel1) anova(baseline_model2_without,baseline_model2) # Predict WTP with this model input_values <- read.xlsx("Meta_data.xlsx",sheet =1) fitted_model <- lm(as.formula(baseline_model2, env = parent.frame()), weights=Rel_weights,data=pmodel1) data_prediction <- find_data(fitted_model) data_prediction$`Payment.interval.(annual)`=1 data_prediction$`Unit.(per.household)`=0 # Set human threat variable to 1 during prediction data_prediction$Human_cause_alt = 1 predicted <- prediction(fitted_model,data_prediction,allow.new.levels = TRUE) # Calculate predicted sample WTP as percentage of GDP per capita (US$ PPP 2017) WTP_Percentage_of_GDP_Mean_Human_Cause = mean(exp(predicted$fitted)/data_prediction$`GDP.(GDP.per.capita.in.current.(2017).US$.PPP.in.the.data.collection.year)`) WTP_Percentage_of_GDP_95LB_Human_Cause = (WTP_Percentage_of_GDP_Mean_Human_Cause-1.96*std.error(exp(predicted$fitted)/data_prediction$`GDP.(GDP.per.capita.in.current.(2017).US$.PPP.in.the.data.collection.year)`)) WTP_Percentage_of_GDP_95UB_Human_Cause = (WTP_Percentage_of_GDP_Mean_Human_Cause+1.96*std.error(exp(predicted$fitted)/data_prediction$`GDP.(GDP.per.capita.in.current.(2017).US$.PPP.in.the.data.collection.year)`)) WTP_Percentage_of_GDP_Mean_Human_Cause WTP_Percentage_of_GDP_95LB_Human_Cause WTP_Percentage_of_GDP_95UB_Human_Cause # Set human threat variable to 0 during prediction data_prediction$Human_cause_alt = 0 predicted <- prediction(fitted_model,data_prediction,allow.new.levels = TRUE) # Calculate predicted sample WTP as percentage of GDP per capita (US$ PPP 2017) WTP_Percentage_of_GDP_Mean_No_Human_Cause = mean(exp(predicted$fitted)/data_prediction$`GDP.(GDP.per.capita.in.current.(2017).US$.PPP.in.the.data.collection.year)`) WTP_Percentage_of_GDP_95LB_No_Human_Cause = (WTP_Percentage_of_GDP_Mean_No_Human_Cause-1.96*std.error(exp(predicted$fitted)/data_prediction$`GDP.(GDP.per.capita.in.current.(2017).US$.PPP.in.the.data.collection.year)`)) WTP_Percentage_of_GDP_95UB_No_Human_Cause = (WTP_Percentage_of_GDP_Mean_No_Human_Cause+1.96*std.error(exp(predicted$fitted)/data_prediction$`GDP.(GDP.per.capita.in.current.(2017).US$.PPP.in.the.data.collection.year)`)) WTP_Percentage_of_GDP_Mean_No_Human_Cause WTP_Percentage_of_GDP_95LB_No_Human_Cause WTP_Percentage_of_GDP_95UB_No_Human_Cause # C. Model I (X=1 if at least one human threat) (OUTLIERS EXCLUDED) data_excluding_outliers <- read.xlsx("Meta_data.xlsx",sheet = 2) pmodel_excluding_outliers <- pdata.frame(data_excluding_outliers, index = c("StudyID", "ObsID")) # Set the first index to DataID to cluster by data collection event attach(pmodel_excluding_outliers) baseline_model_excluding_outliers <- lm(Ln_WTP ~ `Species.abundance`+`Species.richness`+Marine+Wetland+Grassland_Shrubland+Birds+Mammals+Other+`WTP.Measure.is.a.gain.(Compensating.Surplus)`+`Is.the.outcome.uncertain?`+Human_cause+ `Face-to-face` + `Payment.interval.(annual)` + `Payment.vehicle.(voluntary)`+ `Unit.(per.household)`+Protest.responses.are.removed+CE+CE_Mixed.logit+ OE+`CV-Payment.Card`+`CV_Non-parametric`+Study.year+`Non-use.value.isolated.by.allocation.of.TEV/NUV.to.existence.and.bequest.value.components?`+`Non-use.value.isolated.by.excluding.users.based.on.past.or.intended.recreational.use?`+ `North-America`+ ln_GDP, weights = Rel_weights, data=pmodel_excluding_outliers) # Calculate clustered standard errors baseline_model2_pm_excluding_outliers <- plm(as.formula(baseline_model_excluding_outliers), data = pmodel_excluding_outliers, index = c("StudyID", "ObsID"), model = "pooling") G <- length(unique(pmodel1$Author)) N <- nobs(baseline_model_excluding_outliers) dfa <- (G/(G - 1)) * (N - 1)/baseline_model2_pm_excluding_outliers$df.residual pmodel1 <- pdata.frame(data_excluding_outliers, index = c("StudyID", "ObsID")) c_vcov <- dfa * vcovHC(baseline_model2_pm_excluding_outliers, type = "HC1", cluster = "group", drop.index=F, adjust = T) coeftest(baseline_model_excluding_outliers, vcov = c_vcov) summary(baseline_model_excluding_outliers)$adj.r.squared summary(baseline_model_excluding_outliers)$fstatistic pf(summary(baseline_model_excluding_outliers)$fstatistic[1],summary(baseline_model_excluding_outliers)$fstatistic[2],summary(baseline_model_excluding_outliers)$fstatistic[3],lower.tail=FALSE) # Calculate F value baseline_model2_without_excluding_outliers <- lm(Ln_WTP ~ Species.abundance+`Species.richness`+Marine+Wetland+Grassland_Shrubland+Birds+Mammals+Other+`WTP.Measure.is.a.gain.(Compensating.Surplus)`+`Is.the.outcome.uncertain?`+ `Face-to-face` + `Payment.interval.(annual)` + `Payment.vehicle.(voluntary)`+ `Unit.(per.household)`+Protest.responses.are.removed+CE+CE_Mixed.logit+ OE+`CV-Payment.Card`+`CV_Non-parametric`+Study.year+`Non-use.value.isolated.by.allocation.of.TEV/NUV.to.existence.and.bequest.value.components?`+`Non-use.value.isolated.by.excluding.users.based.on.past.or.intended.recreational.use?`+ `North-America`+ ln_GDP, weights = Rel_weights, data=pmodel_excluding_outliers) anova(baseline_model2_without_excluding_outliers,baseline_model_excluding_outliers) # Predict WTP with this model input_values <- read.xlsx("Meta_data.xlsx",sheet =2) fitted_model <- lm(as.formula(baseline_model_excluding_outliers, env = parent.frame()), weights=Rel_weights,data=pmodel1) data_prediction <- find_data(fitted_model) data_prediction$`Payment.interval.(annual)`=1 data_prediction$`Unit.(per.household)`=0 # Set human threat variable to 1 during prediction data_prediction$Human_cause = 1 predicted <- prediction(fitted_model,data_prediction,allow.new.levels = TRUE) # Calculate predicted sample WTP as percentage of GDP per capita (US$ PPP 2017) WTP_Percentage_of_GDP_Mean_Human_Cause = mean(exp(predicted$fitted)/data_prediction$`GDP.(GDP.per.capita.in.current.(2017).US$.PPP.in.the.data.collection.year)`) WTP_Percentage_of_GDP_95LB_Human_Cause = (WTP_Percentage_of_GDP_Mean_Human_Cause-1.96*std.error(exp(predicted$fitted)/data_prediction$`GDP.(GDP.per.capita.in.current.(2017).US$.PPP.in.the.data.collection.year)`)) WTP_Percentage_of_GDP_95UB_Human_Cause = (WTP_Percentage_of_GDP_Mean_Human_Cause+1.96*std.error(exp(predicted$fitted)/data_prediction$`GDP.(GDP.per.capita.in.current.(2017).US$.PPP.in.the.data.collection.year)`)) WTP_Percentage_of_GDP_Mean_Human_Cause WTP_Percentage_of_GDP_95LB_Human_Cause WTP_Percentage_of_GDP_95UB_Human_Cause # Set human threat variable to 0 during prediction data_prediction$Human_cause = 0 predicted <- prediction(fitted_model,data_prediction,allow.new.levels = TRUE) # Calculate predicted sample WTP as percentage of GDP per capita (US$ PPP 2017) WTP_Percentage_of_GDP_Mean_No_Human_Cause = mean(exp(predicted$fitted)/data_prediction$`GDP.(GDP.per.capita.in.current.(2017).US$.PPP.in.the.data.collection.year)`) WTP_Percentage_of_GDP_95LB_No_Human_Cause = (WTP_Percentage_of_GDP_Mean_No_Human_Cause-1.96*std.error(exp(predicted$fitted)/data_prediction$`GDP.(GDP.per.capita.in.current.(2017).US$.PPP.in.the.data.collection.year)`)) WTP_Percentage_of_GDP_95UB_No_Human_Cause = (WTP_Percentage_of_GDP_Mean_No_Human_Cause+1.96*std.error(exp(predicted$fitted)/data_prediction$`GDP.(GDP.per.capita.in.current.(2017).US$.PPP.in.the.data.collection.year)`)) WTP_Percentage_of_GDP_Mean_No_Human_Cause WTP_Percentage_of_GDP_95LB_No_Human_Cause WTP_Percentage_of_GDP_95UB_No_Human_Cause # D. Model II (X=1 if only human threats) (OUTLIERS EXCLUDED) data_excluding_outliers <- read.xlsx("Meta_data.xlsx",sheet = 2) pmodel_excluding_outliers <- pdata.frame(data_excluding_outliers, index = c("StudyID", "ObsID")) # Set the first index to DataID to cluster by data collection event attach(pmodel_excluding_outliers) baseline_model_excluding_outliers <- lm(Ln_WTP ~ `Species.abundance`+`Species.richness`+Marine+Wetland+Grassland_Shrubland+Birds+Mammals+Other+`WTP.Measure.is.a.gain.(Compensating.Surplus)`+`Is.the.outcome.uncertain?`+Human_cause_alt+ `Face-to-face` + `Payment.interval.(annual)` + `Payment.vehicle.(voluntary)`+ `Unit.(per.household)`+Protest.responses.are.removed+CE+CE_Mixed.logit+ OE+`CV-Payment.Card`+`CV_Non-parametric`+Study.year+`Non-use.value.isolated.by.allocation.of.TEV/NUV.to.existence.and.bequest.value.components?`+`Non-use.value.isolated.by.excluding.users.based.on.past.or.intended.recreational.use?`+ `North-America`+ ln_GDP, weights = Rel_weights, data=pmodel_excluding_outliers) # Calculate clustered standard errors baseline_model2_pm_excluding_outliers <- plm(as.formula(baseline_model_excluding_outliers), data = pmodel_excluding_outliers, index = c("StudyID", "ObsID"), model = "pooling") G <- length(unique(pmodel1$Author)) N <- nobs(baseline_model_excluding_outliers) dfa <- (G/(G - 1)) * (N - 1)/baseline_model2_pm_excluding_outliers$df.residual pmodel1 <- pdata.frame(data_excluding_outliers, index = c("StudyID", "ObsID")) c_vcov <- dfa * vcovHC(baseline_model2_pm_excluding_outliers, type = "HC1", cluster = "group", drop.index=F, adjust = T) coeftest(baseline_model_excluding_outliers, vcov = c_vcov) summary(baseline_model_excluding_outliers)$adj.r.squared summary(baseline_model_excluding_outliers)$fstatistic pf(summary(baseline_model_excluding_outliers)$fstatistic[1],summary(baseline_model_excluding_outliers)$fstatistic[2],summary(baseline_model_excluding_outliers)$fstatistic[3],lower.tail=FALSE) # Calculate F value baseline_model2_without_excluding_outliers <- lm(Ln_WTP ~ Species.abundance+`Species.richness`+Marine+Wetland+Grassland_Shrubland+Birds+Mammals+Other+`WTP.Measure.is.a.gain.(Compensating.Surplus)`+`Is.the.outcome.uncertain?`+ `Face-to-face` + `Payment.interval.(annual)` + `Payment.vehicle.(voluntary)`+ `Unit.(per.household)`+Protest.responses.are.removed+CE+CE_Mixed.logit+ OE+`CV-Payment.Card`+`CV_Non-parametric`+Study.year+`Non-use.value.isolated.by.allocation.of.TEV/NUV.to.existence.and.bequest.value.components?`+`Non-use.value.isolated.by.excluding.users.based.on.past.or.intended.recreational.use?`+ `North-America`+ ln_GDP, weights = Rel_weights, data=pmodel_excluding_outliers) anova(baseline_model2_without_excluding_outliers,baseline_model_excluding_outliers) # Predict WTP with this model input_values <- read.xlsx("Meta_data.xlsx",sheet =2) fitted_model <- lm(as.formula(baseline_model_excluding_outliers, env = parent.frame()), weights=Rel_weights,data=pmodel1) data_prediction <- find_data(fitted_model) data_prediction$`Payment.interval.(annual)`=1 data_prediction$`Unit.(per.household)`=0 # Set human threat variable to 1 during prediction data_prediction$Human_cause_alt = 1 predicted <- prediction(fitted_model,data_prediction,allow.new.levels = TRUE) # Calculate predicted sample WTP as percentage of GDP per capita (US$ PPP 2017) WTP_Percentage_of_GDP_Mean_Human_Cause = mean(exp(predicted$fitted)/data_prediction$`GDP.(GDP.per.capita.in.current.(2017).US$.PPP.in.the.data.collection.year)`) WTP_Percentage_of_GDP_95LB_Human_Cause = (WTP_Percentage_of_GDP_Mean_Human_Cause-1.96*std.error(exp(predicted$fitted)/data_prediction$`GDP.(GDP.per.capita.in.current.(2017).US$.PPP.in.the.data.collection.year)`)) WTP_Percentage_of_GDP_95UB_Human_Cause = (WTP_Percentage_of_GDP_Mean_Human_Cause+1.96*std.error(exp(predicted$fitted)/data_prediction$`GDP.(GDP.per.capita.in.current.(2017).US$.PPP.in.the.data.collection.year)`)) WTP_Percentage_of_GDP_Mean_Human_Cause WTP_Percentage_of_GDP_95LB_Human_Cause WTP_Percentage_of_GDP_95UB_Human_Cause # Set human threat variable to 0 during prediction data_prediction$Human_cause_alt = 0 predicted <- prediction(fitted_model,data_prediction,allow.new.levels = TRUE) # Calculate predicted sample WTP as percentage of GDP per capita (US$ PPP 2017) WTP_Percentage_of_GDP_Mean_No_Human_Cause = mean(exp(predicted$fitted)/data_prediction$`GDP.(GDP.per.capita.in.current.(2017).US$.PPP.in.the.data.collection.year)`) WTP_Percentage_of_GDP_95LB_No_Human_Cause = (WTP_Percentage_of_GDP_Mean_No_Human_Cause-1.96*std.error(exp(predicted$fitted)/data_prediction$`GDP.(GDP.per.capita.in.current.(2017).US$.PPP.in.the.data.collection.year)`)) WTP_Percentage_of_GDP_95UB_No_Human_Cause = (WTP_Percentage_of_GDP_Mean_No_Human_Cause+1.96*std.error(exp(predicted$fitted)/data_prediction$`GDP.(GDP.per.capita.in.current.(2017).US$.PPP.in.the.data.collection.year)`)) WTP_Percentage_of_GDP_Mean_No_Human_Cause WTP_Percentage_of_GDP_95LB_No_Human_Cause WTP_Percentage_of_GDP_95UB_No_Human_Cause # E. Model III with control for publication bias (X=1 if at least one human threat) (OUTLIERS EXCLUDED) data_excluding_outliers <- read.xlsx("Meta_data.xlsx",sheet = 2) pmodel_excluding_outliers <- pdata.frame(data_excluding_outliers, index = c("StudyID", "ObsID")) # Set the first index to DataID to cluster by data collection event attach(pmodel_excluding_outliers) baseline_model_excluding_outliers <- lm(Ln_WTP ~ `Species.abundance`+`Species.richness`+Marine+Wetland+Grassland_Shrubland+Birds+Mammals+Other+`WTP.Measure.is.a.gain.(Compensating.Surplus)`+`Is.the.outcome.uncertain?`+Human_cause+ `Face-to-face` + `Payment.interval.(annual)` + `Payment.vehicle.(voluntary)`+ `Unit.(per.household)`+Protest.responses.are.removed+CE+CE_Mixed.logit+ OE+`CV-Payment.Card`+`CV_Non-parametric`+Study.year+`Non-use.value.isolated.by.allocation.of.TEV/NUV.to.existence.and.bequest.value.components?`+`Non-use.value.isolated.by.excluding.users.based.on.past.or.intended.recreational.use?`+ `North-America`+ ln_GDP+`Published.in.peer-reviewed.journal`, weights = Rel_weights, data=pmodel_excluding_outliers) # Calculate clustered standard errors baseline_model2_pm_excluding_outliers <- plm(as.formula(baseline_model_excluding_outliers), data = pmodel_excluding_outliers, index = c("StudyID", "ObsID"), model = "pooling") G <- length(unique(pmodel1$Author)) N <- nobs(baseline_model_excluding_outliers) dfa <- (G/(G - 1)) * (N - 1)/baseline_model2_pm_excluding_outliers$df.residual pmodel1 <- pdata.frame(data_excluding_outliers, index = c("StudyID", "ObsID")) c_vcov <- dfa * vcovHC(baseline_model2_pm_excluding_outliers, type = "HC1", cluster = "group", drop.index=F, adjust = T) coeftest(baseline_model_excluding_outliers, vcov = c_vcov) summary(baseline_model_excluding_outliers)$adj.r.squared summary(baseline_model_excluding_outliers)$fstatistic pf(summary(baseline_model_excluding_outliers)$fstatistic[1],summary(baseline_model_excluding_outliers)$fstatistic[2],summary(baseline_model_excluding_outliers)$fstatistic[3],lower.tail=FALSE) # Calculate F value baseline_model2_without_excluding_outliers <- lm(Ln_WTP ~ Species.abundance+`Species.richness`+Marine+Wetland+Grassland_Shrubland+Birds+Mammals+Other+`WTP.Measure.is.a.gain.(Compensating.Surplus)`+`Is.the.outcome.uncertain?`+ `Face-to-face` + `Payment.interval.(annual)` + `Payment.vehicle.(voluntary)`+ `Unit.(per.household)`+Protest.responses.are.removed+CE+CE_Mixed.logit+ OE+`CV-Payment.Card`+`CV_Non-parametric`+Study.year+`Non-use.value.isolated.by.allocation.of.TEV/NUV.to.existence.and.bequest.value.components?`+`Non-use.value.isolated.by.excluding.users.based.on.past.or.intended.recreational.use?`+ `North-America`+ ln_GDP, weights = Rel_weights, data=pmodel_excluding_outliers) anova(baseline_model2_without_excluding_outliers,baseline_model_excluding_outliers) # Predict WTP with this model input_values <- read.xlsx("Meta_data.xlsx",sheet =2) fitted_model <- lm(as.formula(baseline_model_excluding_outliers, env = parent.frame()), weights=Rel_weights,data=pmodel1) data_prediction <- find_data(fitted_model) data_prediction$`Payment.interval.(annual)`=1 data_prediction$`Unit.(per.household)`=0 # Set human threat variable to 1 during prediction data_prediction$Human_cause = 1 predicted <- prediction(fitted_model,data_prediction,allow.new.levels = TRUE) # Calculate predicted sample WTP as percentage of GDP per capita (US$ PPP 2017) WTP_Percentage_of_GDP_Mean_Human_Cause = mean(exp(predicted$fitted)/data_prediction$`GDP.(GDP.per.capita.in.current.(2017).US$.PPP.in.the.data.collection.year)`) WTP_Percentage_of_GDP_95LB_Human_Cause = (WTP_Percentage_of_GDP_Mean_Human_Cause-1.96*std.error(exp(predicted$fitted)/data_prediction$`GDP.(GDP.per.capita.in.current.(2017).US$.PPP.in.the.data.collection.year)`)) WTP_Percentage_of_GDP_95UB_Human_Cause = (WTP_Percentage_of_GDP_Mean_Human_Cause+1.96*std.error(exp(predicted$fitted)/data_prediction$`GDP.(GDP.per.capita.in.current.(2017).US$.PPP.in.the.data.collection.year)`)) WTP_Percentage_of_GDP_Mean_Human_Cause WTP_Percentage_of_GDP_95LB_Human_Cause WTP_Percentage_of_GDP_95UB_Human_Cause # Set human threat variable to 0 during prediction data_prediction$Human_cause = 0 predicted <- prediction(fitted_model,data_prediction,allow.new.levels = TRUE) # Calculate predicted sample WTP as percentage of GDP per capita (US$ PPP 2017) WTP_Percentage_of_GDP_Mean_No_Human_Cause = mean(exp(predicted$fitted)/data_prediction$`GDP.(GDP.per.capita.in.current.(2017).US$.PPP.in.the.data.collection.year)`) WTP_Percentage_of_GDP_95LB_No_Human_Cause = (WTP_Percentage_of_GDP_Mean_No_Human_Cause-1.96*std.error(exp(predicted$fitted)/data_prediction$`GDP.(GDP.per.capita.in.current.(2017).US$.PPP.in.the.data.collection.year)`)) WTP_Percentage_of_GDP_95UB_No_Human_Cause = (WTP_Percentage_of_GDP_Mean_No_Human_Cause+1.96*std.error(exp(predicted$fitted)/data_prediction$`GDP.(GDP.per.capita.in.current.(2017).US$.PPP.in.the.data.collection.year)`)) WTP_Percentage_of_GDP_Mean_No_Human_Cause WTP_Percentage_of_GDP_95LB_No_Human_Cause WTP_Percentage_of_GDP_95UB_No_Human_Cause