rscript.
R
lenovo pc
2026-04-14
# ============================================================
# AGRO136: Introduction to Statistics - Practical Assignment
# Totonga Community Livelihoods Survey - CORRECTED SCRIPT
# ============================================================
# Load required libraries
library(readxl)
library(ggplot2)
library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(car)
## Loading required package: carData
##
## Attaching package: 'car'
## The following object is masked from 'package:dplyr':
##
## recode
library(corrplot)
## corrplot 0.95 loaded
# Set seed for reproducibility
[Link](2026)
# Create directories
[Link]("plots", showWarnings = FALSE)
[Link]("outputs", showWarnings = FALSE)
# ============================================================
# 1. DATA IMPORT
# ============================================================
# Import the Excel data
raw_data <- read_excel("C:/Users/lenovo pc/Desktop/Totonga_Community_Li
velihoods_Survey_-_all_versions_-_labels_-_2026-[Link]")
cat("============================================================\n")
## ============================================================
cat("AGRO136 - Totonga Community Livelihoods Survey Analysis\n")
## AGRO136 - Totonga Community Livelihoods Survey Analysis
cat("============================================================\n\n")
## ============================================================
cat("Original data:", nrow(raw_data), "rows,", ncol(raw_data), "columns
\n")
## Original data: 39 rows, 66 columns
# Check consent column values
consent_col <- raw_data$`Do you agree to participate in this survey vol
untarily`
cat("\nUnique consent values:", paste(unique(consent_col), collapse=",
"), "\n")
##
## Unique consent values: NA
# ============================================================
# 2. DATA CLEANING - FIXED
# ============================================================
# Filter valid responses (consent = "yes" OR "Yes" OR "Yes.")
clean_data <- raw_data[tolower(consent_col) == "yes", ]
# Also include rows where consent might be in a different column format
if(nrow(clean_data) == 0) {
# Try alternative consent column names
consent_cols <- grep("agree", names(raw_data), value = TRUE, ignore.c
ase = TRUE)
for(col in consent_cols) {
temp_data <- raw_data[tolower(raw_data[[col]]) == "yes", ]
if(nrow(temp_data) > nrow(clean_data)) {
clean_data <- temp_data
}
}
}
cat("After consent filter:", nrow(clean_data), "rows\n")
## After consent filter: 39 rows
# Remove rows with missing critical data (age column)
age_col <- grep("age", names(raw_data), value = TRUE, [Link] = TRU
E)[1]
if(length(age_col) > 0) {
clean_data <- clean_data[, ]
}
cat("After removing missing age:", nrow(clean_data), "rows\n")
## After removing missing age: 0 rows
# If still no data, use all data (the file might not have consent issue
s)
if(nrow(clean_data) == 0) {
cat("WARNING: No valid consent found. Using all data for analysis.\n")
clean_data <- raw_data
}
## WARNING: No valid consent found. Using all data for analysis.
# Create derived variables with clean names
clean_data$hh_size <- [Link](clean_data$`How many people live in th
e household?`)
clean_data$age <- [Link](clean_data$`Age of respondent`)
clean_data$gender <- [Link](clean_data$`Gender of respondent`)
clean_data$education <- [Link](clean_data$`Highest education level`)
clean_data$marital_status <- [Link](clean_data$`Marital status`)
clean_data$land_access <- clean_data$`Does the household have access to
land?`
clean_data$land_ha <- [Link](clean_data$`How many hectares of land
are cultivated?`)
clean_data$irrigation <- [Link](clean_data$`Do you use irrigation?`)
clean_data$water_source <- clean_data$`Main source of farming water`
clean_data$fertilizer_use <- [Link](clean_data$`Did you use fertiliz
er last season?`)
clean_data$seed_type <- clean_data$`Type of seed mostly used`
clean_data$extension <- [Link](clean_data$`Have you received extensi
on services?`)
clean_data$crop_income <- [Link](clean_data$`Estimated income from
crops last season (USD)`)
clean_data$livestock_income <- [Link](clean_data$`Estimated income
from livestock last year (USD)`)
clean_data$total_income <- clean_data$crop_income + clean_data$livestoc
k_income
clean_data$income_stable <- clean_data$`Is your income stable throughou
t the year?`
clean_data$food_shortage <- [Link](clean_data$`Did your household ex
perience food shortage in last 12 months?`)
clean_data$drought <- [Link](clean_data$`Did drought affect your far
ming in last 3 years?`)
clean_data$floods <- [Link](clean_data$`Did floods affect your farmi
ng in last 3 years?`)
clean_data$mobile_access <- [Link](clean_data$`Do you have access to
a mobile phone?`)
clean_data$market_access <- [Link](clean_data$`Do you have access to
agricultural markets?`)
clean_data$market_distance <- [Link](clean_data$`Distance to neares
t market (km)`)
clean_data$credit_access <- [Link](clean_data$`Do you have access to
agricultural credit?`)
clean_data$training_desire <- [Link](clean_data$`Would you like more
agricultural training?`)
clean_data$farming_hours <- [Link](clean_data$`Hours per day spent
farming`)
clean_data$hh_work_hours <- [Link](clean_data$`Hours per day spent
on household work`)
clean_data$challenge <- clean_data$`What is the biggest challenge facin
g your farming activities?`
clean_data$suggestion <- clean_data$`What suggestions do you have to im
prove livelihoods in the community?`
# Fix for pest columns - check if they exist
pest_cols <- grep("Which pests affected", names(clean_data), value = TR
UE)
if(length(pest_cols) > 0) {
clean_data$pest_presence <- "No"
for(pest in pest_cols) {
clean_data$pest_presence[clean_data[[pest]] == 1] <- "Yes"
}
clean_data$pest_presence <- [Link](clean_data$pest_presence)
} else {
clean_data$pest_presence <- [Link]("Unknown")
}
# Create binary variables for regression
clean_data$food_shortage_binary <- ifelse(clean_data$food_shortage == "
Yes", 1, 0)
clean_data$irrigation_binary <- ifelse(clean_data$irrigation == "Yes",
1, 0)
# Remove rows with NA income for plotting
clean_data <- clean_data[, ]
# Remove income outliers (99th percentile) - only if there are enough o
bservations
if(nrow(clean_data) > 10) {
income_99 <- quantile(clean_data$total_income, 0.99, [Link] = TRUE)
clean_data <- clean_data[clean_data$total_income < income_99, ]
}
cat("\nFinal sample size:", nrow(clean_data), "households\n\n")
##
## Final sample size: 36 households
# Check if we have data
if(nrow(clean_data) == 0) {
cat("ERROR: No data available after cleaning. Please check your Excel
file.\n")
cat("The file has", nrow(raw_data), "rows. Here are the first few row
s:\n")
print(head(raw_data[, 1:5]))
stop("Analysis stopped: No valid data.")
}
# Save cleaned data
[Link](clean_data, "cleaned_totonga_data.csv", [Link] = FALSE)
# ============================================================
# 3. CODEBOOK
# ============================================================
codebook <- [Link](
Variable = c("hh_size", "age", "gender", "education", "land_ha", "irr
igation",
"crop_income", "livestock_income", "total_income", "food
_shortage",
"drought", "pest_presence", "farming_hours", "hh_work_ho
urs"),
Description = c("Household size (number of people)", "Age of responde
nt (years)",
"Gender of respondent", "Highest education level",
"Cultivated land (hectares)", "Uses irrigation",
"Income from crops (USD)", "Income from livestock (US
D)",
"Total household income (USD)", "Experienced food sho
rtage",
"Drought affected farming", "Pest problems present",
"Hours per day farming", "Hours per day on household
work"),
Type = c("Numeric", "Numeric", "Categorical", "Categorical", "Numeric
",
"Categorical", "Numeric", "Numeric", "Numeric", "Categorical
",
"Categorical", "Categorical", "Numeric", "Numeric")
)
[Link](codebook, "[Link]", [Link] = FALSE)
cat("Codebook saved with", nrow(codebook), "variables\n\n")
## Codebook saved with 14 variables
# ============================================================
# 4. DESCRIPTIVE STATISTICS
# ============================================================
cat("============================================================\n")
## ============================================================
cat("DESCRIPTIVE STATISTICS\n")
## DESCRIPTIVE STATISTICS
cat("============================================================\n\n")
## ============================================================
# Numeric summaries - only if data exists
if(nrow(clean_data) > 0) {
cat("Household Size:\n")
print(summary(clean_data$hh_size))
cat("\nCultivated Land (hectares):\n")
print(summary(clean_data$land_ha))
cat("\nTotal Income (USD):\n")
print(summary(clean_data$total_income))
cat("\nFarming Hours per Day:\n")
print(summary(clean_data$farming_hours))
cat("\nHousehold Work Hours per Day:\n")
print(summary(clean_data$hh_work_hours))
# Categorical frequencies
cat("\n=== GENDER ===\n")
gender_table <- table(clean_data$gender)
print(gender_table)
if(length(gender_table) > 0) {
cat("Percentages:\n")
print(round([Link](gender_table) * 100, 1))
}
cat("\n=== EDUCATION ===\n")
edu_table <- table(clean_data$education)
print(edu_table)
if(length(edu_table) > 0) {
cat("Percentages:\n")
print(round([Link](edu_table) * 100, 1))
}
cat("\n=== IRRIGATION ACCESS ===\n")
irr_table <- table(clean_data$irrigation)
print(irr_table)
if(length(irr_table) > 0) {
cat("Percentages:\n")
print(round([Link](irr_table) * 100, 1))
}
cat("\n=== FOOD SHORTAGE ===\n")
food_table <- table(clean_data$food_shortage)
print(food_table)
if(length(food_table) > 0) {
cat("Percentages:\n")
print(round([Link](food_table) * 100, 1))
}
cat("\n=== EXTENSION SERVICES ===\n")
ext_table <- table(clean_data$extension)
print(ext_table)
if(length(ext_table) > 0) {
cat("Percentages:\n")
print(round([Link](ext_table) * 100, 1))
}
# Summary table for report
summary_stats <- [Link](
Variable = c("Household Size", "Land (ha)", "Total Income (USD)",
"Farming Hours", "Household Work Hours"),
Mean = round(c(mean(clean_data$hh_size, [Link] = TRUE),
mean(clean_data$land_ha, [Link] = TRUE),
mean(clean_data$total_income, [Link] = TRUE),
mean(clean_data$farming_hours, [Link] = TRUE),
mean(clean_data$hh_work_hours, [Link] = TRUE)), 2),
Median = round(c(median(clean_data$hh_size, [Link] = TRUE),
median(clean_data$land_ha, [Link] = TRUE),
median(clean_data$total_income, [Link] = TRUE),
median(clean_data$farming_hours, [Link] = TRUE),
median(clean_data$hh_work_hours, [Link] = TRUE)),
2),
SD = round(c(sd(clean_data$hh_size, [Link] = TRUE),
sd(clean_data$land_ha, [Link] = TRUE),
sd(clean_data$total_income, [Link] = TRUE),
sd(clean_data$farming_hours, [Link] = TRUE),
sd(clean_data$hh_work_hours, [Link] = TRUE)), 2)
)
print(summary_stats)
[Link](summary_stats, "outputs/descriptive_stats.csv", [Link] =
FALSE)
}
## Household Size:
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 3.000 5.000 6.000 9.028 9.250 56.000
##
## Cultivated Land (hectares):
## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 0.5 2.0 4.0 186.6 10.0 2500.0 3
##
## Total Income (USD):
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 130 1075 3500 4155 6250 11000
##
## Farming Hours per Day:
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1.000 4.000 5.000 5.222 6.000 8.000
##
## Household Work Hours per Day:
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1.000 2.000 3.000 3.083 3.000 8.000
##
## === GENDER ===
##
## Female Male
## 18 18
## Percentages:
##
## Female Male
## 50 50
##
## === EDUCATION ===
##
## No formal education Primary Secondary
Tertiary
## 2 1 14
19
## Percentages:
##
## No formal education Primary Secondary
Tertiary
## 5.6 2.8 38.9
52.8
##
## === IRRIGATION ACCESS ===
##
## No Yes
## 6 30
## Percentages:
##
## No Yes
## 16.7 83.3
##
## === FOOD SHORTAGE ===
##
## No Yes
## 25 10
## Percentages:
##
## No Yes
## 71.4 28.6
##
## === EXTENSION SERVICES ===
##
## No Yes
## 13 23
## Percentages:
##
## No Yes
## 36.1 63.9
## Variable Mean Median SD
## 1 Household Size 9.03 6 9.04
## 2 Land (ha) 186.62 4 595.61
## 3 Total Income (USD) 4155.06 3500 3268.77
## 4 Farming Hours 5.22 5 1.66
## 5 Household Work Hours 3.08 3 1.63
# ============================================================
# 5. VISUALISATIONS (Plots appear in RStudio Plots pane)
# ============================================================
cat("\n============================================================\n")
##
## ============================================================
cat("GENERATING PLOTS\n")
## GENERATING PLOTS
cat("============================================================\n")
## ============================================================
cat("Plots will appear in the Plots pane.\n\n")
## Plots will appear in the Plots pane.
# PLOT 1: Income Distribution
cat("Plot 1: Income Distribution\n")
## Plot 1: Income Distribution
if(nrow(clean_data) > 0 && sum() > 0) {
p1 <- ggplot(clean_data[ & clean_data
$total_income > 0, ],
aes(x = total_income)) +
geom_histogram(bins = 15, fill = "#3498DB", color = "black", alpha
= 0.7) +
labs(title = "Distribution of Total Household Income",
subtitle = paste("Based on", nrow(clean_data), "households"),
x = "Total Income (USD)",
y = "Number of Households") +
theme_minimal() +
theme([Link] = element_text(hjust = 0.5, face = "bold", size =
14),
[Link] = element_text(hjust = 0.5, size = 10))
print(p1)
ggsave("plots/income_distribution.png", p1, width = 8, height = 6)
cat("✓ Saved: plots/income_distribution.png\n")
} else {
cat("✗ Not enough income data for plot\n")
}
## ✓ Saved: plots/income_distribution.png
readline(prompt = "Press Enter for next plot...")
## Press Enter for next plot...
## [1] ""
# PLOT 2: Food Shortage by Irrigation
cat("\nPlot 2: Food Shortage by Irrigation Access\n")
##
## Plot 2: Food Shortage by Irrigation Access
if(nrow(clean_data) > 0 && sum() > 0 &&
sum() > 0) {
plot_data <- clean_data[ & , ]
p2 <- ggplot(plot_data, aes(x = irrigation, fill = food_shortage)) +
geom_bar(position = "fill") +
scale_fill_manual(values = c("Yes" = "#E74C3C", "No" = "#2ECC71"),
[Link] = "grey") +
labs(title = "Food Shortage by Irrigation Access",
subtitle = "Proportion of households reporting food shortage",
x = "Uses Irrigation",
y = "Proportion of Households",
fill = "Food Shortage") +
theme_minimal() +
theme([Link] = element_text(hjust = 0.5, face = "bold", size =
14))
print(p2)
ggsave("plots/food_shortage_by_irrigation.png", p2, width = 8, height
= 6)
cat("✓ Saved: plots/food_shortage_by_irrigation.png\n")
} else {
cat("✗ Not enough data for irrigation/food shortage plot\n")
}
## ✓ Saved: plots/food_shortage_by_irrigation.png
readline(prompt = "Press Enter for next plot...")
## Press Enter for next plot...
## [1] ""
# PLOT 3: Income by Education Level
cat("\nPlot 3: Income by Education Level\n")
##
## Plot 3: Income by Education Level
if(nrow(clean_data) > 0 && sum() > 0 &&
sum() > 0) {
plot_data <- clean_data[ &  &
clean_data$total_income > 0, ]
if(nrow(plot_data) > 0) {
p3 <- ggplot(plot_data, aes(x = education, y = total_income, fill =
education)) +
geom_boxplot(alpha = 0.7) +
labs(title = "Total Income by Education Level",
subtitle = "Boxplots show median, IQR, and outliers",
x = "Education Level",
y = "Total Income (USD)") +
theme_minimal() +
theme([Link] = element_text(hjust = 0.5, face = "bold", size
= 14),
[Link].x = element_text(angle = 45, hjust = 1),
[Link] = "none")
print(p3)
ggsave("plots/income_by_education.png", p3, width = 10, height = 6)
cat("✓ Saved: plots/income_by_education.png\n")
} else {
cat("✗ Not enough data for income by education plot\n")
}
}
## ✓ Saved: plots/income_by_education.png
readline(prompt = "Press Enter for next plot...")
## Press Enter for next plot...
## [1] ""
# PLOT 4: Time Allocation
cat("\nPlot 4: Time Allocation\n")
##
## Plot 4: Time Allocation
if(nrow(clean_data) > 0 && sum() > 0 &&
sum() > 0) {
plot_data <- clean_data[ & , ]
if(nrow(plot_data) > 1) {
cor_value <- cor(plot_data$farming_hours, plot_data$hh_work_hours,
use = "[Link]")
p4 <- ggplot(plot_data, aes(x = farming_hours, y = hh_work_hours))
+
geom_point(alpha = 0.6, color = "#2C3E50", size = 3) +
geom_smooth(method = "lm", se = TRUE, color = "#E74C3C", fill = "
#F1948A") +
labs(title = "Time Allocation: Farming vs Household Work",
subtitle = paste("Correlation:", round(cor_value, 2)),
x = "Hours per Day Spent Farming",
y = "Hours per Day on Household Work") +
theme_minimal() +
theme([Link] = element_text(hjust = 0.5, face = "bold", size
= 14))
print(p4)
ggsave("plots/time_allocation.png", p4, width = 8, height = 6)
cat("✓ Saved: plots/time_allocation.png\n")
} else {
cat("✗ Not enough data for time allocation plot\n")
}
}
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## ✓ Saved: plots/time_allocation.png
readline(prompt = "Press Enter for next plot...")
## Press Enter for next plot...
## [1] ""
# PLOT 5: Income by Gender
cat("\nPlot 5: Income by Gender\n")
##
## Plot 5: Income by Gender
if(nrow(clean_data) > 0 && sum() > 0 &&
sum() > 0) {
plot_data <- clean_data[ &  &
clean_data$total_income > 0, ]
if(nrow(plot_data) > 0) {
p5 <- ggplot(plot_data, aes(x = gender, y = total_income, fill = ge
nder)) +
geom_boxplot(alpha = 0.7) +
labs(title = "Income Distribution by Gender",
subtitle = paste("Female n=", sum(plot_data$gender == "Femal
e", [Link] = TRUE),
"| Male n=", sum(plot_data$gender == "Male",
[Link] = TRUE)),
x = "Gender",
y = "Total Income (USD)") +
theme_minimal() +
theme([Link] = element_text(hjust = 0.5, face = "bold", size
= 14),
[Link] = "none")
print(p5)
ggsave("plots/income_by_gender.png", p5, width = 8, height = 6)
cat("✓ Saved: plots/income_by_gender.png\n")
} else {
cat("✗ Not enough data for income by gender plot\n")
}
}
## ✓ Saved: plots/income_by_gender.png
cat("\n✓ Plots saved to 'plots/' folder\n\n")
##
## ✓ Plots saved to 'plots/' folder
# ============================================================
# 6. STATISTICAL COMPARISONS (Only if enough data)
# ============================================================
if(nrow(clean_data) >= 10) {
cat("============================================================\n")
cat("STATISTICAL COMPARISONS\n")
cat("============================================================\n\n
")
# T-test: Income by food shortage
if(length(unique(clean_data$food_shortage)) >= 2) {
t_test_income <- [Link](total_income ~ food_shortage, data = clean_
data)
cat("T-TEST: Income by Food Shortage Status\n")
cat("t =", round(t_test_income$statistic, 3),
", df =", round(t_test_income$parameter, 1),
", p-value =", format(t_test_income$[Link], scientific = TRUE,
digits = 3), "\n\n")
}
# Chi-square: Irrigation and food shortage
if(length(unique(clean_data$irrigation)) >= 2 && length(unique(clean_
data$food_shortage)) >= 2) {
chi_table <- table(clean_data$irrigation, clean_data$food_shortage)
if(sum(chi_table) > 0) {
chi_square <- [Link](chi_table)
cat("CHI-SQUARE: Irrigation vs Food Shortage\n")
cat("X-squared =", round(chi_square$statistic, 3),
", df =", chi_square$parameter,
", p-value =", format(chi_square$[Link], scientific = TRUE,
digits = 3), "\n")
print(chi_table)
cat("\n")
}
}
# ============================================================
# 7. LINEAR REGRESSION MODEL
# ============================================================
cat("============================================================\n")
cat("LINEAR REGRESSION MODEL\n")
cat("============================================================\n")
cat("Dependent Variable: Total Household Income (USD)\n\n")
# Prepare data for regression
reg_data <- clean_data[ &
 &
 &
, ]
if(nrow(reg_data) >= 10) {
lm_model <- lm(total_income ~ land_ha + irrigation + hh_size, data
= reg_data)
summary_lm <- summary(lm_model)
print(summary_lm)
lm_coef <- [Link](
Variable = rownames(summary_lm$coefficients),
Coefficient = round(summary_lm$coefficients[, 1], 2),
Std_Error = round(summary_lm$coefficients[, 2], 2),
p_value = format(summary_lm$coefficients[, 4], scientific = TRUE,
digits = 3)
)
[Link](lm_coef, "outputs/linear_regression.csv", [Link] = FAL
SE)
cat("\nR-squared =", round(summary_lm$[Link], 3), "\n")
} else {
cat("Not enough data for linear regression\n")
}
# ============================================================
# 8. LOGISTIC REGRESSION MODEL
# ============================================================
cat("\n============================================================\n
")
cat("LOGISTIC REGRESSION MODEL\n")
cat("============================================================\n")
cat("Dependent Variable: Food Shortage (1=Yes, 0=No)\n\n")
# Prepare data for logistic regression
log_data <- clean_data[ &
 &
 &
, ]
if(nrow(log_data) >= 10 && length(unique(log_data$food_shortage_binar
y)) >= 2) {
log_model <- glm(food_shortage_binary ~ total_income + irrigation +
drought,
data = log_data, family = binomial)
summary_log <- summary(log_model)
print(summary_log)
odds_ratios <- exp(coef(log_model))
odds_table <- [Link](
Variable = names(odds_ratios),
Odds_Ratio = round(odds_ratios, 3),
P_value = format(summary_log$coefficients[, 4], scientific = TRUE,
digits = 3)
)
[Link](odds_table, "outputs/logistic_regression.csv", [Link]
= FALSE)
} else {
cat("Not enough data for logistic regression\n")
}
}
## ============================================================
## STATISTICAL COMPARISONS
## ============================================================
##
## T-TEST: Income by Food Shortage Status
## t = 1.427 , df = 23 , p-value = 1.67e-01
## Warning in [Link](chi_table): Chi-squared approximation may be i
ncorrect
## CHI-SQUARE: Irrigation vs Food Shortage
## X-squared = 0.045 , df = 1 , p-value = 8.32e-01
##
## No Yes
## No 5 1
## Yes 20 9
##
## ============================================================
## LINEAR REGRESSION MODEL
## ============================================================
## Dependent Variable: Total Household Income (USD)
##
##
## Call:
## lm(formula = total_income ~ land_ha + irrigation + hh_size, data = r
eg_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5028.2 -2214.1 -204.7 1523.7 6057.8
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5453.0079 1610.4766 3.386 0.00205 **
## land_ha -1.3894 0.9524 -1.459 0.15535
## irrigationYes -1873.4659 1705.2814 -1.099 0.28097
## hh_size 75.3092 60.4576 1.246 0.22286
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3138 on 29 degrees of freedom
## Multiple R-squared: 0.1646, Adjusted R-squared: 0.07823
## F-statistic: 1.905 on 3 and 29 DF, p-value: 0.1508
##
##
## R-squared = 0.165
##
## ============================================================
## LOGISTIC REGRESSION MODEL
## ============================================================
## Dependent Variable: Food Shortage (1=Yes, 0=No)
##
##
## Call:
## glm(formula = food_shortage_binary ~ total_income + irrigation +
## drought, family = binomial, data = log_data)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.7122758 1.5386402 -1.113 0.266
## total_income -0.0001476 0.0001455 -1.014 0.310
## irrigationYes 0.8920794 1.2246657 0.728 0.466
## droughtYes 0.7774154 0.9316821 0.834 0.404
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 41.879 on 34 degrees of freedom
## Residual deviance: 39.088 on 31 degrees of freedom
## AIC: 47.088
##
## Number of Fisher Scoring iterations: 4
# ============================================================
# 9. QUALITATIVE CODING
# ============================================================
cat("\n============================================================\n")
##
## ============================================================
cat("QUALITATIVE ANALYSIS\n")
## QUALITATIVE ANALYSIS
cat("============================================================\n")
## ============================================================
# Extract challenges
challenges <- clean_data$challenge
challenges <- challenges[ & challenges != ""]
cat("Number of responses:", length(challenges), "\n")
## Number of responses: 36
if(length(challenges) > 0) {
# Coding function
code_challenge <- function(text) {
text_lower <- tolower(text)
if (grepl("pest|disease|armyworm|aphid|borer|locust", text_lower))
{
return("Pests/Diseases")
} else if (grepl("water|drought|irrigation|borehole|water scarcity",
text_lower)) {
return("Water Scarcity")
} else if (grepl("capital|loan|finance|credit|fund|resource", text_
lower)) {
return("Lack of Capital")
} else if (grepl("theft|security", text_lower)) {
return("Theft/Security")
} else if (grepl("climate|flood|climate change", text_lower)) {
return("Climate Change")
} else if (grepl("machinery|equipment|input|tool", text_lower)) {
return("Lack of Equipment")
} else if (grepl("market|price", text_lower)) {
return("Market Access")
} else if (grepl("extension|training|knowledge", text_lower)) {
return("Extension Gap")
} else if (grepl("electricity|power", text_lower)) {
return("Electricity")
} else if (grepl("soil fertility|soil", text_lower)) {
return("Soil Fertility")
} else {
return("Other")
}
}
# Apply coding
challenge_categories <- sapply(challenges, code_challenge)
# Frequency table
category_table <- table(challenge_categories)
category_df <- [Link](
Category = names(category_table),
Frequency = [Link](category_table),
Percentage = round([Link](category_table) / length(challenges)
* 100, 1)
)
category_df <- category_df[order(-category_df$Frequency), ]
print(category_df)
# Sample quotes
cat("\n=== SAMPLE QUOTES ===\n")
for(cat in head(unique(category_df$Category), 5)) {
cat("\n---", cat, "---\n")
quotes <- challenges[challenge_categories == cat]
if(length(quotes) > 0) {
cat(paste('"', quotes[1], '"', sep = ""), "\n")
}
}
# Save qualitative results
[Link](category_df, "outputs/qualitative_categories.csv", [Link]
s = FALSE)
}
## Category Frequency Percentage
## 8 Pests/Diseases 9 25.0
## 7 Other 7 19.4
## 4 Lack of Capital 6 16.7
## 5 Lack of Equipment 5 13.9
## 10 Theft/Security 2 5.6
## 11 Water Scarcity 2 5.6
## 1 Climate Change 1 2.8
## 2 Electricity 1 2.8
## 3 Extension Gap 1 2.8
## 6 Market Access 1 2.8
## 9 Soil Fertility 1 2.8
##
## === SAMPLE QUOTES ===
##
## --- Pests/Diseases ---
## "pests and diseases"
##
## --- Other ---
## "Predators"
##
## --- Lack of Capital ---
## "Capital"
##
## --- Lack of Equipment ---
## "Lack of Morden machinery"
##
## --- Theft/Security ---
## "Theft"
# ============================================================
# 10. FINAL SUMMARY
# ============================================================
cat("\n\n")
cat("============================================================\n")
## ============================================================
cat(" ANALYSIS COMPLETE \n")
## ANALYSIS COMPLETE
cat("============================================================\n")
## ============================================================
cat("\n✓ DATA FILES:\n")
##
## ✓ DATA FILES:
cat(" - cleaned_totonga_data.csv (cleaned dataset)\n")
## - cleaned_totonga_data.csv (cleaned dataset)
cat(" - [Link] (variable descriptions)\n")
## - [Link] (variable descriptions)
cat("\n✓ PLOTS FOLDER:\n")
##
## ✓ PLOTS FOLDER:
cat(" - income_distribution.png\n")
## - income_distribution.png
cat(" - food_shortage_by_irrigation.png\n")
## - food_shortage_by_irrigation.png
cat(" - income_by_education.png\n")
## - income_by_education.png
cat(" - time_allocation.png\n")
## - time_allocation.png
cat(" - income_by_gender.png\n")
## - income_by_gender.png
cat("\n✓ OUTPUTS FOLDER:\n")
##
## ✓ OUTPUTS FOLDER:
cat(" - descriptive_stats.csv\n")
## - descriptive_stats.csv
cat(" - linear_regression.csv (if enough data)\n")
## - linear_regression.csv (if enough data)
cat(" - logistic_regression.csv (if enough data)\n")
## - logistic_regression.csv (if enough data)
cat(" - qualitative_categories.csv\n")
## - qualitative_categories.csv
cat("\n✓ FINAL SAMPLE SIZE:", nrow(clean_data), "households\n")
##
## ✓ FINAL SAMPLE SIZE: 36 households
cat("============================================================\n")
## ============================================================