R Notebook
Stat 129 Final Project
Comparative Analysis of Rainfed and Irrigated Palay Volume of Production: A Multiple Linear Regression Approach
Macitas, G.I., Cayot, J.G., & Javiero, F.Y.
Division of Professional Education, College of Arts and Sciences, UP Visayas
Description
This study delves into the comparative analysis of rainfed and irrigated palay (rice) volume of production in the Philippines, employing a multiple linear regression approach. The research utilizes data extracted from the Philippine Statistics Authority’s 2012 Census of Agriculture and Fisheries. Palay production is a vital component of the Philippine agricultural sector, with both rainfed and irrigated cultivation methods playing significant roles in the country’s food security and economic stability. By employing multiple linear regression, this study aims to identify the key factors influencing palay production under rainfed and irrigated conditions, and to discern the relative impacts of these factors on production volume.
I. Load libraries
# Load the libraries
library(dplyr)
library(ggplot2)
library(corrplot)
library(caret)
library(dummy)
library(leaps)
library(car)
library(olsrr)
library(stargazer)
library(lmtest)
library(performance)
library(lm.beta)II. Load the Data
# Load the Data
data <- read.csv("Stat 129 Working Data.csv")
head(data)III. Rename the variables
# Rename the columns
data <- data %>%
rename(
palay = Palay.Total,
drought = Number.of.Barangays.that.Experienced.Drought.in.the.Past.Five.Years,
train_org = Number.of.Barangays.with.Organizations.Offering.Trainings.and.or.Extension.Services..Farmers..Fishermen.s.Organization.,
pest_dealers = Number.of.Barangays.with.Pesticides.Dealers,
wage_rate = Annual.Real.Wage.Rate.of.Palay.Workers,
hh_agri_act = Number.of.Household.Members.of.Agricultural.Operators.who.Engaged.in.Agricultural.Activities,
hh_over_10 = Number.of.Household.Members.of.Agricultural.Operators.who.were.10.Years.Old.and.Over,
infestation = Number.of.Barangays.that.Experienced.Infestation.Animal.Disease.in.the.Past.Five.Years,
rice_consumption = Annual.Per.Capita.Consumption.of.Rice,
island_group = Island.Group,
farm_area = Area.of.Farms.Under.Permanent.and.Temporary.Crops,
num_farms = Number.of.Farms.Under.Permanent.and.Temporary.Crops,
rainfed = Rainfed.Palay,
irrigated = Irrigated.Palay
)IV. Exploratory Data Analysis
summary(data)## X palay drought train_org
## Length:83 Min. : 0 Min. : 0.00 Min. : 3.0
## Class :character 1st Qu.: 51291 1st Qu.: 30.50 1st Qu.: 83.5
## Mode :character Median : 120953 Median : 66.00 Median :154.0
## Mean : 217259 Mean : 93.66 Mean :189.4
## 3rd Qu.: 291704 3rd Qu.:114.50 3rd Qu.:241.0
## Max. :1587163 Max. :691.00 Max. :743.0
##
## pest_dealers wage_rate hh_agri_act hh_over_10
## Min. : 2.00 Min. :139.0 Min. : 12 Min. : 8144
## 1st Qu.: 30.00 1st Qu.:169.1 1st Qu.: 503 1st Qu.: 118864
## Median : 55.00 Median :178.5 Median : 2444 Median : 215901
## Mean : 71.94 Mean :185.8 Mean : 7185 Mean : 376756
## 3rd Qu.: 93.00 3rd Qu.:203.2 3rd Qu.: 6941 3rd Qu.: 506914
## Max. :318.00 Max. :234.8 Max. :83597 Max. :1866445
##
## infestation rice_consumption island_group farm_area
## Min. : 0.0 Min. : 72.52 Length:83 Min. : 934
## 1st Qu.: 61.5 1st Qu.:110.86 Class :character 1st Qu.: 27597
## Median :113.0 Median :121.42 Mode :character Median : 60068
## Mean :138.8 Mean :120.23 Mean : 81808
## 3rd Qu.:180.0 3rd Qu.:130.86 3rd Qu.:116458
## Max. :676.0 Max. :164.59 Max. :328451
## NA's :1
## num_farms irrigated rainfed
## Min. : 4606 Min. : 0 Min. : 0
## 1st Qu.: 35475 1st Qu.: 24652 1st Qu.: 5270
## Median : 58070 Median : 85142 Median : 27683
## Mean : 68074 Mean : 161403 Mean : 55856
## 3rd Qu.: 86269 3rd Qu.: 170479 3rd Qu.: 84762
## Max. :198913 Max. :1452825 Max. :542556
## NA's :1# Check total production of irrigated vs. rainfed
irrigated_total <- sum(data$irrigated)
rainfed_total <- sum(data$rainfed)
print('Irrigated Palay Production in MT: ')## [1] "Irrigated Palay Production in MT: "irrigated_total## [1] 13396483print("Rainfed Palay Production in MT: ")## [1] "Rainfed Palay Production in MT: "rainfed_total## [1] 4636042# Histogram for Palay Production 2012
ggplot(data, aes(x = palay)) +
geom_histogram(binwidth = 100000, fill = "blue", color = "black", alpha = 0.5) +
labs(x = "Total Palay Production", y = "Frequency") +
ggtitle("Histogram of Palay Production 2012")
# Histogram for Irrigated Palay Production 2012
ggplot(data, aes(x = irrigated)) +
geom_histogram(binwidth = 100000, fill = "blue", color = "black", alpha = 0.5) +
labs(x = "Irrigated Palay Production", y = "Frequency") +
ggtitle("Histogram of Palay Production 2012")
# Histogram for Irrigated Palay Production 2012
ggplot(data, aes(x = rainfed)) +
geom_histogram(binwidth = 10000, fill = "blue", color = "black", alpha = 0.5) +
labs(x = "Rainfed Palay Production", y = "Frequency") +
ggtitle("Histogram of Palay Production 2012")
# Boxplot Batch Comparison
# Combine the data into a single dataframe
combined_data <- data.frame(
Production_Type = rep(c("Total", "Rainfed", "Irrigated"), each = nrow(data)),
Production_Value = c(data$palay, data$rainfed, data$irrigated)
)
# Boxplot
ggplot(combined_data, aes(x = Production_Type, y = Production_Value, fill = Production_Type)) +
geom_boxplot() +
labs(x = "Production Type", y = "Production (MT)", fill = "Production Type") +
ggtitle("Boxplot of Palay Production Types")
# Check missing values
missing_counts <- colSums(is.na(data))
total_missing <- sum(missing_counts)
print(total_missing)## [1] 2# Find rows with missing values
rows_with_missing <- which(!complete.cases(data))
# Display rows with missing values
print(data[rows_with_missing, ])## X palay drought train_org pest_dealers wage_rate hh_agri_act
## 69 ....City of Davao 15847 23 79 29 180.05 41715
## hh_over_10 infestation rice_consumption island_group farm_area num_farms
## 69 1293912 46 116.702 Mindanao NA NA
## irrigated rainfed
## 69 10275 5572# Impute missing values with column means
for (col in names(data)) {
if (anyNA(data[[col]])) {
data[[col]][is.na(data[[col]])] <- mean(data[[col]], na.rm = TRUE)
}
}
# Check missing values
missing_counts <- colSums(is.na(data))
total_missing <- sum(missing_counts)
print(total_missing)## [1] 0# Calculate correlation matrix
cor_matrix <- cor(data %>% select(- c(X, island_group)))
# Convert correlation matrix to tidy format for ggplot
cor_df <- reshape2::melt(cor_matrix)
# Plot heatmap using ggplot2
ggplot(cor_df, aes(x = Var1, y = Var2, fill = value)) +
geom_tile() +
scale_fill_gradient2(low = "blue", mid = "white", high = "red",
midpoint = 0, limits = c(-1, 1), name = "Correlation") +
labs(x = "", y = "") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
coord_fixed() # Optional: use fixed aspect ratio
V. Feature Engineering
# Create dummy variables for island_group
data$Luzon <- as.numeric(data$island_group == "Luzon")
data$Visayas <- as.numeric(data$island_group == "Visayas")
data[30:40,]# Create Area/Farm variable
data$area_per_farm <- data$farm_area/data$num_farmscolnames(data)## [1] "X" "palay" "drought" "train_org"
## [5] "pest_dealers" "wage_rate" "hh_agri_act" "hh_over_10"
## [9] "infestation" "rice_consumption" "island_group" "farm_area"
## [13] "num_farms" "irrigated" "rainfed" "Luzon"
## [17] "Visayas" "area_per_farm"data1 <- data_clean <- subset(data, select = -c(X, island_group))
colnames(data1)## [1] "palay" "drought" "train_org" "pest_dealers"
## [5] "wage_rate" "hh_agri_act" "hh_over_10" "infestation"
## [9] "rice_consumption" "farm_area" "num_farms" "irrigated"
## [13] "rainfed" "Luzon" "Visayas" "area_per_farm"VI. Split the dataset
# Split the dataset
# Set seed for reproducibility
set.seed(42)
# Split the data
trainIndex <- createDataPartition(data1$palay, p = .8,
list = FALSE,
times = 1)
data1Train <- data1[trainIndex,]
data1Test <- data1[-trainIndex,]
# Check dimensions of both sets
print(dim(data1Train))## [1] 67 16print(dim(data1Test))## [1] 16 16VII. Multiple Linear Regression
Total Palay Production
# Prepare the Training set for Palay
# Remove the irrigated and rainfed columns
palayTrain <- subset(data1Train, select = -c(irrigated, rainfed))
str(palayTrain)## 'data.frame': 67 obs. of 14 variables:
## $ palay : num 82333 105891 16696 62465 168854 ...
## $ drought : int 70 45 2 76 67 20 0 166 117 319 ...
## $ train_org : int 107 50 82 70 33 58 298 429 280 559 ...
## $ pest_dealers : int 29 39 55 33 26 27 119 93 71 220 ...
## $ wage_rate : num 179 179 179 179 179 ...
## $ hh_agri_act : int 4804 5715 5987 7219 8711 17028 7408 9783 6829 3362 ...
## $ hh_over_10 : int 109042 133703 103504 96817 84535 1144133 208417 228992 504579 1390728 ...
## $ infestation : int 104 74 44 107 103 62 0 172 120 262 ...
## $ rice_consumption: num 137 160 107 146 165 ...
## $ farm_area : num 17454 21370 17372 20825 27533 ...
## $ num_farms : num 46866 25505 52940 38237 37547 ...
## $ Luzon : num 1 1 1 1 1 1 1 1 1 1 ...
## $ Visayas : num 0 0 0 0 0 0 0 0 0 0 ...
## $ area_per_farm : num 0.372 0.838 0.328 0.545 0.733 ...# Create the full model
palaymodel1 <- lm(palay ~ ., data = palayTrain)
# Perform ANOVA F-test
summary(palaymodel1)##
## Call:
## lm(formula = palay ~ ., data = palayTrain)
##
## Residuals:
## Min 1Q Median 3Q Max
## -196994 -43085 -12258 41933 224962
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.649e+05 1.653e+05 -0.997 0.323152
## drought 1.011e+03 2.644e+02 3.823 0.000349 ***
## train_org -2.311e+02 2.025e+02 -1.141 0.258878
## pest_dealers 5.260e+03 4.116e+02 12.781 < 2e-16 ***
## wage_rate -3.623e+02 7.040e+02 -0.515 0.608914
## hh_agri_act -3.398e+00 1.343e+00 -2.529 0.014443 *
## hh_over_10 4.919e-02 2.984e-02 1.648 0.105223
## infestation -7.741e+02 2.739e+02 -2.826 0.006637 **
## rice_consumption 1.933e+03 6.834e+02 2.829 0.006583 **
## farm_area 2.593e-01 8.000e-01 0.324 0.747125
## num_farms -1.116e+00 1.089e+00 -1.025 0.310013
## Luzon -9.035e+04 3.281e+04 -2.754 0.008057 **
## Visayas -2.415e+04 4.572e+04 -0.528 0.599617
## area_per_farm 3.243e+03 4.805e+04 0.067 0.946452
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 86430 on 53 degrees of freedom
## Multiple R-squared: 0.9361, Adjusted R-squared: 0.9205
## F-statistic: 59.75 on 13 and 53 DF, p-value: < 2.2e-16# Both-direction stepwise regression and obtain the model with least AIC
stepwise.palay.model <- step(palaymodel1, direction = "both")## Start: AIC=1535.49
## palay ~ drought + train_org + pest_dealers + wage_rate + hh_agri_act +
## hh_over_10 + infestation + rice_consumption + farm_area +
## num_farms + Luzon + Visayas + area_per_farm
##
## Df Sum of Sq RSS AIC
## - area_per_farm 1 3.4019e+07 3.9597e+11 1533.5
## - farm_area 1 7.8479e+08 3.9672e+11 1533.6
## - wage_rate 1 1.9789e+09 3.9791e+11 1533.8
## - Visayas 1 2.0836e+09 3.9802e+11 1533.8
## - num_farms 1 7.8489e+09 4.0378e+11 1534.8
## - train_org 1 9.7307e+09 4.0567e+11 1535.1
## <none> 3.9594e+11 1535.5
## - hh_over_10 1 2.0295e+10 4.1623e+11 1536.8
## - hh_agri_act 1 4.7788e+10 4.4372e+11 1541.1
## - Luzon 1 5.6642e+10 4.5258e+11 1542.5
## - infestation 1 5.9652e+10 4.5559e+11 1542.9
## - rice_consumption 1 5.9779e+10 4.5571e+11 1542.9
## - drought 1 1.0916e+11 5.0509e+11 1549.8
## - pest_dealers 1 1.2203e+12 1.6162e+12 1627.7
##
## Step: AIC=1533.49
## palay ~ drought + train_org + pest_dealers + wage_rate + hh_agri_act +
## hh_over_10 + infestation + rice_consumption + farm_area +
## num_farms + Luzon + Visayas
##
## Df Sum of Sq RSS AIC
## - wage_rate 1 1.9511e+09 3.9792e+11 1531.8
## - Visayas 1 2.0539e+09 3.9802e+11 1531.8
## - farm_area 1 4.6943e+09 4.0066e+11 1532.3
## - train_org 1 1.1125e+10 4.0709e+11 1533.3
## <none> 3.9597e+11 1533.5
## - num_farms 1 1.7817e+10 4.1379e+11 1534.4
## - hh_over_10 1 2.0570e+10 4.1654e+11 1534.9
## + area_per_farm 1 3.4019e+07 3.9594e+11 1535.5
## - hh_agri_act 1 4.8784e+10 4.4475e+11 1539.3
## - Luzon 1 5.7953e+10 4.5392e+11 1540.7
## - rice_consumption 1 6.0621e+10 4.5659e+11 1541.0
## - infestation 1 6.4320e+10 4.6029e+11 1541.6
## - drought 1 1.0962e+11 5.0559e+11 1547.9
## - pest_dealers 1 1.2430e+12 1.6390e+12 1626.7
##
## Step: AIC=1531.82
## palay ~ drought + train_org + pest_dealers + hh_agri_act + hh_over_10 +
## infestation + rice_consumption + farm_area + num_farms +
## Luzon + Visayas
##
## Df Sum of Sq RSS AIC
## - Visayas 1 6.6527e+08 3.9859e+11 1529.9
## - farm_area 1 5.7201e+09 4.0364e+11 1530.8
## <none> 3.9792e+11 1531.8
## - train_org 1 1.2342e+10 4.1026e+11 1531.9
## - num_farms 1 1.6312e+10 4.1423e+11 1532.5
## - hh_over_10 1 1.9496e+10 4.1742e+11 1533.0
## + wage_rate 1 1.9511e+09 3.9597e+11 1533.5
## + area_per_farm 1 6.2021e+06 3.9791e+11 1533.8
## - hh_agri_act 1 4.7079e+10 4.4500e+11 1537.3
## - infestation 1 6.2776e+10 4.6070e+11 1539.6
## - rice_consumption 1 6.3406e+10 4.6133e+11 1539.7
## - Luzon 1 6.8332e+10 4.6625e+11 1540.4
## - drought 1 1.0772e+11 5.0564e+11 1545.9
## - pest_dealers 1 1.5537e+12 1.9516e+12 1636.4
##
## Step: AIC=1529.94
## palay ~ drought + train_org + pest_dealers + hh_agri_act + hh_over_10 +
## infestation + rice_consumption + farm_area + num_farms +
## Luzon
##
## Df Sum of Sq RSS AIC
## - farm_area 1 8.3979e+09 4.0698e+11 1529.3
## <none> 3.9859e+11 1529.9
## - train_org 1 1.2626e+10 4.1121e+11 1530.0
## - num_farms 1 1.6484e+10 4.1507e+11 1530.7
## - hh_over_10 1 1.8937e+10 4.1752e+11 1531.0
## + Visayas 1 6.6527e+08 3.9792e+11 1531.8
## + wage_rate 1 5.6244e+08 3.9802e+11 1531.8
## + area_per_farm 1 3.8083e+06 3.9858e+11 1531.9
## - hh_agri_act 1 4.7090e+10 4.4568e+11 1535.4
## - rice_consumption 1 6.4445e+10 4.6303e+11 1538.0
## - infestation 1 6.9430e+10 4.6802e+11 1538.7
## - Luzon 1 7.6728e+10 4.7531e+11 1539.7
## - drought 1 1.0794e+11 5.0653e+11 1544.0
## - pest_dealers 1 1.5571e+12 1.9557e+12 1634.5
##
## Step: AIC=1529.33
## palay ~ drought + train_org + pest_dealers + hh_agri_act + hh_over_10 +
## infestation + rice_consumption + num_farms + Luzon
##
## Df Sum of Sq RSS AIC
## - num_farms 1 8.2709e+09 4.1525e+11 1528.7
## <none> 4.0698e+11 1529.3
## + farm_area 1 8.3979e+09 3.9859e+11 1529.9
## + area_per_farm 1 6.4504e+09 4.0053e+11 1530.3
## - hh_over_10 1 1.9357e+10 4.2634e+11 1530.5
## + Visayas 1 3.3430e+09 4.0364e+11 1530.8
## + wage_rate 1 2.6948e+08 4.0671e+11 1531.3
## - train_org 1 4.0506e+10 4.4749e+11 1533.7
## - rice_consumption 1 5.7339e+10 4.6432e+11 1536.2
## - infestation 1 6.1151e+10 4.6813e+11 1536.7
## - hh_agri_act 1 6.5762e+10 4.7275e+11 1537.4
## - Luzon 1 8.7028e+10 4.9401e+11 1540.3
## - drought 1 1.0310e+11 5.1008e+11 1542.5
## - pest_dealers 1 1.6969e+12 2.1039e+12 1637.4
##
## Step: AIC=1528.68
## palay ~ drought + train_org + pest_dealers + hh_agri_act + hh_over_10 +
## infestation + rice_consumption + Luzon
##
## Df Sum of Sq RSS AIC
## <none> 4.1525e+11 1528.7
## + num_farms 1 8.2709e+09 4.0698e+11 1529.3
## + area_per_farm 1 5.2235e+09 4.1003e+11 1529.8
## - hh_over_10 1 2.0852e+10 4.3611e+11 1530.0
## + Visayas 1 1.0219e+09 4.1423e+11 1530.5
## + farm_area 1 1.8456e+08 4.1507e+11 1530.7
## + wage_rate 1 9.4186e+06 4.1524e+11 1530.7
## - rice_consumption 1 6.2024e+10 4.7728e+11 1536.0
## - train_org 1 7.3159e+10 4.8841e+11 1537.5
## - infestation 1 7.8820e+10 4.9407e+11 1538.3
## - Luzon 1 8.0598e+10 4.9585e+11 1538.6
## - hh_agri_act 1 1.0178e+11 5.1704e+11 1541.4
## - drought 1 1.2039e+11 5.3565e+11 1543.7
## - pest_dealers 1 2.2857e+12 2.7010e+12 1652.1stepwise.palay.model##
## Call:
## lm(formula = palay ~ drought + train_org + pest_dealers + hh_agri_act +
## hh_over_10 + infestation + rice_consumption + Luzon, data = palayTrain)
##
## Coefficients:
## (Intercept) drought train_org pest_dealers
## -2.301e+05 1.020e+03 -4.207e+02 5.030e+03
## hh_agri_act hh_over_10 infestation rice_consumption
## -4.238e+00 4.855e-02 -7.457e+02 1.889e+03
## Luzon
## -9.083e+04palaymodel2 <- lm(palay ~ .-area_per_farm, data = palayTrain)
palaymodel3 <- lm(palay ~ .-area_per_farm -wage_rate, data = palayTrain)
palaymodel4 <- lm(palay ~ .-area_per_farm -wage_rate -Visayas, data = palayTrain)
palaymodel5 <- lm(palay ~ .-area_per_farm -wage_rate -Visayas -farm_area, data = palayTrain)
bestpalaymodel <- lm(palay ~ .-area_per_farm -wage_rate -Visayas -farm_area -num_farms, data = palayTrain)# Summarize the models in one table
palaymodels <- list(palaymodel1,palaymodel2,palaymodel3, palaymodel4, palaymodel5, bestpalaymodel)
resultspalay <- data.frame(stargazer(palaymodels,
type = "text",
title = "Regression Model Results",
dep.var.labels = "Palay Production",
covariate.labels = c("Drought", "Training Organizations", "Pesticide Dealers",
"Wage Rate", "HH Agricultural Activities", "HH Over 10",
"Infestation", "Rice Consumption", "Farm Area",
"Number of Farms", "Luzon", "Visayas", "Area per Farm"),
omit.stat = c("f", "ser"),
align = TRUE))##
## Regression Model Results
## =======================================================================================================================
## Dependent variable:
## --------------------------------------------------------------------------------------------
## Palay Production
## (1) (2) (3) (4) (5) (6)
## -----------------------------------------------------------------------------------------------------------------------
## Drought 1,010.576*** 1,011.465*** 990.735*** 991.675*** 964.500*** 1,020.292***
## (264.370) (261.596) (256.760) (254.651) (253.819) (248.808)
##
## Training Organizations -231.090 -225.505 -236.037 -238.495 -350.007** -420.731***
## (202.481) (183.080) (180.720) (179.068) (146.949) (131.618)
##
## Pesticide Dealers 5,260.111*** 5,256.209*** 5,156.322*** 5,148.665*** 5,234.511*** 5,029.533***
## (411.563) (403.708) (351.864) (348.100) (339.544) (281.486)
##
## Wage Rate -362.321 -353.480
## (703.970) (685.264)
##
## HH Agricultural Activities -3.398** -3.382** -3.295** -3.295** -3.714*** -4.238***
## (1.343) (1.311) (1.292) (1.281) (1.224) (1.124)
##
## HH Over 10 0.049 0.049* 0.047 0.046 0.047 0.049*
## (0.030) (0.029) (0.029) (0.028) (0.028) (0.028)
##
## Infestation -774.086*** -778.672*** -744.600*** -762.184*** -680.421*** -745.724***
## (273.938) (262.916) (252.780) (244.035) (232.502) (224.752)
##
## Rice Consumption 1,933.235*** 1,937.690*** 1,972.021*** 1,984.387*** 1,824.657*** 1,889.434***
## (683.415) (673.920) (666.136) (659.474) (643.884) (641.939)
##
## Farm Area 0.259 0.307 0.335 0.377
## (0.800) (0.383) (0.377) (0.347)
##
## Number of Farms -1.116 -1.169 -1.101 -1.107 -0.608
## (1.089) (0.750) (0.733) (0.727) (0.565)
##
## Luzon -90,350.350*** -90,634.540*** -94,986.550*** -90,976.900*** -95,693.920*** -90,830.900***
## (32,812.250) (32,239.560) (30,907.730) (27,709.040) (27,409.780) (27,071.690)
##
## Visayas -24,146.920 -23,892.780 -11,519.690
## (45,722.040) (45,144.730) (37,988.970)
##
## Area per Farm 3,242.645
## (48,051.790)
##
## Constant -164,897.800 -163,528.200 -234,349.900*** -238,927.100*** -213,895.600** -230,097.100***
## (165,344.600) (162,574.700) (86,481.130) (84,460.430) (81,383.930) (80,088.770)
##
## -----------------------------------------------------------------------------------------------------------------------
## Observations 67 67 67 67 67 67
## R2 0.936 0.936 0.936 0.936 0.934 0.933
## Adjusted R2 0.920 0.922 0.923 0.924 0.924 0.924
## =======================================================================================================================
## Note: *p<0.1; **p<0.05; ***p<0.01palayTest <- subset(data1Test, select = -c(irrigated, rainfed))
test_palay <- subset(palayTest, select = -c(palay))
predictions_palay1 <- predict(palaymodel1, newdata = test_palay)
predictions_palay2 <- predict(palaymodel2, newdata = test_palay)
predictions_palay3 <- predict(palaymodel3, newdata = test_palay)
predictions_palay4 <- predict(palaymodel4, newdata = test_palay)
predictions_palay5 <- predict(palaymodel5, newdata = test_palay)
predictions_bestpalaymodel <- predict(bestpalaymodel, newdata = test_palay)
actual_palay <- palayTest$palay
rmsepalay1 <- sqrt(mean((predictions_palay1 - actual_palay)^2))
rmsepalay2 <- sqrt(mean((predictions_palay2 - actual_palay)^2))
rmsepalay3 <- sqrt(mean((predictions_palay3 - actual_palay)^2))
rmsepalay4 <- sqrt(mean((predictions_palay4 - actual_palay)^2))
rmsepalay5 <- sqrt(mean((predictions_palay5 - actual_palay)^2))
rmsepalay6 <- sqrt(mean((predictions_bestpalaymodel - actual_palay)^2))
rownames <- c('Model1', 'Model2', 'Model3', 'Model4', 'Model5', 'Model6')
rmsepalay <- data.frame(row.names = rownames, RMSE = c(rmsepalay1,rmsepalay2,rmsepalay3,rmsepalay4,rmsepalay5,rmsepalay6))
rmsepalayBest on the results, Model 5 better predicts the model better. Let’s check their other metrics. We will adhere to the principle of parsimony.
# R-squared and Adjusted R-squared
rsquared5 <- summary(palaymodel5)$r.squared
adjrsquared5 <- summary(palaymodel5)$adj.r.squared
rsquared6 <- summary(bestpalaymodel)$r.squared
adjrsquared6 <- summary(bestpalaymodel)$adj.r.squared
# AIC and BIC
aic5 <- AIC(palaymodel5)
bic5 <- BIC(palaymodel5)
aic6 <- AIC(bestpalaymodel)
bic6 <- BIC(bestpalaymodel)
# Mallows' Cp
cp5 <- ols_mallows_cp(palaymodel5, palaymodel1)
cp6 <- ols_mallows_cp(bestpalaymodel, palaymodel1)
model_comparison <- data.frame(
Model = c("Model 5", "Best Total Palay Model"),
R_squared = c(rsquared5, rsquared6),
Adjusted_R_squared = c(adjrsquared5, adjrsquared6),
AIC = c(aic5, aic6),
BIC = c(bic5, bic6),
Mallows_Cp = c(cp5, cp6),
RMSE = c(rmsepalay5, rmsepalay6)
)
print(model_comparison)## Model R_squared Adjusted_R_squared AIC BIC
## 1 Model 5 0.9343461 0.9239797 1721.470 1745.722
## 2 Best Total Palay Model 0.9330119 0.9237722 1720.818 1742.865
## Mallows_Cp RMSE
## 1 7.478926 96909.00
## 2 6.586070 96920.51While Model 5 has higher R-square and Adjusted R-square, Model 6 has the lower AIC, BIC and Mallows’ Cp. We can run a general linear test approach to determine if the full model (Model 5) is better than the reduced model (Model 6).
# General Linear Test Approach
anova_result1 <- anova(bestpalaymodel, palaymodel5)
anova_result1Since p-value > 0.05, we fail to reject the null hypothesis that the full model is better than the reduced model. Thus, for the total palay production, we adopt Model 6. Thus our final model for the Palay Production is given by:
bestpalaymodel$coefficients## (Intercept) drought train_org pest_dealers
## -2.300971e+05 1.020292e+03 -4.207310e+02 5.029533e+03
## hh_agri_act hh_over_10 infestation rice_consumption
## -4.238047e+00 4.854889e-02 -7.457244e+02 1.889434e+03
## Luzon
## -9.083090e+04# Calculate standardized coefficients
standardized_palaymodel <- lm.beta(bestpalaymodel)
summary(standardized_palaymodel)##
## Call:
## lm(formula = palay ~ . - area_per_farm - wage_rate - Visayas -
## farm_area - num_farms, data = palayTrain)
##
## Residuals:
## Min 1Q Median 3Q Max
## -169767 -54188 -10901 40616 239191
##
## Coefficients:
## Estimate Standardized Std. Error t value Pr(>|t|)
## (Intercept) -2.301e+05 NA 8.009e+04 -2.873 0.005672 **
## drought 1.020e+03 3.871e-01 2.488e+02 4.101 0.000130 ***
## train_org -4.207e+02 -2.066e-01 1.316e+02 -3.197 0.002252 **
## pest_dealers 5.030e+03 1.096e+00 2.815e+02 17.868 < 2e-16 ***
## hh_agri_act -4.238e+00 -1.905e-01 1.124e+00 -3.770 0.000384 ***
## hh_over_10 4.855e-02 6.298e-02 2.845e-02 1.707 0.093246 .
## infestation -7.457e+02 -2.967e-01 2.248e+02 -3.318 0.001570 **
## rice_consumption 1.889e+03 1.131e-01 6.419e+02 2.943 0.004663 **
## Luzon -9.083e+04 -1.492e-01 2.707e+04 -3.355 0.001403 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 84610 on 58 degrees of freedom
## Multiple R-squared: 0.933, Adjusted R-squared: 0.9238
## F-statistic: 101 on 8 and 58 DF, p-value: < 2.2e-16# Extract the standardized coefficients
standardized_palaycoefficients <- summary(standardized_palaymodel)$coefficients
# Create a data frame for plotting
coefficients_df <- as.data.frame(standardized_palaycoefficients)
coefficients_df$Predictor <- rownames(coefficients_df)
rownames(coefficients_df) <- NULL# Extract only the Estimate column and rename for clarity
coefficients_df <- coefficients_df[, c("Predictor", "Standardized")]
names(coefficients_df)[2] <- "Estimate"
# Keep only the coefficients for predictors (ignore the intercept)
coefficients_df <- coefficients_df[-1, ]
# Plot the standardized coefficients
ggplot(coefficients_df, aes(x = reorder(Predictor, Estimate), y = Estimate, fill=Predictor)) +
geom_bar(stat = "identity", show.legend = FALSE) +
coord_flip() + # Flip coordinates to make it horizontal
xlab("Predictor") +
ylab("Standardized Coefficient") +
ggtitle("Predictor Importance in Total Palay Production Model") 
Let’s perform diagnostic checking.
plot(bestpalaymodel)
High-leverage Observation
For the residual plot, only 1 observation is considered a high-leverage observation. That is observation 19 or the province of Nueva Ecija. This is not surprising since Nueva Ecija is known to be a major producer of rice in the country and is often called the “Rice Granary of the Philippines.” We opt not to remove this case in the model.
Multicollinearity
We can check the VIF of variables to see if there’s a presence of multicollinearity.
# Compute VIF for the model
vif_values <- vif(bestpalaymodel)
# Print VIF values
print(vif_values)## drought train_org pest_dealers hh_agri_act
## 7.714330 3.617180 3.260340 2.211310
## hh_over_10 infestation rice_consumption Luzon
## 1.178994 6.922649 1.278965 1.711148Normality
Let’s check for the normality of the residuals.
# Perform Shapiro-Wilk test on residuals
shapiro.test(residuals(bestpalaymodel))##
## Shapiro-Wilk normality test
##
## data: residuals(bestpalaymodel)
## W = 0.97056, p-value = 0.1127Thus, we fail to reject the null hypothesis that the model is normal.
Homoscedasticity
# Perform Breusch-Pagan Test
bptest(bestpalaymodel)##
## studentized Breusch-Pagan test
##
## data: bestpalaymodel
## BP = 25.661, df = 8, p-value = 0.0012Independence of Errors
# Perform Durbin-Watson
dwtest(bestpalaymodel)##
## Durbin-Watson test
##
## data: bestpalaymodel
## DW = 2.305, p-value = 0.8355
## alternative hypothesis: true autocorrelation is greater than 0Thus, there is no serial autocorrelation in the model.
Linearity
plot(bestpalaymodel, which = 1)
abline(h = 0, col = "red", lty = 2) # Add horizontal line at y = 0
Rainfed Palay Production
# Prepare the Training set for Rainfed Palays
# Remove the irrigated and rainfed columns
rainfedTrain <- subset(data1Train, select = -c(palay, irrigated ))
str(rainfedTrain)## 'data.frame': 67 obs. of 14 variables:
## $ drought : int 70 45 2 76 67 20 0 166 117 319 ...
## $ train_org : int 107 50 82 70 33 58 298 429 280 559 ...
## $ pest_dealers : int 29 39 55 33 26 27 119 93 71 220 ...
## $ wage_rate : num 179 179 179 179 179 ...
## $ hh_agri_act : int 4804 5715 5987 7219 8711 17028 7408 9783 6829 3362 ...
## $ hh_over_10 : int 109042 133703 103504 96817 84535 1144133 208417 228992 504579 1390728 ...
## $ infestation : int 104 74 44 107 103 62 0 172 120 262 ...
## $ rice_consumption: num 137 160 107 146 165 ...
## $ farm_area : num 17454 21370 17372 20825 27533 ...
## $ num_farms : num 46866 25505 52940 38237 37547 ...
## $ rainfed : num 27114 29571 2136 1805 4893 ...
## $ Luzon : num 1 1 1 1 1 1 1 1 1 1 ...
## $ Visayas : num 0 0 0 0 0 0 0 0 0 0 ...
## $ area_per_farm : num 0.372 0.838 0.328 0.545 0.733 ...# Create the full model
rainfedmodel1 <- lm(rainfed ~ ., data = rainfedTrain)
# Perform ANOVA F-test
summary(rainfedmodel1)##
## Call:
## lm(formula = rainfed ~ ., data = rainfedTrain)
##
## Residuals:
## Min 1Q Median 3Q Max
## -75547 -18544 652 18058 79430
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.693e+05 6.775e+04 -2.499 0.015600 *
## drought 6.342e+02 1.083e+02 5.855 3.08e-07 ***
## train_org 2.224e+01 8.296e+01 0.268 0.789722
## pest_dealers 3.327e+02 1.686e+02 1.973 0.053693 .
## wage_rate 9.119e+01 2.884e+02 0.316 0.753120
## hh_agri_act -2.775e+00 5.504e-01 -5.040 5.76e-06 ***
## hh_over_10 3.202e-02 1.223e-02 2.618 0.011492 *
## infestation -1.263e+02 1.122e+02 -1.125 0.265529
## rice_consumption 1.142e+03 2.800e+02 4.078 0.000153 ***
## farm_area 7.022e-02 3.278e-01 0.214 0.831181
## num_farms 1.413e-01 4.462e-01 0.317 0.752801
## Luzon -3.574e+04 1.344e+04 -2.658 0.010361 *
## Visayas -1.389e+04 1.873e+04 -0.741 0.461688
## area_per_farm 8.466e+03 1.969e+04 0.430 0.668933
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 35410 on 53 degrees of freedom
## Multiple R-squared: 0.8576, Adjusted R-squared: 0.8227
## F-statistic: 24.56 on 13 and 53 DF, p-value: < 2.2e-16# Both-direction stepwise regression and obtain the model with least AIC
stepwise.rainfed.model <- step(rainfedmodel1, direction = "both")## Start: AIC=1415.93
## rainfed ~ drought + train_org + pest_dealers + wage_rate + hh_agri_act +
## hh_over_10 + infestation + rice_consumption + farm_area +
## num_farms + Luzon + Visayas + area_per_farm
##
## Df Sum of Sq RSS AIC
## - farm_area 1 5.7564e+07 6.6527e+10 1414.0
## - train_org 1 9.0093e+07 6.6560e+10 1414.0
## - wage_rate 1 1.2536e+08 6.6595e+10 1414.0
## - num_farms 1 1.2570e+08 6.6595e+10 1414.0
## - area_per_farm 1 2.3190e+08 6.6702e+10 1414.2
## - Visayas 1 6.8948e+08 6.7159e+10 1414.6
## - infestation 1 1.5881e+09 6.8058e+10 1415.5
## <none> 6.6470e+10 1415.9
## - pest_dealers 1 4.8832e+09 7.1353e+10 1418.7
## - hh_over_10 1 8.5986e+09 7.5068e+10 1422.1
## - Luzon 1 8.8623e+09 7.5332e+10 1422.3
## - rice_consumption 1 2.0860e+10 8.7330e+10 1432.2
## - hh_agri_act 1 3.1862e+10 9.8332e+10 1440.2
## - drought 1 4.2995e+10 1.0946e+11 1447.3
##
## Step: AIC=1413.98
## rainfed ~ drought + train_org + pest_dealers + wage_rate + hh_agri_act +
## hh_over_10 + infestation + rice_consumption + num_farms +
## Luzon + Visayas + area_per_farm
##
## Df Sum of Sq RSS AIC
## - train_org 1 3.9326e+07 6.6567e+10 1412.0
## - wage_rate 1 9.4212e+07 6.6621e+10 1412.1
## - Visayas 1 8.5289e+08 6.7380e+10 1412.8
## - num_farms 1 9.0632e+08 6.7433e+10 1412.9
## - infestation 1 1.5854e+09 6.8113e+10 1413.6
## <none> 6.6527e+10 1414.0
## - area_per_farm 1 2.0469e+09 6.8574e+10 1414.0
## + farm_area 1 5.7564e+07 6.6470e+10 1415.9
## - pest_dealers 1 5.5135e+09 7.2041e+10 1417.3
## - Luzon 1 8.8656e+09 7.5393e+10 1420.4
## - hh_over_10 1 9.1303e+09 7.5657e+10 1420.6
## - rice_consumption 1 2.1276e+10 8.7803e+10 1430.6
## - hh_agri_act 1 3.5633e+10 1.0216e+11 1440.7
## - drought 1 4.3012e+10 1.0954e+11 1445.4
##
## Step: AIC=1412.02
## rainfed ~ drought + pest_dealers + wage_rate + hh_agri_act +
## hh_over_10 + infestation + rice_consumption + num_farms +
## Luzon + Visayas + area_per_farm
##
## Df Sum of Sq RSS AIC
## - wage_rate 1 1.2450e+08 6.6691e+10 1410.2
## - Visayas 1 8.1515e+08 6.7382e+10 1410.8
## - num_farms 1 1.4773e+09 6.8044e+10 1411.5
## - infestation 1 1.5482e+09 6.8115e+10 1411.6
## <none> 6.6567e+10 1412.0
## - area_per_farm 1 2.0290e+09 6.8596e+10 1412.0
## + train_org 1 3.9326e+07 6.6527e+10 1414.0
## + farm_area 1 6.7972e+06 6.6560e+10 1414.0
## - pest_dealers 1 5.4809e+09 7.2047e+10 1415.3
## - Luzon 1 8.8290e+09 7.5396e+10 1418.4
## - hh_over_10 1 9.5045e+09 7.6071e+10 1419.0
## - rice_consumption 1 2.2093e+10 8.8659e+10 1429.2
## - drought 1 4.5445e+10 1.1201e+11 1444.9
## - hh_agri_act 1 5.0722e+10 1.1729e+11 1448.0
##
## Step: AIC=1410.15
## rainfed ~ drought + pest_dealers + hh_agri_act + hh_over_10 +
## infestation + rice_consumption + num_farms + Luzon + Visayas +
## area_per_farm
##
## Df Sum of Sq RSS AIC
## - num_farms 1 1.3555e+09 6.8047e+10 1409.5
## - Visayas 1 1.4970e+09 6.8188e+10 1409.6
## - infestation 1 1.9530e+09 6.8644e+10 1410.1
## - area_per_farm 1 1.9585e+09 6.8649e+10 1410.1
## <none> 6.6691e+10 1410.2
## + wage_rate 1 1.2450e+08 6.6567e+10 1412.0
## + train_org 1 6.9614e+07 6.6621e+10 1412.1
## + farm_area 1 7.2284e+05 6.6690e+10 1412.2
## - pest_dealers 1 7.9613e+09 7.4652e+10 1415.7
## - Luzon 1 9.0349e+09 7.5726e+10 1416.7
## - hh_over_10 1 9.9350e+09 7.6626e+10 1417.5
## - rice_consumption 1 2.2023e+10 8.8714e+10 1427.3
## - drought 1 4.8680e+10 1.1537e+11 1444.9
## - hh_agri_act 1 5.5079e+10 1.2177e+11 1448.5
##
## Step: AIC=1409.5
## rainfed ~ drought + pest_dealers + hh_agri_act + hh_over_10 +
## infestation + rice_consumption + Luzon + Visayas + area_per_farm
##
## Df Sum of Sq RSS AIC
## - infestation 1 1.0439e+09 6.9090e+10 1408.5
## - area_per_farm 1 1.5462e+09 6.9593e+10 1409.0
## <none> 6.8047e+10 1409.5
## - Visayas 1 2.1482e+09 7.0195e+10 1409.6
## + num_farms 1 1.3555e+09 6.6691e+10 1410.2
## + farm_area 1 6.8923e+08 6.7357e+10 1410.8
## + train_org 1 6.1285e+08 6.7434e+10 1410.9
## + wage_rate 1 2.6699e+06 6.8044e+10 1411.5
## - hh_over_10 1 1.0408e+10 7.8455e+10 1417.0
## - Luzon 1 1.1009e+10 7.9056e+10 1417.5
## - rice_consumption 1 2.0669e+10 8.8715e+10 1425.3
## - pest_dealers 1 2.0745e+10 8.8791e+10 1425.3
## - drought 1 4.7633e+10 1.1568e+11 1443.0
## - hh_agri_act 1 5.3805e+10 1.2185e+11 1446.5
##
## Step: AIC=1408.52
## rainfed ~ drought + pest_dealers + hh_agri_act + hh_over_10 +
## rice_consumption + Luzon + Visayas + area_per_farm
##
## Df Sum of Sq RSS AIC
## - area_per_farm 1 1.3230e+09 7.0413e+10 1407.8
## <none> 6.9090e+10 1408.5
## - Visayas 1 2.6064e+09 7.1697e+10 1409.0
## + infestation 1 1.0439e+09 6.8047e+10 1409.5
## + num_farms 1 4.4643e+08 6.8644e+10 1410.1
## + wage_rate 1 1.9071e+08 6.8900e+10 1410.3
## + train_org 1 1.6948e+08 6.8921e+10 1410.3
## + farm_area 1 6.3215e+07 6.9027e+10 1410.5
## - Luzon 1 9.9818e+09 7.9072e+10 1415.6
## - hh_over_10 1 1.1248e+10 8.0339e+10 1416.6
## - pest_dealers 1 1.9747e+10 8.8838e+10 1423.4
## - rice_consumption 1 2.0330e+10 8.9421e+10 1423.8
## - hh_agri_act 1 5.5330e+10 1.2442e+11 1445.9
## - drought 1 9.9529e+10 1.6862e+11 1466.3
##
## Step: AIC=1407.79
## rainfed ~ drought + pest_dealers + hh_agri_act + hh_over_10 +
## rice_consumption + Luzon + Visayas
##
## Df Sum of Sq RSS AIC
## <none> 7.0413e+10 1407.8
## + area_per_farm 1 1.3230e+09 6.9090e+10 1408.5
## + farm_area 1 9.4664e+08 6.9467e+10 1408.9
## + infestation 1 8.2081e+08 6.9593e+10 1409.0
## + num_farms 1 3.0620e+08 7.0107e+10 1409.5
## + wage_rate 1 1.1753e+08 7.0296e+10 1409.7
## + train_org 1 1.3609e+07 7.0400e+10 1409.8
## - Visayas 1 4.9019e+09 7.5315e+10 1410.3
## - hh_over_10 1 1.0667e+10 8.1080e+10 1415.2
## - Luzon 1 1.4892e+10 8.5305e+10 1418.6
## - rice_consumption 1 1.9701e+10 9.0114e+10 1422.3
## - pest_dealers 1 2.1425e+10 9.1838e+10 1423.6
## - hh_agri_act 1 5.8115e+10 1.2853e+11 1446.1
## - drought 1 9.8747e+10 1.6916e+11 1464.5rainfedmodel2 <- lm(rainfed ~ .-farm_area, data = rainfedTrain)
rainfedmodel3 <- lm(rainfed ~ .-train_org -farm_area, data = rainfedTrain)
rainfedmodel4 <- lm(rainfed ~ .-wage_rate -train_org -farm_area, data = rainfedTrain)
rainfedmodel5 <- lm(rainfed ~ .-num_farms -farm_area -train_org -wage_rate, data = rainfedTrain)
rainfedmodel6 <- lm(rainfed ~ .-infestation -num_farms -wage_rate -train_org -farm_area, data = rainfedTrain)
bestrainfedmodel <- lm(rainfed ~ .-area_per_farm -farm_area -infestation -num_farms -wage_rate -train_org, data = rainfedTrain)# Summarize the models in one table
rainfedmodels <- list(rainfedmodel1,rainfedmodel2,rainfedmodel3, rainfedmodel4, rainfedmodel5, rainfedmodel6, bestrainfedmodel)
resultsrainfed <- data.frame(stargazer(rainfedmodels,
type = "text",
title = "Regression Model Results",
dep.var.labels = "Rainfed Palay Production",
covariate.labels = c("Drought", "Training Organizations", "Pesticide Dealers",
"Wage Rate", "HH Agricultural Activities", "HH Over 10",
"Infestation", "Rice Consumption", "Farm Area",
"Number of Farms", "Luzon", "Visayas", "Area per Farm"),
omit.stat = c("f", "ser"),
align = TRUE))##
## Regression Model Results
## =======================================================================================================================================
## Dependent variable:
## ------------------------------------------------------------------------------------------------------------
## Rainfed Palay Production
## (1) (2) (3) (4) (5) (6) (7)
## ---------------------------------------------------------------------------------------------------------------------------------------
## Drought 634.231*** 632.416*** 636.325*** 643.429*** 634.032*** 559.706*** 557.019***
## (108.321) (107.031) (103.844) (100.639) (100.374) (61.232) (61.237)
##
## Training Organizations 22.236 12.024
## (82.963) (67.302)
##
## Pesticide Dealers 332.748* 341.938** 340.508** 364.168** 458.278*** 440.027*** 454.577***
## (168.630) (161.636) (160.010) (140.848) (109.937) (108.073) (107.288)
##
## Wage Rate 91.194 76.917 86.693
## (288.439) (278.144) (270.299)
##
## HH Agricultural Activities -2.775*** -2.809*** -2.857*** -2.891*** -2.803*** -2.669*** -2.716***
## (0.550) (0.522) (0.441) (0.425) (0.418) (0.392) (0.389)
##
## HH Over 10 0.032** 0.032*** 0.033*** 0.033*** 0.034*** 0.035*** 0.034***
## (0.012) (0.012) (0.012) (0.012) (0.012) (0.011) (0.011)
##
## Infestation -126.306 -117.565 -113.986 -122.975 -82.652
## (112.241) (103.636) (100.783) (96.029) (88.386)
##
## Rice Consumption 1,142.009*** 1,130.005*** 1,118.449*** 1,108.004*** 1,041.822*** 1,032.424*** 1,013.747***
## (280.017) (271.919) (261.780) (257.657) (250.381) (249.908) (249.513)
##
## Farm Area 0.070
## (0.328)
##
## Number of Farms 0.141 0.219 0.242 0.225
## (0.446) (0.256) (0.219) (0.211)
##
## Luzon -35,738.330** -35,744.950*** -35,628.290*** -34,345.910*** -37,094.940*** -33,162.090*** -37,629.650***
## (13,444.220) (13,324.880) (13,191.230) (12,469.580) (12,215.230) (11,455.980) (10,652.770)
##
## Visayas -13,890.320 -14,925.850 -14,103.890 -16,755.680 -19,721.670 -21,533.790 -27,337.060**
## (18,733.760) (17,938.880) (17,185.650) (14,944.900) (14,701.830) (14,557.780) (13,488.780)
##
## Area per Farm 8,466.159 12,158.660 11,757.580 11,509.470 10,117.810 9,315.661
## (19,688.340) (9,432.798) (9,080.792) (8,975.016) (8,890.474) (8,839.414)
##
## Constant -169,279.700** -168,836.800** -168,885.900** -151,146.100*** -137,103.500*** -141,126.100*** -125,108.400***
## (67,746.930) (67,114.490) (66,520.650) (36,658.230) (34,255.120) (33,947.170) (30,383.120)
##
## ---------------------------------------------------------------------------------------------------------------------------------------
## Observations 67 67 67 67 67 67 67
## R2 0.858 0.858 0.857 0.857 0.854 0.852 0.849
## Adjusted R2 0.823 0.826 0.829 0.832 0.831 0.832 0.831
## =======================================================================================================================================
## Note: *p<0.1; **p<0.05; ***p<0.01rainfedTest <- subset(data1Test, select = -c(palay, irrigated))
test_rainfed <- subset(rainfedTest, select = -c(rainfed))
predictions_rainfed1 <- predict(rainfedmodel1, newdata = test_rainfed)
predictions_rainfed2 <- predict(rainfedmodel2, newdata = test_rainfed)
predictions_rainfed3 <- predict(rainfedmodel3, newdata = test_rainfed)
predictions_rainfed4 <- predict(rainfedmodel4, newdata = test_rainfed)
predictions_rainfed5 <- predict(rainfedmodel5, newdata = test_rainfed)
predictions_rainfed6 <- predict(rainfedmodel6, newdata = test_rainfed)
predictions_bestrainfedmodel <- predict(bestrainfedmodel, newdata = test_rainfed)
actual_rainfed <- rainfedTest$rainfed
rmserainfed1 <- sqrt(mean((predictions_rainfed1 - actual_rainfed)^2))
rmserainfed2 <- sqrt(mean((predictions_rainfed2 - actual_rainfed)^2))
rmserainfed3 <- sqrt(mean((predictions_rainfed3 - actual_rainfed)^2))
rmserainfed4 <- sqrt(mean((predictions_rainfed4 - actual_rainfed)^2))
rmserainfed5 <- sqrt(mean((predictions_rainfed5 - actual_rainfed)^2))
rmserainfed6 <- sqrt(mean((predictions_rainfed6 - actual_rainfed)^2))
rmserainfed7 <- sqrt(mean((predictions_bestrainfedmodel - actual_rainfed)^2))
rownames <- c('Model1', 'Model2', 'Model3', 'Model4', 'Model5', 'Model6', 'Model7')
rmserainfed <- data.frame(row.names = rownames, RMSE = c(rmserainfed1,rmserainfed2,rmserainfed3,rmserainfed4,rmserainfed5,rmserainfed6,rmserainfed7))
rmserainfedModel 5 has the lowest RMSE which means it is the best predictor.
# R-squared and Adjusted R-squared
rsquared5 <- summary(rainfedmodel5)$r.squared
adjrsquared5 <- summary(rainfedmodel5)$adj.r.squared
rsquared6 <- summary(rainfedmodel6)$r.squared
adjrsquared6 <- summary(rainfedmodel6)$adj.r.squared
rsquared7 <- summary(bestrainfedmodel)$r.squared
adjrsquared7 <- summary(bestrainfedmodel)$adj.r.squared
# AIC and BIC
aic5 <- AIC(rainfedmodel5)
bic5 <- BIC(rainfedmodel5)
aic6 <- AIC(rainfedmodel6)
bic6 <- BIC(rainfedmodel6)
aic7 <- AIC(bestrainfedmodel)
bic7 <- BIC(bestrainfedmodel)
# Mallows' Cp
cp5 <- ols_mallows_cp(rainfedmodel5, rainfedmodel1)
cp6 <- ols_mallows_cp(rainfedmodel6, rainfedmodel1)
cp7 <- ols_mallows_cp(bestrainfedmodel, rainfedmodel1)
model_comparison2 <- data.frame(
Model = c("Model 5", "Model 6", "Best Rainfed Model"),
R_squared = c(rsquared5, rsquared6, rsquared7),
Adjusted_R_squared = c(adjrsquared5, adjrsquared6, adjrsquared7),
AIC = c(aic5, aic6, aic7),
BIC = c(bic5, bic6, bic7),
Mallows_Cp = c(cp5, cp6, cp7),
RMSE = c(rmserainfed5, rmserainfed6, rmserainfed7)
)
print(model_comparison2)## Model R_squared Adjusted_R_squared AIC BIC Mallows_Cp
## 1 Model 5 0.8542579 0.8312460 1601.635 1625.887 7.257343
## 2 Model 6 0.8520220 0.8316113 1600.655 1622.702 6.089736
## 3 Best Rainfed Model 0.8491883 0.8312954 1599.926 1619.768 5.144665
## RMSE
## 1 65637.57
## 2 65995.40
## 3 66291.27While Model 5 has higher R-square and Adjusted R-square, Model 7 has the lowest AIC, BIC and Mallows’ Cp. Our final model for the Palay Production is given by:
bestrainfedmodel$coefficients## (Intercept) drought pest_dealers hh_agri_act
## -1.251084e+05 5.570192e+02 4.545766e+02 -2.716421e+00
## hh_over_10 rice_consumption Luzon Visayas
## 3.407467e-02 1.013747e+03 -3.762965e+04 -2.733705e+04# Calculate standardized coefficients
standardized_rainfedmodel <- lm.beta(bestrainfedmodel)
summary(standardized_rainfedmodel)##
## Call:
## lm(formula = rainfed ~ . - area_per_farm - farm_area - infestation -
## num_farms - wage_rate - train_org, data = rainfedTrain)
##
## Residuals:
## Min 1Q Median 3Q Max
## -67673 -22567 11 18183 94697
##
## Coefficients:
## Estimate Standardized Std. Error t value Pr(>|t|)
## (Intercept) -1.251e+05 NA 3.038e+04 -4.118 0.000121 ***
## drought 5.570e+02 7.700e-01 6.124e+01 9.096 7.90e-13 ***
## pest_dealers 4.546e+02 3.611e-01 1.073e+02 4.237 8.06e-05 ***
## hh_agri_act -2.716e+00 -4.450e-01 3.893e-01 -6.978 2.95e-09 ***
## hh_over_10 3.407e-02 1.610e-01 1.140e-02 2.990 0.004069 **
## rice_consumption 1.014e+03 2.212e-01 2.495e+02 4.063 0.000145 ***
## Luzon -3.763e+04 -2.252e-01 1.065e+04 -3.532 0.000807 ***
## Visayas -2.734e+04 -1.256e-01 1.349e+04 -2.027 0.047225 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 34550 on 59 degrees of freedom
## Multiple R-squared: 0.8492, Adjusted R-squared: 0.8313
## F-statistic: 47.46 on 7 and 59 DF, p-value: < 2.2e-16# Extract the standardized coefficients
standardized_rainfedcoefficients <- summary(standardized_rainfedmodel)$coefficients
# Create a data frame for plotting
coefficients_df <- as.data.frame(standardized_rainfedcoefficients)
coefficients_df$Predictor <- rownames(coefficients_df)
rownames(coefficients_df) <- NULL
# Extract only the Estimate column and rename for clarity
coefficients_df <- coefficients_df[, c("Predictor", "Standardized")]
names(coefficients_df)[2] <- "Estimate"
# Keep only the coefficients for predictors (ignore the intercept)
coefficients_df <- coefficients_df[-1, ]
# Plot the standardized coefficients
# Plot the standardized coefficients
ggplot(coefficients_df, aes(x = reorder(Predictor, Estimate), y = Estimate, fill=Predictor)) +
geom_bar(stat = "identity", show.legend = FALSE) +
coord_flip() + # Flip coordinates to make it horizontal
xlab("Predictor") +
ylab("Standardized Coefficient") +
ggtitle("Predictor Importance in Rainfed Palay Production Model") 
Diagnostic Checking
plot(bestrainfedmodel)
High-leverage Observation
For the residual plot of rainfed palays, observation 43 or the province of Iloilo is the only obvservation considered as a high-leverage observation in rainfed palay.
Multicollinearity
# Compute VIF for the model
vif_values <- vif(bestrainfedmodel)
print(vif_values)## drought pest_dealers hh_agri_act hh_over_10
## 2.803320 2.841432 1.591113 1.135358
## rice_consumption Luzon Visayas
## 1.159148 1.589514 1.501796Normality
shapiro.test(residuals(bestrainfedmodel))##
## Shapiro-Wilk normality test
##
## data: residuals(bestrainfedmodel)
## W = 0.98559, p-value = 0.6319The residuals are normally distributed.
Homoscedasticity
bptest(bestrainfedmodel)##
## studentized Breusch-Pagan test
##
## data: bestrainfedmodel
## BP = 15.946, df = 7, p-value = 0.02562There is some presence of heteroscedasticity.
Independence of Errors
dwtest(bestrainfedmodel)##
## Durbin-Watson test
##
## data: bestrainfedmodel
## DW = 2.3298, p-value = 0.854
## alternative hypothesis: true autocorrelation is greater than 0The residuals are independent.
Linearity
plot(bestrainfedmodel, which = 1)
abline(h = 0, col = "darkgreen", lty = 2) # Add horizontal line at y = 0
The random cloud indicates presence of linear relationship.
print(bestrainfedmodel)##
## Call:
## lm(formula = rainfed ~ . - area_per_farm - farm_area - infestation -
## num_farms - wage_rate - train_org, data = rainfedTrain)
##
## Coefficients:
## (Intercept) drought pest_dealers hh_agri_act
## -1.251e+05 5.570e+02 4.546e+02 -2.716e+00
## hh_over_10 rice_consumption Luzon Visayas
## 3.407e-02 1.014e+03 -3.763e+04 -2.734e+04Irrigated Palay Production
# Preparing the training set for Palay
irrigatedTrain <- subset(data1Train, select = -c(rainfed, palay))
str(irrigatedTrain)## 'data.frame': 67 obs. of 14 variables:
## $ drought : int 70 45 2 76 67 20 0 166 117 319 ...
## $ train_org : int 107 50 82 70 33 58 298 429 280 559 ...
## $ pest_dealers : int 29 39 55 33 26 27 119 93 71 220 ...
## $ wage_rate : num 179 179 179 179 179 ...
## $ hh_agri_act : int 4804 5715 5987 7219 8711 17028 7408 9783 6829 3362 ...
## $ hh_over_10 : int 109042 133703 103504 96817 84535 1144133 208417 228992 504579 1390728 ...
## $ infestation : int 104 74 44 107 103 62 0 172 120 262 ...
## $ rice_consumption: num 137 160 107 146 165 ...
## $ farm_area : num 17454 21370 17372 20825 27533 ...
## $ num_farms : num 46866 25505 52940 38237 37547 ...
## $ irrigated : num 55219 76320 14560 60660 163961 ...
## $ Luzon : num 1 1 1 1 1 1 1 1 1 1 ...
## $ Visayas : num 0 0 0 0 0 0 0 0 0 0 ...
## $ area_per_farm : num 0.372 0.838 0.328 0.545 0.733 ...# Creating the full model
irrigatedmodel1 <- lm(irrigated ~ ., data = irrigatedTrain)
summary(irrigatedmodel1)##
## Call:
## lm(formula = irrigated ~ ., data = irrigatedTrain)
##
## Residuals:
## Min 1Q Median 3Q Max
## -208279 -43913 -5612 39272 187871
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.382e+03 1.572e+05 0.028 0.9779
## drought 3.763e+02 2.513e+02 1.498 0.1402
## train_org -2.533e+02 1.925e+02 -1.316 0.1938
## pest_dealers 4.927e+03 3.912e+02 12.595 <2e-16 ***
## wage_rate -4.535e+02 6.692e+02 -0.678 0.5009
## hh_agri_act -6.233e-01 1.277e+00 -0.488 0.6275
## hh_over_10 1.717e-02 2.837e-02 0.605 0.5476
## infestation -6.478e+02 2.604e+02 -2.488 0.0160 *
## rice_consumption 7.912e+02 6.496e+02 1.218 0.2286
## farm_area 1.891e-01 7.604e-01 0.249 0.8046
## num_farms -1.258e+00 1.035e+00 -1.215 0.2298
## Luzon -5.461e+04 3.119e+04 -1.751 0.0857 .
## Visayas -1.026e+04 4.346e+04 -0.236 0.8144
## area_per_farm -5.224e+03 4.568e+04 -0.114 0.9094
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 82160 on 53 degrees of freedom
## Multiple R-squared: 0.919, Adjusted R-squared: 0.8992
## F-statistic: 46.28 on 13 and 53 DF, p-value: < 2.2e-16# Perform stepwise regression
stepwise.irrigated.model <- step(irrigatedmodel1, direction = "both")## Start: AIC=1528.7
## irrigated ~ drought + train_org + pest_dealers + wage_rate +
## hh_agri_act + hh_over_10 + infestation + rice_consumption +
## farm_area + num_farms + Luzon + Visayas + area_per_farm
##
## Df Sum of Sq RSS AIC
## - area_per_farm 1 8.8278e+07 3.5785e+11 1526.7
## - Visayas 1 3.7593e+08 3.5814e+11 1526.8
## - farm_area 1 4.1726e+08 3.5818e+11 1526.8
## - hh_agri_act 1 1.6083e+09 3.5937e+11 1527.0
## - hh_over_10 1 2.4731e+09 3.6024e+11 1527.2
## - wage_rate 1 3.1004e+09 3.6087e+11 1527.3
## - num_farms 1 9.9611e+09 3.6773e+11 1528.5
## - rice_consumption 1 1.0013e+10 3.6778e+11 1528.5
## <none> 3.5777e+11 1528.7
## - train_org 1 1.1693e+10 3.6946e+11 1528.8
## - drought 1 1.5139e+10 3.7290e+11 1529.5
## - Luzon 1 2.0694e+10 3.7846e+11 1530.5
## - infestation 1 4.1773e+10 3.9954e+11 1534.1
## - pest_dealers 1 1.0708e+12 1.4286e+12 1619.5
##
## Step: AIC=1526.71
## irrigated ~ drought + train_org + pest_dealers + wage_rate +
## hh_agri_act + hh_over_10 + infestation + rice_consumption +
## farm_area + num_farms + Luzon + Visayas
##
## Df Sum of Sq RSS AIC
## - Visayas 1 4.0931e+08 3.5826e+11 1524.8
## - farm_area 1 6.3717e+08 3.5849e+11 1524.8
## - hh_agri_act 1 1.7926e+09 3.5965e+11 1525.0
## - hh_over_10 1 2.6990e+09 3.6055e+11 1525.2
## - wage_rate 1 3.4166e+09 3.6127e+11 1525.3
## - rice_consumption 1 9.9252e+09 3.6778e+11 1526.5
## <none> 3.5785e+11 1526.7
## - train_org 1 1.5054e+10 3.7291e+11 1527.5
## - drought 1 1.5061e+10 3.7291e+11 1527.5
## - num_farms 1 1.7918e+10 3.7577e+11 1528.0
## - Luzon 1 2.0690e+10 3.7854e+11 1528.5
## + area_per_farm 1 8.8278e+07 3.5777e+11 1528.7
## - infestation 1 4.3504e+10 4.0136e+11 1532.4
## - pest_dealers 1 1.0951e+12 1.4530e+12 1618.6
##
## Step: AIC=1524.79
## irrigated ~ drought + train_org + pest_dealers + wage_rate +
## hh_agri_act + hh_over_10 + infestation + rice_consumption +
## farm_area + num_farms + Luzon
##
## Df Sum of Sq RSS AIC
## - farm_area 1 1.2834e+09 3.5955e+11 1523.0
## - hh_agri_act 1 1.6858e+09 3.5995e+11 1523.1
## - hh_over_10 1 2.4502e+09 3.6071e+11 1523.2
## - wage_rate 1 3.1706e+09 3.6143e+11 1523.4
## - rice_consumption 1 1.0463e+10 3.6873e+11 1524.7
## <none> 3.5826e+11 1524.8
## - drought 1 1.4784e+10 3.7305e+11 1525.5
## - train_org 1 1.5688e+10 3.7395e+11 1525.7
## - num_farms 1 1.7620e+10 3.7588e+11 1526.0
## - Luzon 1 2.0387e+10 3.7865e+11 1526.5
## + Visayas 1 4.0931e+08 3.5785e+11 1526.7
## + area_per_farm 1 1.2166e+08 3.5814e+11 1526.8
## - infestation 1 4.4097e+10 4.0236e+11 1530.6
## - pest_dealers 1 1.1957e+12 1.5540e+12 1621.1
##
## Step: AIC=1523.03
## irrigated ~ drought + train_org + pest_dealers + wage_rate +
## hh_agri_act + hh_over_10 + infestation + rice_consumption +
## num_farms + Luzon
##
## Df Sum of Sq RSS AIC
## - hh_over_10 1 2.5009e+09 3.6205e+11 1521.5
## - wage_rate 1 2.8783e+09 3.6242e+11 1521.6
## - hh_agri_act 1 2.8858e+09 3.6243e+11 1521.6
## - rice_consumption 1 9.3714e+09 3.6892e+11 1522.8
## <none> 3.5955e+11 1523.0
## - drought 1 1.4023e+10 3.7357e+11 1523.6
## - num_farms 1 2.0250e+10 3.7980e+11 1524.7
## + farm_area 1 1.2834e+09 3.5826e+11 1524.8
## + Visayas 1 1.0555e+09 3.5849e+11 1524.8
## + area_per_farm 1 7.1226e+08 3.5883e+11 1524.9
## - Luzon 1 2.2918e+10 3.8246e+11 1525.2
## - train_org 1 3.1827e+10 3.9137e+11 1526.7
## - infestation 1 4.3854e+10 4.0340e+11 1528.7
## - pest_dealers 1 1.2489e+12 1.6084e+12 1621.4
##
## Step: AIC=1521.49
## irrigated ~ drought + train_org + pest_dealers + wage_rate +
## hh_agri_act + infestation + rice_consumption + num_farms +
## Luzon
##
## Df Sum of Sq RSS AIC
## - hh_agri_act 1 2.1216e+09 3.6417e+11 1519.9
## - wage_rate 1 2.7240e+09 3.6477e+11 1520.0
## - rice_consumption 1 8.7934e+09 3.7084e+11 1521.1
## <none> 3.6205e+11 1521.5
## - drought 1 1.5968e+10 3.7801e+11 1522.4
## + hh_over_10 1 2.5009e+09 3.5955e+11 1523.0
## + farm_area 1 1.3341e+09 3.6071e+11 1523.2
## - num_farms 1 2.1031e+10 3.8308e+11 1523.3
## - Luzon 1 2.1497e+10 3.8354e+11 1523.4
## + Visayas 1 6.5620e+08 3.6139e+11 1523.4
## + area_per_farm 1 5.6266e+08 3.6148e+11 1523.4
## - train_org 1 2.9461e+10 3.9151e+11 1524.7
## - infestation 1 4.7130e+10 4.0918e+11 1527.7
## - pest_dealers 1 1.2761e+12 1.6381e+12 1620.6
##
## Step: AIC=1519.88
## irrigated ~ drought + train_org + pest_dealers + wage_rate +
## infestation + rice_consumption + num_farms + Luzon
##
## Df Sum of Sq RSS AIC
## - wage_rate 1 2.1461e+09 3.6631e+11 1518.3
## - rice_consumption 1 9.6968e+09 3.7387e+11 1519.7
## <none> 3.6417e+11 1519.9
## - drought 1 1.4285e+10 3.7845e+11 1520.5
## + farm_area 1 2.3372e+09 3.6183e+11 1521.5
## + hh_agri_act 1 2.1216e+09 3.6205e+11 1521.5
## + hh_over_10 1 1.7367e+09 3.6243e+11 1521.6
## + area_per_farm 1 1.0407e+09 3.6313e+11 1521.7
## + Visayas 1 8.3632e+08 3.6333e+11 1521.7
## - Luzon 1 2.2142e+10 3.8631e+11 1521.8
## - train_org 1 2.9484e+10 3.9365e+11 1523.1
## - num_farms 1 3.0302e+10 3.9447e+11 1523.2
## - infestation 1 4.5672e+10 4.0984e+11 1525.8
## - pest_dealers 1 1.2876e+12 1.6517e+12 1619.2
##
## Step: AIC=1518.28
## irrigated ~ drought + train_org + pest_dealers + infestation +
## rice_consumption + num_farms + Luzon
##
## Df Sum of Sq RSS AIC
## - rice_consumption 1 1.0299e+10 3.7661e+11 1518.1
## <none> 3.6631e+11 1518.3
## - drought 1 1.3529e+10 3.7984e+11 1518.7
## + wage_rate 1 2.1461e+09 3.6417e+11 1519.9
## + farm_area 1 1.8176e+09 3.6450e+11 1520.0
## + hh_over_10 1 1.7113e+09 3.6460e+11 1520.0
## + hh_agri_act 1 1.5438e+09 3.6477e+11 1520.0
## + area_per_farm 1 5.7076e+08 3.6574e+11 1520.2
## + Visayas 1 1.0168e+05 3.6631e+11 1520.3
## - num_farms 1 2.8244e+10 3.9456e+11 1521.2
## - train_org 1 3.2020e+10 3.9833e+11 1521.9
## - Luzon 1 3.6142e+10 4.0246e+11 1522.6
## - infestation 1 4.3918e+10 4.1023e+11 1523.9
## - pest_dealers 1 1.4891e+12 1.8554e+12 1625.0
##
## Step: AIC=1518.14
## irrigated ~ drought + train_org + pest_dealers + infestation +
## num_farms + Luzon
##
## Df Sum of Sq RSS AIC
## <none> 3.7661e+11 1518.1
## + rice_consumption 1 1.0299e+10 3.6631e+11 1518.3
## - drought 1 1.7190e+10 3.9380e+11 1519.1
## + wage_rate 1 2.7480e+09 3.7387e+11 1519.7
## + hh_agri_act 1 2.2717e+09 3.7434e+11 1519.7
## + hh_over_10 1 1.0979e+09 3.7552e+11 1519.9
## + farm_area 1 5.0317e+08 3.7611e+11 1520.0
## + area_per_farm 1 9.0669e+07 3.7652e+11 1520.1
## + Visayas 1 2.5148e+05 3.7661e+11 1520.1
## - Luzon 1 2.8358e+10 4.0497e+11 1521.0
## - num_farms 1 3.3631e+10 4.1024e+11 1521.9
## - infestation 1 3.7668e+10 4.1428e+11 1522.5
## - train_org 1 4.4034e+10 4.2065e+11 1523.5
## - pest_dealers 1 1.4859e+12 1.8626e+12 1623.2stepwise.irrigated.model##
## Call:
## lm(formula = irrigated ~ drought + train_org + pest_dealers +
## infestation + num_farms + Luzon, data = irrigatedTrain)
##
## Coefficients:
## (Intercept) drought train_org pest_dealers infestation
## 13979.466 335.039 -296.733 4775.120 -495.925
## num_farms Luzon
## -1.114 -51880.660# Create the models produced by stepwise regression
irrigatedmodel2 <- lm(irrigated ~ .-area_per_farm, data=irrigatedTrain)
irrigatedmodel3 <- lm(irrigated ~ .-area_per_farm -Visayas, data=irrigatedTrain)
irrigatedmodel4 <- lm(irrigated ~ .-area_per_farm -Visayas -farm_area, data=irrigatedTrain)
irrigatedmodel5 <- lm(irrigated ~ .-area_per_farm -Visayas -farm_area -hh_over_10, data=irrigatedTrain)
irrigatedmodel6 <- lm(irrigated ~ .-area_per_farm -Visayas -farm_area -hh_over_10 -hh_agri_act, data=irrigatedTrain)
irrigatedmodel7 <- lm(irrigated ~ .-area_per_farm -Visayas -farm_area -hh_over_10 -hh_agri_act -wage_rate, data=irrigatedTrain)
bestirrigatedmodel <- lm(irrigated ~ .-area_per_farm -Visayas -farm_area -hh_over_10 -hh_agri_act -wage_rate -rice_consumption, data=irrigatedTrain)# Summarize models in one table
irrigatedmodels <- list(irrigatedmodel1, irrigatedmodel2, irrigatedmodel3, irrigatedmodel4,
irrigatedmodel5, irrigatedmodel6, irrigatedmodel7, bestirrigatedmodel)
resultsirrigated <- data.frame(stargazer(irrigatedmodels,
type = "text",
title = "Regression Model Results",
dep.var.labels = "Irrigated Palay Production",
covariate.labels = c("Drought", "Training Organizations", "Pesticide Dealers",
"Wage Rate", "HH Agricultural Activities", "HH Over 10",
"Infestation", "Rice Consumption", "Farm Area",
"Number of Farms", "Luzon", "Visayas", "Area per Farm"),
omit.stat = c("f", "ser"),
align = TRUE))##
## Regression Model Results
## ==========================================================================================================================================
## Dependent variable:
## ---------------------------------------------------------------------------------------------------------------
## Irrigated Palay Production
## (1) (2) (3) (4) (5) (6) (7) (8)
## ------------------------------------------------------------------------------------------------------------------------------------------
## Drought 376.345 374.912 370.492 358.724 379.263 309.120 299.939 335.039
## (251.304) (248.687) (245.925) (242.733) (239.200) (204.937) (203.192) (202.458)
##
## Training Organizations -253.326 -262.324 -266.520 -310.405** -291.544** -252.287** -260.832** -296.733**
## (192.473) (174.045) (171.740) (139.417) (135.372) (116.423) (114.856) (112.032)
##
## Pesticide Dealers 4,927.363*** 4,933.648*** 4,904.251*** 4,932.740*** 4,889.501*** 4,860.993*** 4,780.700*** 4,775.120***
## (391.222) (383.786) (361.976) (353.677) (344.963) (339.452) (308.693) (310.353)
##
## Wage Rate -453.515 -467.758 -381.734 -362.584 -352.573 -310.008
## (669.177) (651.448) (547.156) (541.528) (538.384) (530.252)
##
## HH Agricultural Activities -0.623 -0.648 -0.627 -0.786 -0.664
## (1.277) (1.247) (1.233) (1.172) (1.149)
##
## HH Over 10 0.017 0.018 0.017 0.017
## (0.028) (0.028) (0.027) (0.027)
##
## Infestation -647.780** -640.392** -643.785** -609.052** -626.718*** -583.573*** -544.994** -495.925**
## (260.399) (249.942) (247.431) (233.043) (230.075) (216.374) (204.914) (202.443)
##
## Rice Consumption 791.226 784.050 800.622 739.514 714.851 747.435 768.852
## (649.638) (640.664) (631.726) (612.106) (607.550) (601.445) (596.970)
##
## Farm Area 0.189 0.113 0.148
## (0.760) (0.364) (0.333)
##
## Number of Farms -1.257 -1.172 -1.160 -0.961* -0.978* -1.092** -1.030** -1.114**
## (1.035) (0.713) (0.705) (0.541) (0.538) (0.497) (0.483) (0.481)
##
## Luzon -54,612.020* -54,154.240* -52,548.930* -54,856.280* -52,776.810* -53,511.020* -61,022.220** -51,880.660**
## (31,190.540) (30,648.620) (29,703.590) (29,034.770) (28,688.080) (28,494.920) (25,292.010) (24,408.510)
##
## Visayas -10,256.610 -10,665.990
## (43,462.270) (42,916.950)
##
## Area per Farm -5,223.514
## (45,676.880)
##
## Constant 4,381.904 2,175.607 -18,101.310 -11,916.970 -5,223.424 -18,129.860 -76,484.610 13,979.470
## (157,172.600) (154,552.100) (130,137.300) (128,458.400) (127,322.500) (124,627.100) (74,210.470) (24,078.970)
##
## ------------------------------------------------------------------------------------------------------------------------------------------
## Observations 67 67 67 67 67 67 67 67
## R2 0.919 0.919 0.919 0.919 0.918 0.918 0.917 0.915
## Adjusted R2 0.899 0.901 0.903 0.904 0.905 0.906 0.907 0.906
## ==========================================================================================================================================
## Note: *p<0.1; **p<0.05; ***p<0.01# Validate using RMSE
irrigatedTest <- subset(data1Test, select = -c(rainfed, palay))
test_irrigated <- subset(irrigatedTest, select = -c(irrigated))
predictions_irrigated1 <- predict(irrigatedmodel1, newdata = test_irrigated)
predictions_irrigated2 <- predict(irrigatedmodel2, newdata = test_irrigated)
predictions_irrigated3 <- predict(irrigatedmodel3, newdata = test_irrigated)
predictions_irrigated4 <- predict(irrigatedmodel4, newdata = test_irrigated)
predictions_irrigated5 <- predict(irrigatedmodel5, newdata = test_irrigated)
predictions_irrigated6 <- predict(irrigatedmodel6, newdata = test_irrigated)
predictions_irrigated7 <- predict(irrigatedmodel7, newdata = test_irrigated)
predictions_bestirrigatedmodel <- predict(bestirrigatedmodel, newdata= test_irrigated)
actual_irrigated <- irrigatedTest$irrigated
rmseirrigated1 <- sqrt(mean((predictions_irrigated1 - actual_irrigated)^2))
rmseirrigated2 <- sqrt(mean((predictions_irrigated2 - actual_irrigated)^2))
rmseirrigated3 <- sqrt(mean((predictions_irrigated3 - actual_irrigated)^2))
rmseirrigated4 <- sqrt(mean((predictions_irrigated4 - actual_irrigated)^2))
rmseirrigated5 <- sqrt(mean((predictions_irrigated5 - actual_irrigated)^2))
rmseirrigated6 <- sqrt(mean((predictions_irrigated6 - actual_irrigated)^2))
rmseirrigated7 <- sqrt(mean((predictions_irrigated7 - actual_irrigated)^2))
rmseirrigated8 <- sqrt(mean((predictions_bestirrigatedmodel - actual_irrigated)^2))
rownames <- c('Model1', 'Model2', 'Model3', 'Model4', 'Model5', 'Model6','Model7','Model8')
rmseirrigated <- data.frame(row.names = rownames, RMSE = c(rmseirrigated1, rmseirrigated2, rmseirrigated3, rmseirrigated4, rmseirrigated5,
rmseirrigated6, rmseirrigated7, rmseirrigated8))
rmseirrigated# R-squared and Adjusted R-squared
rsquared8 <- summary(bestirrigatedmodel)$r.squared
adjrsquared8 <- summary(bestirrigatedmodel)$adj.r.squared
# AIC and BIC
aic8 <- AIC(bestrainfedmodel)
bic8 <- BIC(bestrainfedmodel)
# Mallows' Cp
cp8 <- ols_mallows_cp(bestirrigatedmodel, rainfedmodel1)
model_comparison3 <- data.frame(
Model = "Best Irrigated Model",
R_squared = rsquared8,
Adjusted_R_squared = adjrsquared8,
AIC = aic8,
BIC = bic8,
Mallows_Cp = cp8,
RMSE = rmseirrigated8)
print(model_comparison3)## Model R_squared Adjusted_R_squared AIC BIC
## 1 Best Irrigated Model 0.914781 0.9062591 1599.926 1619.768
## Mallows_Cp RMSE
## 1 247.2952 95400.42# Recall that Model 8 is the best model from the stepwise regression.
bestirrigatedmodel$coefficients## (Intercept) drought train_org pest_dealers infestation
## 13979.465732 335.038805 -296.732969 4775.119653 -495.924811
## num_farms Luzon
## -1.113845 -51880.660349# Calculate standardized coefficients
standardized_irrigatedmodel <- lm.beta(bestirrigatedmodel)
summary(standardized_irrigatedmodel)##
## Call:
## lm(formula = irrigated ~ . - area_per_farm - Visayas - farm_area -
## hh_over_10 - hh_agri_act - wage_rate - rice_consumption,
## data = irrigatedTrain)
##
## Residuals:
## Min 1Q Median 3Q Max
## -207646 -45867 -7142 44686 202900
##
## Coefficients:
## Estimate Standardized Std. Error t value Pr(>|t|)
## (Intercept) 1.398e+04 NA 2.408e+04 0.581 0.5637
## drought 3.350e+02 1.505e-01 2.025e+02 1.655 0.1032
## train_org -2.967e+02 -1.726e-01 1.120e+02 -2.649 0.0103 *
## pest_dealers 4.775e+03 1.233e+00 3.104e+02 15.386 <2e-16 ***
## infestation -4.959e+02 -2.337e-01 2.024e+02 -2.450 0.0172 *
## num_farms -1.114e+00 -2.023e-01 4.812e-01 -2.315 0.0241 *
## Luzon -5.188e+04 -1.009e-01 2.441e+04 -2.126 0.0377 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 79230 on 60 degrees of freedom
## Multiple R-squared: 0.9148, Adjusted R-squared: 0.9063
## F-statistic: 107.3 on 6 and 60 DF, p-value: < 2.2e-16# Extract the standardized coefficients
standardized_irrigatedcoefficients <- summary(standardized_irrigatedmodel)$coefficients
# Create a data frame for plotting
coefficients_df <- as.data.frame(standardized_irrigatedcoefficients)
coefficients_df$Predictor <- rownames(coefficients_df)
rownames(coefficients_df) <- NULL
# Extract only the Estimate column and rename for clarity
coefficients_df <- coefficients_df[, c("Predictor", "Standardized")]
names(coefficients_df)[2] <- "Estimate"
# Keep only the coefficients for predictors (ignore the intercept)
coefficients_df <- coefficients_df[-1, ]
# Plot the standardized coefficients
# Plot the standardized coefficients
ggplot(coefficients_df, aes(x = reorder(Predictor, Estimate), y = Estimate, fill=Predictor)) +
geom_bar(stat = "identity", show.legend = FALSE) +
coord_flip() + # Flip coordinates to make it horizontal
xlab("Predictor") +
ylab("Standardized Coefficient") +
ggtitle("Predictor Importance in Irrigated Palay Production Model") 
Diagnostic Checking
# Diagnostic Checking
plot(bestirrigatedmodel)
Multicollinearity
# Compute VIF for the model
vif_values <- vif(bestirrigatedmodel)
# Print VIF values
print(vif_values)## drought train_org pest_dealers infestation num_farms Luzon
## 5.826135 2.989273 4.520651 6.406371 5.375617 1.586651Normality
# Perform Shapiro-Wilk Test
shapiro.test(residuals(bestirrigatedmodel))##
## Shapiro-Wilk normality test
##
## data: residuals(bestirrigatedmodel)
## W = 0.98207, p-value = 0.4454The residuals are normal.
Homoscedasticity
# Perform Breusch-Pagan Test
bptest(bestirrigatedmodel)##
## studentized Breusch-Pagan test
##
## data: bestirrigatedmodel
## BP = 28.87, df = 6, p-value = 6.438e-05There is presence of heteroscedasticity.
Independence of Error Terms
# Perform Durbin-Watson Test
dwtest(bestirrigatedmodel)##
## Durbin-Watson test
##
## data: bestirrigatedmodel
## DW = 2.1603, p-value = 0.626
## alternative hypothesis: true autocorrelation is greater than 0There is no serial autocorrelation.
Linearity
plot(bestirrigatedmodel, which = 1)
abline(h=0, col="red", lty=2) # Add horizontal line at y = 0
summary(bestirrigatedmodel)##
## Call:
## lm(formula = irrigated ~ . - area_per_farm - Visayas - farm_area -
## hh_over_10 - hh_agri_act - wage_rate - rice_consumption,
## data = irrigatedTrain)
##
## Residuals:
## Min 1Q Median 3Q Max
## -207646 -45867 -7142 44686 202900
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.398e+04 2.408e+04 0.581 0.5637
## drought 3.350e+02 2.025e+02 1.655 0.1032
## train_org -2.967e+02 1.120e+02 -2.649 0.0103 *
## pest_dealers 4.775e+03 3.104e+02 15.386 <2e-16 ***
## infestation -4.959e+02 2.024e+02 -2.450 0.0172 *
## num_farms -1.114e+00 4.812e-01 -2.315 0.0241 *
## Luzon -5.188e+04 2.441e+04 -2.126 0.0377 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 79230 on 60 degrees of freedom
## Multiple R-squared: 0.9148, Adjusted R-squared: 0.9063
## F-statistic: 107.3 on 6 and 60 DF, p-value: < 2.2e-16# Summarize three models in one table
allpalaymodels <- list(bestpalaymodel, bestrainfedmodel, bestirrigatedmodel)
resultsallpalaymodels <- stargazer(allpalaymodels,
type = "text",
title = "Regression Model Results",
dep.var.labels = c("Total Palay Production", "Rainfed Palay Production", "Irrigated Palay Production"),
covariate.labels = c("Drought", "Training Organizations", "Pesticide Dealers",
"Wage Rate", "HH Agricultural Activities", "HH Over 10",
"Infestation", "Rice Consumption", "Farm Area",
"Number of Farms", "Luzon", "Visayas", "Area per Farm"),
omit.stat = c("f", "ser"),
align = TRUE)##
## Regression Model Results
## =====================================================================================================
## Dependent variable:
## --------------------------------------------------------------------------
## Total Palay Production Rainfed Palay Production Irrigated Palay Production
## (1) (2) (3)
## -----------------------------------------------------------------------------------------------------
## Drought 1,020.292*** 557.019*** 335.039
## (248.808) (61.237) (202.458)
##
## Training Organizations -420.731*** -296.733**
## (131.618) (112.032)
##
## Pesticide Dealers 5,029.533*** 454.577*** 4,775.120***
## (281.486) (107.288) (310.353)
##
## Wage Rate -4.238*** -2.716***
## (1.124) (0.389)
##
## HH Agricultural Activities 0.049* 0.034***
## (0.028) (0.011)
##
## HH Over 10 -745.724*** -495.925**
## (224.752) (202.443)
##
## Infestation 1,889.434*** 1,013.747***
## (641.939) (249.513)
##
## Rice Consumption -1.114**
## (0.481)
##
## Farm Area -90,830.900*** -37,629.650*** -51,880.660**
## (27,071.690) (10,652.770) (24,408.510)
##
## Number of Farms -27,337.060**
## (13,488.780)
##
## Luzon -230,097.100*** -125,108.400*** 13,979.470
## (80,088.770) (30,383.120) (24,078.970)
##
## -----------------------------------------------------------------------------------------------------
## Observations 67 67 67
## R2 0.933 0.849 0.915
## Adjusted R2 0.924 0.831 0.906
## =====================================================================================================
## Note: *p<0.1; **p<0.05; ***p<0.01# Combine all models
model_comparison_final <- rbind(model_comparison[model_comparison$Model == "Best Total Palay Model", ], model_comparison2[model_comparison2$Model == "Best Rainfed Model", ], model_comparison3[model_comparison3$Model == "Best Irrigated Model", ])
# Reset row names
rownames(model_comparison_final) <- NULL
model_comparison_final