# Global options for R chunks
knitr::opts_chunk$set(
echo = TRUE,
fig.align = "center",
fig.height = 7,
fig.width = 10,
message = TRUE,
warning = TRUE
)
# Load necessary packages
# Ensure these are installed by running install.packages(c("readxl", "MASS", "dplyr", "ggplot2", "broom", "stringr", "forcats", "effects", "tidyr", "tibble", "stargazer")) in your R Console once.
library(readxl)
library(MASS)
library(dplyr)
##
## Attaching package: 'dplyr'
## The following object is masked from 'package:MASS':
##
## select
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(ggplot2)
library(broom)
library(stringr)
library(forcats)
library(effects)
## Loading required package: carData
## lattice theme set by effectsTheme()
## See ?effectsTheme for details.
library(tidyr)
library(tibble)
library(stargazer) # Load for professional tables
##
## Please cite as:
## Hlavac, Marek (2022). stargazer: Well-Formatted Regression and Summary Statistics Tables.
## R package version 5.2.3. https://CRAN.R-project.org/package=stargazer
# --- Load the Data ---
# IMPORTANT: Ensure the file path is correct and the Excel file is NOT open.
df <- read_excel("C:/Users/dusti/OneDrive/Desktop/Disertation Proposal/Book, Survey, Interviews Data/Survey/Stats-Quant/RQ2ProduceOffersStats.xlsx")
# --- Print Column Names (for verification) ---
cat("### Original Column Names in DataFrame:\n")
## ### Original Column Names in DataFrame:
print(colnames(df))
## [1] "Participant#" "PermanentWorkers"
## [3] "ProduceOffersRecieved" "MarketPricingImpact"
## [5] "GovPoliticiesImpactProd" "VineyardFinancialStability"
## [7] "RemittanceIncomeTotal" "RemittanceInvestmentAmtIntoVineyard"
## [9] "ProduceLostPerc" "HealthyProduceUnsoldPerc"
## [11] "MarketAccessRestricted" "FearFinancialStability"
# --- Data Cleaning and Type Coercion (Revised for Separation/Sparsity) ---
# Dependent variable: ProduceOffersRecieved
produce_offers_levels <- sort(unique(df$ProduceOffersRecieved[!is.na(df$ProduceOffersRecieved)]))
df$ProduceOffersRecieved <- ordered(df$ProduceOffersRecieved, levels = produce_offers_levels)
cat("\n### Levels of ProduceOffersRecieved (after conversion to ordered factor):\n")
##
## ### Levels of ProduceOffersRecieved (after conversion to ordered factor):
print(levels(df$ProduceOffersRecieved))
## [1] "0" "1" "2" "3" "4"
# --- RESTRUCTURE PROBLEMatic LIKERT SCALE VARIABLES TO BINARY-LIKE ---
# If value is 5 (Strongly Agree/High Impact), map to 2. Else (1,2,3,4), map to 1.
df$MarketPricingImpact_Bin <- df$MarketPricingImpact %>%
recode(`1` = 1, `2` = 1, `3` = 1, `4` = 1, `5` = 2) %>%
ordered(levels = c(1, 2))
cat("\n### Levels of MarketPricingImpact_Bin:\n")
##
## ### Levels of MarketPricingImpact_Bin:
print(levels(df$MarketPricingImpact_Bin))
## [1] "1" "2"
df$GovPoliticiesImpactProd_Bin <- df$GovPoliticiesImpactProd %>%
recode(`1` = 1, `2` = 1, `3` = 1, `4` = 1, `5` = 2) %>%
ordered(levels = c(1, 2))
cat("\n### Levels of GovPoliticiesImpactProd_Bin:\n")
##
## ### Levels of GovPoliticiesImpactProd_Bin:
print(levels(df$GovPoliticiesImpactProd_Bin))
## [1] "1" "2"
df$MarketAccessRestricted_Bin <- df$MarketAccessRestricted %>%
recode(`1` = 1, `2` = 1, `3` = 1, `4` = 1, `5` = 2) %>%
ordered(levels = c(1, 2))
cat("\n### Levels of MarketAccessRestricted_Bin:\n")
##
## ### Levels of MarketAccessRestricted_Bin:
print(levels(df$MarketAccessRestricted_Bin))
## [1] "1" "2"
df$FearFinancialStability_Bin <- df$FearFinancialStability %>%
recode(`1` = 1, `2` = 1, `3` = 1, `4` = 1, `5` = 2) %>%
ordered(levels = c(1, 2))
cat("\n### Levels of FearFinancialStability_Bin:\n")
##
## ### Levels of FearFinancialStability_Bin:
print(levels(df$FearFinancialStability_Bin))
## [1] "1" "2"
df$VineyardFinancialStability_Bin <- df$VineyardFinancialStability %>%
recode(`1` = 1, `2` = 1, `3` = 1, `4` = 1, `5` = 2) %>%
ordered(levels = c(1, 2))
cat("\n### Levels of VineyardFinancialStability_Bin:\n")
##
## ### Levels of VineyardFinancialStability_Bin:
print(levels(df$VineyardFinancialStability_Bin))
## [1] "1" "2"
df$RemittanceIncome_Binned <- df$RemittanceIncomeTotal %>%
recode(`0` = 1, `1` = 2, `2` = 3, `3` = 3) %>%
ordered(levels = c(1, 2, 3))
cat("\n### Levels of RemittanceIncome_Binned:\n")
##
## ### Levels of RemittanceIncome_Binned:
print(levels(df$RemittanceIncome_Binned))
## [1] "1" "2" "3"
# --- RESTRUCTURE SPARSE NUMERICAL VARIABLES BY BINNING ---
df$RemittanceInvestment_Binned <- df$RemittanceInvestmentAmtIntoVineyard %>%
cut(breaks = c(-Inf, 0, 1000, 3000, Inf),
labels = c("No Investment", "Low Investment", "Medium Investment", "High Investment"),
right = TRUE, include.lowest = TRUE) %>%
ordered(levels = c("No Investment", "Low Investment", "Medium Investment", "High Investment"))
cat("\n### Levels of RemittanceInvestment_Binned:\n")
##
## ### Levels of RemittanceInvestment_Binned:
print(levels(df$RemittanceInvestment_Binned))
## [1] "No Investment" "Low Investment" "Medium Investment"
## [4] "High Investment"
df$HealthyProduceUnsold_Binned <- df$HealthyProduceUnsoldPerc %>%
cut(breaks = c(-Inf, 0, 10, 25, Inf),
labels = c("None", "Low (0-10%)", "Medium (10-25%)", "High (>25%)"),
right = TRUE, include.lowest = TRUE) %>%
ordered(levels = c("None", "Low (0-10%)", "Medium (10-25%)", "High (>25%)"))
cat("\n### Levels of HealthyProduceUnsold_Binned:\n")
##
## ### Levels of HealthyProduceUnsold_Binned:
print(levels(df$HealthyProduceUnsold_Binned))
## [1] "None" "Low (0-10%)" "Medium (10-25%)" "High (>25%)"
# Numerical (ratio) variables:
df$PermanentWorkers <- as.numeric(df$PermanentWorkers)
df$ProduceLostPerc <- as.numeric(df$ProduceLostPerc)
# --- Select Relevant Columns for the Model and Handle Missing Data ---
selected_cols <- c(
"ProduceOffersRecieved",
"PermanentWorkers",
"MarketPricingImpact_Bin",
"GovPoliticiesImpactProd_Bin",
"VineyardFinancialStability_Bin",
"RemittanceIncome_Binned",
"RemittanceInvestment_Binned",
"ProduceLostPerc",
"HealthyProduceUnsold_Binned",
"MarketAccessRestricted_Bin",
"FearFinancialStability_Bin"
)
model_data <- df %>%
dplyr::select(all_of(selected_cols)) %>%
na.omit()
cat(paste("\n### Number of observations for OLR model after selection and NA removal:", nrow(model_data), "\n"))
##
## ### Number of observations for OLR model after selection and NA removal: 61
cat(paste("### Number of variables for OLR model:", ncol(model_data), "\n"))
## ### Number of variables for OLR model: 11
cat("\n### Column names in model_data (should show ALL _Bin and _Binned variables):\n")
##
## ### Column names in model_data (should show ALL _Bin and _Binned variables):
print(colnames(model_data))
## [1] "ProduceOffersRecieved" "PermanentWorkers"
## [3] "MarketPricingImpact_Bin" "GovPoliticiesImpactProd_Bin"
## [5] "VineyardFinancialStability_Bin" "RemittanceIncome_Binned"
## [7] "RemittanceInvestment_Binned" "ProduceLostPerc"
## [9] "HealthyProduceUnsold_Binned" "MarketAccessRestricted_Bin"
## [11] "FearFinancialStability_Bin"
# --- Run the Ordinal Logistic Regression Model ---
olr_model <- polr(
ProduceOffersRecieved ~ .,
data = model_data,
Hess = TRUE # Needed to calculate standard errors and p-values
)
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
# --- View Model Summary ---
summary_model <- summary(olr_model)
cat("### Model Summary:\n")
## ### Model Summary:
print(summary_model)
## Call:
## polr(formula = ProduceOffersRecieved ~ ., data = model_data,
## Hess = TRUE)
##
## Coefficients:
## Value Std. Error t value
## PermanentWorkers 0.46345 0.21631 2.14254
## MarketPricingImpact_Bin.L -0.44050 1.54479 -0.28515
## GovPoliticiesImpactProd_Bin.L 0.81322 0.72717 1.11834
## VineyardFinancialStability_Bin.L -0.13370 0.54970 -0.24322
## RemittanceIncome_Binned.L -0.06837 1.87891 -0.03639
## RemittanceIncome_Binned.Q 0.18385 1.09684 0.16762
## RemittanceInvestment_Binned.L 1.83100 1.97796 0.92570
## RemittanceInvestment_Binned.Q 0.12460 1.44879 0.08600
## RemittanceInvestment_Binned.C 0.40723 1.04364 0.39020
## ProduceLostPerc 0.02750 0.01541 1.78461
## HealthyProduceUnsold_Binned.L -1.08853 0.97265 -1.11914
## HealthyProduceUnsold_Binned.Q -0.84048 0.90722 -0.92644
## HealthyProduceUnsold_Binned.C -2.31922 0.89319 -2.59655
## MarketAccessRestricted_Bin.L 0.50671 0.50935 0.99482
## FearFinancialStability_Bin.L -0.47242 0.73669 -0.64127
##
## Intercepts:
## Value Std. Error t value
## 0|1 -0.1790 1.2669 -0.1413
## 1|2 3.4327 1.3487 2.5451
## 2|3 5.0580 1.4890 3.3970
## 3|4 6.5886 1.7781 3.7053
##
## Residual Deviance: 115.855
## AIC: 153.855
# --- Calculate p-values (polr doesn't automatically provide them) ---
coefs <- summary_model$coefficients
std_errs <- sqrt(diag(vcov(olr_model)))
coef_table <- cbind(coefs, "Std. Error" = std_errs)
p_values <- pnorm(abs(coef_table[, "t value"]), lower.tail = FALSE) * 2
coef_table <- cbind(coef_table, "p-value" = p_values)
cat("\n### Model Coefficients with p-values:\n")
##
## ### Model Coefficients with p-values:
print(coef_table)
## Value Std. Error t value Std. Error
## PermanentWorkers 0.46345427 0.21631060 2.14254070 0.21631060
## MarketPricingImpact_Bin.L -0.44049716 1.54479166 -0.28514988 1.54479166
## GovPoliticiesImpactProd_Bin.L 0.81322301 0.72716653 1.11834494 0.72716653
## VineyardFinancialStability_Bin.L -0.13369728 0.54969683 -0.24322003 0.54969683
## RemittanceIncome_Binned.L -0.06836761 1.87890737 -0.03638690 1.87890737
## RemittanceIncome_Binned.Q 0.18384944 1.09684164 0.16761712 1.09684164
## RemittanceInvestment_Binned.L 1.83100194 1.97796051 0.92570197 1.97796051
## RemittanceInvestment_Binned.Q 0.12459905 1.44878989 0.08600216 1.44878989
## RemittanceInvestment_Binned.C 0.40722879 1.04363568 0.39020206 1.04363568
## ProduceLostPerc 0.02749973 0.01540942 1.78460505 0.01540942
## HealthyProduceUnsold_Binned.L -1.08852774 0.97265078 -1.11913521 0.97265078
## HealthyProduceUnsold_Binned.Q -0.84048237 0.90721676 -0.92644052 0.90721676
## HealthyProduceUnsold_Binned.C -2.31921913 0.89319340 -2.59654755 0.89319340
## MarketAccessRestricted_Bin.L 0.50670821 0.50934696 0.99481934 0.50934696
## FearFinancialStability_Bin.L -0.47242299 0.73669415 -0.64127426 0.73669415
## 0|1 -0.17902091 1.26691713 -0.14130436 1.26691713
## 1|2 3.43270286 1.34872491 2.54514678 1.34872491
## 2|3 5.05803058 1.48896330 3.39701495 1.48896330
## 3|4 6.58856169 1.77813160 3.70532850 1.77813160
## p-value
## PermanentWorkers 0.0321499972
## MarketPricingImpact_Bin.L 0.7755293135
## GovPoliticiesImpactProd_Bin.L 0.2634197005
## VineyardFinancialStability_Bin.L 0.8078349482
## RemittanceIncome_Binned.L 0.9709738606
## RemittanceIncome_Binned.Q 0.8668845014
## RemittanceInvestment_Binned.L 0.3546008828
## RemittanceInvestment_Binned.Q 0.9314647028
## RemittanceInvestment_Binned.C 0.6963871395
## ProduceLostPerc 0.0743254018
## HealthyProduceUnsold_Binned.L 0.2630824592
## HealthyProduceUnsold_Binned.Q 0.3542170943
## HealthyProduceUnsold_Binned.C 0.0094165871
## MarketAccessRestricted_Bin.L 0.3198241378
## FearFinancialStability_Bin.L 0.5213445125
## 0|1 0.8876295069
## 1|2 0.0109231811
## 2|3 0.0006812525
## 3|4 0.0002111169
# --- Calculate Odds Ratios ---
cat("\n### Odds Ratios:\n")
##
## ### Odds Ratios:
print(exp(coef(olr_model)))
## PermanentWorkers MarketPricingImpact_Bin.L
## 1.58955526 0.64371631
## GovPoliticiesImpactProd_Bin.L VineyardFinancialStability_Bin.L
## 2.25516471 0.87485486
## RemittanceIncome_Binned.L RemittanceIncome_Binned.Q
## 0.93391709 1.20183486
## RemittanceInvestment_Binned.L RemittanceInvestment_Binned.Q
## 6.24013580 1.13269421
## RemittanceInvestment_Binned.C ProduceLostPerc
## 1.50264786 1.02788134
## HealthyProduceUnsold_Binned.L HealthyProduceUnsold_Binned.Q
## 0.33671186 0.43150233
## HealthyProduceUnsold_Binned.C MarketAccessRestricted_Bin.L
## 0.09835035 1.65981842
## FearFinancialStability_Bin.L
## 0.62348973
# Prepare data for stargazer (ensure correct variable names for labels)
stargazer(olr_model, type = "html",
title = "Ordinal Logistic Regression Results: Produce Offers Received",
dep.var.labels = "Produce Offers Received (Higher = More Offers)",
covariate.labels = c(
"Permanent Workers",
"Market Pricing Impact (Strong Agree)",
"Gov Policies Impact (Strong Agree)",
"Vineyard Financial Stability (Strong Agree)",
"Remittance Income: Low",
"Remittance Income: Med/High",
"Remittance Investment: Low",
"Remittance Investment: Medium",
"Remittance Investment: High",
"Produce Lost (%)",
"Healthy Produce Unsold: Low",
"Healthy Produce Unsold: Medium",
"Healthy Produce Unsold: High",
"Market Access Restricted (Strong Agree)",
"Fear Financial Stability (Strong Agree)"
),
ci = TRUE, ci.level = 0.90, # Display 90% CIs
notes = "Odds Ratios (exp(Coef)) are interpreted as the change in odds of receiving more produce offers for a one-unit increase in the predictor. Significance at p < 0.10.",
notes.align = "l"
)
##
## <table style="text-align:center"><caption><strong>Ordinal Logistic Regression Results: Produce Offers Received</strong></caption>
## <tr><td colspan="2" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left"></td><td><em>Dependent variable:</em></td></tr>
## <tr><td></td><td colspan="1" style="border-bottom: 1px solid black"></td></tr>
## <tr><td style="text-align:left"></td><td>Produce Offers Received (Higher = More Offers)</td></tr>
## <tr><td colspan="2" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left">Permanent Workers</td><td>0.463<sup>**</sup></td></tr>
## <tr><td style="text-align:left"></td><td>(0.108, 0.819)</td></tr>
## <tr><td style="text-align:left"></td><td></td></tr>
## <tr><td style="text-align:left">Market Pricing Impact (Strong Agree)</td><td>-0.440</td></tr>
## <tr><td style="text-align:left"></td><td>(-2.981, 2.100)</td></tr>
## <tr><td style="text-align:left"></td><td></td></tr>
## <tr><td style="text-align:left">Gov Policies Impact (Strong Agree)</td><td>0.813</td></tr>
## <tr><td style="text-align:left"></td><td>(-0.383, 2.009)</td></tr>
## <tr><td style="text-align:left"></td><td></td></tr>
## <tr><td style="text-align:left">Vineyard Financial Stability (Strong Agree)</td><td>-0.134</td></tr>
## <tr><td style="text-align:left"></td><td>(-1.038, 0.770)</td></tr>
## <tr><td style="text-align:left"></td><td></td></tr>
## <tr><td style="text-align:left">Remittance Income: Low</td><td>-0.068</td></tr>
## <tr><td style="text-align:left"></td><td>(-3.159, 3.022)</td></tr>
## <tr><td style="text-align:left"></td><td></td></tr>
## <tr><td style="text-align:left">Remittance Income: Med/High</td><td>0.184</td></tr>
## <tr><td style="text-align:left"></td><td>(-1.620, 1.988)</td></tr>
## <tr><td style="text-align:left"></td><td></td></tr>
## <tr><td style="text-align:left">Remittance Investment: Low</td><td>1.831</td></tr>
## <tr><td style="text-align:left"></td><td>(-1.422, 5.084)</td></tr>
## <tr><td style="text-align:left"></td><td></td></tr>
## <tr><td style="text-align:left">Remittance Investment: Medium</td><td>0.125</td></tr>
## <tr><td style="text-align:left"></td><td>(-2.258, 2.508)</td></tr>
## <tr><td style="text-align:left"></td><td></td></tr>
## <tr><td style="text-align:left">Remittance Investment: High</td><td>0.407</td></tr>
## <tr><td style="text-align:left"></td><td>(-1.309, 2.124)</td></tr>
## <tr><td style="text-align:left"></td><td></td></tr>
## <tr><td style="text-align:left">Produce Lost (%)</td><td>0.027<sup>*</sup></td></tr>
## <tr><td style="text-align:left"></td><td>(0.002, 0.053)</td></tr>
## <tr><td style="text-align:left"></td><td></td></tr>
## <tr><td style="text-align:left">Healthy Produce Unsold: Low</td><td>-1.089</td></tr>
## <tr><td style="text-align:left"></td><td>(-2.688, 0.511)</td></tr>
## <tr><td style="text-align:left"></td><td></td></tr>
## <tr><td style="text-align:left">Healthy Produce Unsold: Medium</td><td>-0.840</td></tr>
## <tr><td style="text-align:left"></td><td>(-2.333, 0.652)</td></tr>
## <tr><td style="text-align:left"></td><td></td></tr>
## <tr><td style="text-align:left">Healthy Produce Unsold: High</td><td>-2.319<sup>***</sup></td></tr>
## <tr><td style="text-align:left"></td><td>(-3.788, -0.850)</td></tr>
## <tr><td style="text-align:left"></td><td></td></tr>
## <tr><td style="text-align:left">Market Access Restricted (Strong Agree)</td><td>0.507</td></tr>
## <tr><td style="text-align:left"></td><td>(-0.331, 1.345)</td></tr>
## <tr><td style="text-align:left"></td><td></td></tr>
## <tr><td style="text-align:left">Fear Financial Stability (Strong Agree)</td><td>-0.472</td></tr>
## <tr><td style="text-align:left"></td><td>(-1.684, 0.739)</td></tr>
## <tr><td style="text-align:left"></td><td></td></tr>
## <tr><td colspan="2" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left">Observations</td><td>61</td></tr>
## <tr><td colspan="2" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left"><em>Note:</em></td><td style="text-align:left"><sup>*</sup>p<0.1; <sup>**</sup>p<0.05; <sup>***</sup>p<0.01</td></tr>
## <tr><td style="text-align:left"></td><td style="text-align:left">Odds Ratios (exp(Coef)) are interpreted as the change in odds of receiving more produce offers for a one-unit increase in the predictor. Significance at p < 0.10.</td></tr>
## </table>
# --- Code for Forest Plot ---
# Prepare Data for Forest Plot
odds_ratios <- exp(coef(olr_model))
conf_intervals <- confint(olr_model, level = 0.90) # Using 90% Confidence Intervals
## Waiting for profiling to be done...
odds_ratios_ci <- exp(conf_intervals)
forest_data <- data.frame(
term = names(odds_ratios),
odds_ratio = odds_ratios,
lower_ci = odds_ratios_ci[, 1],
upper_ci = odds_ratios_ci[, 2]
)
if (exists("coef_table")) {
forest_data <- forest_data %>%
filter(!str_detect(term, "\\|")) # Remove intercept terms first
p_value_lookup <- setNames(coef_table[forest_data$term, "p-value"], forest_data$term)
forest_data$`p-value` <- p_value_lookup[forest_data$term]
forest_data$sig <- ifelse(forest_data$`p-value` < 0.10, "Significant (p < 0.10)", "Not Significant") # Significance at p < 0.10
} else {
warning("coef_table not found. P-values missing for significance indication in plot.")
forest_data$sig <- "Unknown Significance"
}
forest_data <- forest_data %>%
mutate(
term = case_when(
term == "MarketPricingImpact_Bin2" ~ "Market Pricing Impact (Strong Agree)",
term == "GovPoliticiesImpactProd_Bin2" ~ "Gov Policies Impact (Strong Agree)",
term == "VineyardFinancialStability_Bin2" ~ "Vineyard Financial Stability (Strong Agree)",
term == "RemittanceIncome_Binned2" ~ "Remittance Income: Low",
term == "RemittanceIncome_Binned3" ~ "Remittance Income: Med/High",
term == "RemittanceInvestment_BinnedLow Investment" ~ "Remittance Investment: Low",
term == "RemittanceInvestment_BinnedMedium Investment" ~ "Remittance Investment: Medium",
term == "RemittanceInvestment_BinnedHigh Investment" ~ "Remittance Investment: High",
term == "HealthyProduceUnsold_BinnedLow (0-10%)" ~ "Healthy Produce Unsold: Low",
term == "HealthyProduceUnsold_BinnedMedium (10-25%)" ~ "Healthy Produce Unsold: Medium",
term == "HealthyProduceUnsold_BinnedHigh (>25%)" ~ "Healthy Produce Unsold: High",
term == "MarketAccessRestricted_Bin2" ~ "Market Access Restricted (Strong Agree)",
term == "FearFinancialStability_Bin2" ~ "Fear Financial Stability (Strong Agree)",
TRUE ~ term
),
term = fct_inorder(term)
) %>%
arrange(desc(odds_ratio))
max_finite_upper_ci <- max(forest_data$upper_ci[!is.infinite(forest_data$upper_ci)], na.rm = TRUE)
max_x_limit <- ifelse(is.finite(max_finite_upper_ci), max_finite_upper_ci * 1.5, 1000)
forest_data$lower_ci[is.infinite(forest_data$lower_ci)] <- 0.001
forest_data$upper_ci[is.infinite(forest_data$upper_ci)] <- max_x_limit
forest_plot <- ggplot(forest_data, aes(y = term, x = odds_ratio, xmin = lower_ci, xmax = upper_ci, color = sig)) +
geom_point(size = 3) +
geom_errorbarh(height = 0.2) +
geom_vline(xintercept = 1, linetype = "dashed", color = "black") +
scale_x_log10(limits = c(min(forest_data$lower_ci, na.rm=TRUE), max_x_limit)) +
labs(
title = "Odds Ratios with 90% Confidence Intervals (Final Transformed Model)",
x = "Odds Ratio (Log Scale)",
y = "Predictor Variable",
color = "Significance"
) +
theme_minimal() +
theme(
plot.title = element_text(hjust = 0.5, face = "bold"),
axis.text.y = element_text(size = 10, face = "bold"),
axis.title.x = element_text(size = 12),
axis.title.y = element_text(size = 12),
legend.position = "bottom"
) +
scale_color_manual(values = c("Significant (p < 0.10)" = "darkgreen", "Not Significant" = "darkred", "Unknown Significance" = "gray"))
print(forest_plot)
ggsave("olr_forest_plot_final_transformed_90CI.png", plot = forest_plot, width = 10, height = 8, dpi = 300)
# Example 1: MarketPricingImpact_Bin
effect_mp_impact_bin <- effects::allEffects(olr_model, latent = FALSE)
predicted_probs_mp_bin <- as.data.frame(effect_mp_impact_bin$`MarketPricingImpact_Bin`)
predicted_probs_long_mp_bin <- predicted_probs_mp_bin %>%
pivot_longer(cols = starts_with("prob."), names_to = "ProduceOffersRecieved_Category", values_to = "Predicted_Probability") %>%
mutate(ProduceOffersRecieved_Category = str_replace(ProduceOffersRecieved_Category, "prob.", "Offers = ") %>% fct_inorder())
probs_plot_mp_bin <- ggplot(predicted_probs_long_mp_bin, aes(x = MarketPricingImpact_Bin, y = Predicted_Probability, color = ProduceOffersRecieved_Category, group = ProduceOffersRecieved_Category)) +
geom_line(linewidth = 1.2) + geom_point(size = 3) +
labs(title = "Predicted Probability of Produce Offers by Market Pricing Impact (Binned)", x = "Market Pricing Impact (1=Less than Strong Agree, 2=Strong Agree)", y = "Predicted Probability", color = "Produce Offers Received") +
scale_y_continuous(labels = scales::percent, limits = c(0, 1)) + scale_x_discrete(labels = c("Less than Strong Agree", "Strong Agree")) + theme_minimal() +
theme(plot.title = element_text(hjust = 0.5, face = "bold"), axis.title.x = element_text(size = 12), axis.title.y = element_text(size = 12), legend.position = "bottom")
print(probs_plot_mp_bin)
ggsave("olr_predicted_probs_market_pricing_binned_90CI.png", plot = probs_plot_mp_bin, width = 10, height = 6, dpi = 300)
# Example 2: GovPoliticiesImpactProd_Bin
effect_gp_impact_bin <- effects::allEffects(olr_model, latent = FALSE)
predicted_probs_gp_bin <- as.data.frame(effect_gp_impact_bin$`GovPoliticiesImpactProd_Bin`)
predicted_probs_long_gp_bin <- predicted_probs_gp_bin %>%
pivot_longer(cols = starts_with("prob."), names_to = "ProduceOffersRecieved_Category", values_to = "Predicted_Probability") %>%
mutate(ProduceOffersRecieved_Category = str_replace(ProduceOffersRecieved_Category, "prob.", "Offers = ") %>% fct_inorder())
probs_plot_gp_bin <- ggplot(predicted_probs_long_gp_bin, aes(x = GovPoliticiesImpactProd_Bin, y = Predicted_Probability, color = ProduceOffersRecieved_Category, group = ProduceOffersRecieved_Category)) +
geom_line(linewidth = 1.2) + geom_point(size = 3) +
labs(title = "Predicted Probability of Produce Offers by Gov Policies Impact (Binned)", x = "Gov Policies Impact (1=Less than Strong Agree, 2=Strong Agree)", y = "Predicted Probability", color = "Produce Offers Received") +
scale_y_continuous(labels = scales::percent, limits = c(0, 1)) + scale_x_discrete(labels = c("Less than Strong Agree", "Strong Agree")) + theme_minimal() +
theme(plot.title = element_text(hjust = 0.5, face = "bold"), axis.title.x = element_text(size = 12), axis.title.y = element_text(size = 12), legend.position = "bottom")
print(probs_plot_gp_bin)
ggsave("olr_predicted_probs_gov_policies_binned_90CI.png", plot = probs_plot_gp_bin, width = 10, height = 6, dpi = 300)
# Example 3: VineyardFinancialStability_Bin
effect_vfs_bin <- effects::allEffects(olr_model, latent = FALSE)
predicted_probs_vfs_bin <- as.data.frame(effect_vfs_bin$`VineyardFinancialStability_Bin`)
predicted_probs_long_vfs_bin <- predicted_probs_vfs_bin %>%
pivot_longer(cols = starts_with("prob."), names_to = "ProduceOffersRecieved_Category", values_to = "Predicted_Probability") %>%
mutate(ProduceOffersRecieved_Category = str_replace(ProduceOffersRecieved_Category, "prob.", "Offers = ") %>% fct_inorder())
probs_plot_vfs_bin <- ggplot(predicted_probs_long_vfs_bin, aes(x = VineyardFinancialStability_Bin, y = Predicted_Probability, color = ProduceOffersRecieved_Category, group = ProduceOffersRecieved_Category)) +
geom_line(linewidth = 1.2) + geom_point(size = 3) +
labs(title = "Predicted Probability of Produce Offers by Vineyard Financial Stability (Binned)", x = "Financial Stability (1=Less than Strong Agree, 2=Strong Agree)", y = "Predicted Probability", color = "Produce Offers Received") +
scale_y_continuous(labels = scales::percent, limits = c(0, 1)) + scale_x_discrete(labels = c("Less than Strong Agree", "Strong Agree")) + theme_minimal() +
theme(plot.title = element_text(hjust = 0.5, face = "bold"), axis.title.x = element_text(size = 12), axis.title.y = element_text(size = 12), legend.position = "bottom")
print(probs_plot_vfs_bin)
ggsave("olr_predicted_probs_vineyard_stability_binned_90CI.png", plot = probs_plot_vfs_bin, width = 10, height = 6, dpi = 300)
# Example 4: RemittanceIncome_Binned
effect_rib <- effects::allEffects(olr_model, latent = FALSE)
predicted_probs_rib <- as.data.frame(effect_rib$`RemittanceIncome_Binned`)
predicted_probs_long_rib <- predicted_probs_rib %>%
pivot_longer(cols = starts_with("prob."), names_to = "ProduceOffersRecieved_Category", values_to = "Predicted_Probability") %>%
mutate(ProduceOffersRecieved_Category = str_replace(ProduceOffersRecieved_Category, "prob.", "Offers = ") %>% fct_inorder())
probs_plot_rib <- ggplot(predicted_probs_long_rib, aes(x = RemittanceIncome_Binned, y = Predicted_Probability, color = ProduceOffersRecieved_Category, group = ProduceOffersRecieved_Category)) +
geom_line(linewidth = 1.2) + geom_point(size = 3) +
labs(title = "Predicted Probability of Produce Offers by Remittance Income (Binned)", x = "Remittance Income", y = "Predicted Probability", color = "Produce Offers Received") +
scale_y_continuous(labels = scales::percent, limits = c(0, 1)) + scale_x_discrete() + theme_minimal() +
theme(plot.title = element_text(hjust = 0.5, face = "bold"), axis.title.x = element_text(size = 12), axis.title.y = element_text(size = 12), legend.position = "bottom")
print(probs_plot_rib)
ggsave("olr_predicted_probs_remittance_income_binned_90CI.png", plot = probs_plot_rib, width = 10, height = 6, dpi = 300)
# Example 5: RemittanceInvestment_Binned
effect_riib <- effects::allEffects(olr_model, latent = FALSE)
predicted_probs_riib <- as.data.frame(effect_riib$`RemittanceInvestment_Binned`)
predicted_probs_long_riib <- predicted_probs_riib %>%
pivot_longer(cols = starts_with("prob."), names_to = "ProduceOffersRecieved_Category", values_to = "Predicted_Probability") %>%
mutate(ProduceOffersRecieved_Category = str_replace(ProduceOffersRecieved_Category, "prob.", "Offers = ") %>% fct_inorder())
probs_plot_riib <- ggplot(predicted_probs_long_riib, aes(x = RemittanceInvestment_Binned, y = Predicted_Probability, color = ProduceOffersRecieved_Category, group = ProduceOffersRecieved_Category)) +
geom_line(linewidth = 1.2) + geom_point(size = 3) +
labs(title = "Predicted Probability of Produce Offers by Remittance Investment (Binned)", x = "Remittance Investment", y = "Predicted Probability", color = "Produce Offers Received") +
scale_y_continuous(labels = scales::percent, limits = c(0, 1)) + scale_x_discrete() + theme_minimal() +
theme(plot.title = element_text(hjust = 0.5, face = "bold"), axis.title.x = element_text(size = 12), axis.title.y = element_text(size = 12), legend.position = "bottom")
print(probs_plot_riib)
ggsave("olr_predicted_probs_remittance_investment_binned_90CI.png", plot = probs_plot_riib, width = 10, height = 6, dpi = 300)
# Example 6: HealthyProduceUnsold_Binned
effect_hpu_binned <- effects::allEffects(olr_model, latent = FALSE)
predicted_probs_hpu_bin <- as.data.frame(effect_hpu_binned$`HealthyProduceUnsold_Binned`)
predicted_probs_long_hpu_bin <- predicted_probs_hpu_bin %>%
pivot_longer(cols = starts_with("prob."), names_to = "ProduceOffersRecieved_Category", values_to = "Predicted_Probability") %>%
mutate(ProduceOffersRecieved_Category = str_replace(ProduceOffersRecieved_Category, "prob.", "Offers = ") %>% fct_inorder())
probs_plot_hpu_bin <- ggplot(predicted_probs_long_hpu_bin, aes(x = HealthyProduceUnsold_Binned, y = Predicted_Probability, color = ProduceOffersRecieved_Category, group = ProduceOffersRecieved_Category)) +
geom_line(linewidth = 1.2) + geom_point(size = 3) +
labs(title = "Predicted Probability of Produce Offers by Healthy Produce Unsold (Binned)", x = "Healthy Produce Unsold", y = "Predicted Probability", color = "Produce Offers Received") +
scale_y_continuous(labels = scales::percent, limits = c(0, 1)) + scale_x_discrete() + theme_minimal() +
theme(plot.title = element_text(hjust = 0.5, face = "bold"), axis.title.x = element_text(size = 12), axis.title.y = element_text(size = 12), legend.position = "bottom")
print(probs_plot_hpu_bin)
ggsave("olr_predicted_probs_healthy_unsold_binned_90CI.png", plot = probs_plot_hpu_bin, width = 10, height = 6, dpi = 300)
# Example 7: MarketAccessRestricted_Bin
effect_ma_restricted_bin <- effects::allEffects(olr_model, latent = FALSE)
predicted_probs_ma_bin <- as.data.frame(effect_ma_restricted_bin$`MarketAccessRestricted_Bin`)
predicted_probs_long_ma_bin <- predicted_probs_ma_bin %>%
pivot_longer(cols = starts_with("prob."), names_to = "ProduceOffersRecieved_Category", values_to = "Predicted_Probability") %>%
mutate(ProduceOffersRecieved_Category = str_replace(ProduceOffersRecieved_Category, "prob.", "Offers = ") %>% fct_inorder())
probs_plot_ma_bin <- ggplot(predicted_probs_long_ma_bin, aes(x = MarketAccessRestricted_Bin, y = Predicted_Probability, color = ProduceOffersRecieved_Category, group = ProduceOffersRecieved_Category)) +
geom_line(linewidth = 1.2) + geom_point(size = 3) +
labs(title = "Predicted Probability of Produce Offers by Market Access Restriction (Binned)", x = "Market Access Restriction (1=Less than Strong Agree, 2=Strong Agree)", y = "Predicted Probability", color = "Produce Offers Received") +
scale_y_continuous(labels = scales::percent, limits = c(0, 1)) + scale_x_discrete(labels = c("Less than Strong Agree", "Strong Agree")) + theme_minimal() +
theme(plot.title = element_text(hjust = 0.5, face = "bold"), axis.title.x = element_text(size = 12), axis.title.y = element_text(size = 12), legend.position = "bottom")
print(probs_plot_ma_bin)
ggsave("olr_predicted_probs_market_access_binned_90CI.png", plot = probs_plot_ma_bin, width = 10, height = 6, dpi = 300)
# Example 8: FearFinancialStability_Bin
effect_ffs_bin <- effects::allEffects(olr_model, latent = FALSE)
predicted_probs_ffs_bin <- as.data.frame(effect_ffs_bin$`FearFinancialStability_Bin`)
predicted_probs_long_ffs_bin <- predicted_probs_ffs_bin %>%
pivot_longer(cols = starts_with("prob."), names_to = "ProduceOffersRecieved_Category", values_to = "Predicted_Probability") %>%
mutate(ProduceOffersRecieved_Category = str_replace(ProduceOffersRecieved_Category, "prob.", "Offers = ") %>% fct_inorder())
probs_plot_ffs_bin <- ggplot(predicted_probs_long_ffs_bin, aes(x = FearFinancialStability_Bin, y = Predicted_Probability, color = ProduceOffersRecieved_Category, group = ProduceOffersRecieved_Category)) +
geom_line(linewidth = 1.2) + geom_point(size = 3) +
labs(title = "Predicted Probability of Produce Offers by Fear Financial Stability (Binned)", x = "Fear Financial Stability (1=Less than Strong Agree, 2=Strong Agree)", y = "Predicted Probability", color = "Produce Offers Received") +
scale_y_continuous(labels = scales::percent, limits = c(0, 1)) + scale_x_discrete(labels = c("Less than Strong Agree", "Strong Agree")) + theme_minimal() +
theme(plot.title = element_text(hjust = 0.5, face = "bold"), axis.title.x = element_text(size = 12), axis.title.y = element_text(size = 12), legend.position = "bottom")
print(probs_plot_ffs_bin)
ggsave("olr_predicted_probs_fear_financial_stability_binned_90CI.png", plot = probs_plot_ffs_bin, width = 10, height = 6, dpi = 300)
# --- Final Check: Session Information (Useful for debugging environment issues) ---
sessionInfo()
## R version 4.5.1 (2025-06-13 ucrt)
## Platform: x86_64-w64-mingw32/x64
## Running under: Windows 11 x64 (build 26100)
##
## Matrix products: default
## LAPACK version 3.12.1
##
## locale:
## [1] LC_COLLATE=English_United States.utf8
## [2] LC_CTYPE=English_United States.utf8
## [3] LC_MONETARY=English_United States.utf8
## [4] LC_NUMERIC=C
## [5] LC_TIME=English_United States.utf8
##
## time zone: America/Chicago
## tzcode source: internal
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] stargazer_5.2.3 tibble_3.3.0 tidyr_1.3.1 effects_4.2-2
## [5] carData_3.0-5 forcats_1.0.0 stringr_1.5.1 broom_1.0.8
## [9] ggplot2_3.5.2 dplyr_1.1.4 MASS_7.3-65 readxl_1.4.5
##
## loaded via a namespace (and not attached):
## [1] sass_0.4.10 generics_0.1.4 stringi_1.8.7 lattice_0.22-7
## [5] lme4_1.1-37 digest_0.6.37 magrittr_2.0.3 evaluate_1.0.4
## [9] grid_4.5.1 RColorBrewer_1.1-3 fastmap_1.2.0 cellranger_1.1.0
## [13] jsonlite_2.0.0 Matrix_1.7-3 nnet_7.3-20 backports_1.5.0
## [17] DBI_1.2.3 survival_3.8-3 purrr_1.1.0 scales_1.4.0
## [21] textshaping_1.0.1 jquerylib_0.1.4 Rdpack_2.6.4 reformulas_0.4.1
## [25] cli_3.6.5 mitools_2.4 rlang_1.1.6 rbibutils_2.3
## [29] splines_4.5.1 withr_3.0.2 cachem_1.1.0 yaml_2.3.10
## [33] tools_4.5.1 nloptr_2.2.1 minqa_1.2.8 colorspace_2.1-1
## [37] boot_1.3-31 vctrs_0.6.5 R6_2.6.1 lifecycle_1.0.4
## [41] ragg_1.4.0 insight_1.3.1 pkgconfig_2.0.3 pillar_1.11.0
## [45] bslib_0.9.0 gtable_0.3.6 glue_1.8.0 Rcpp_1.1.0
## [49] systemfonts_1.2.3 xfun_0.52 tidyselect_1.2.1 rstudioapi_0.17.1
## [53] knitr_1.50 farver_2.1.2 survey_4.4-2 htmltools_0.5.8.1
## [57] nlme_3.1-168 labeling_0.4.3 rmarkdown_2.29 compiler_4.5.1