R Notebook
Stat 129 Final Project
Comparative Analysis of Rainfed and Irrigated Palay Volume of Production Using Multiple Linear Regression
I. Load libraries
# Load the libraries
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, unionlibrary(ggplot2)
library(corrplot)## corrplot 0.92 loadedlibrary(caret)## Loading required package: latticelibrary(dummy)## dummy 0.1.3## dummyNews()library(leaps)
library(car)## Loading required package: carData##
## Attaching package: 'car'## The following object is masked from 'package:dplyr':
##
## recodelibrary(olsrr)##
## Attaching package: 'olsrr'## The following object is masked from 'package:datasets':
##
## riverslibrary(stargazer)##
## Please cite as:## Hlavac, Marek (2022). stargazer: Well-Formatted Regression and Summary Statistics Tables.## R package version 5.2.3. https://CRAN.R-project.org/package=stargazerlibrary(lmtest)## Loading required package: zoo##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numericII. 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# Correlation Analysis
cor_matrix <- cor(data %>% select(- c(X, island_group)))
corrplot(cor_matrix, method = "circle")
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", "Model 6"),
R_squared = c(rsquared5, rsquared6),
Adjusted_R_squared = c(adjrsquared5, adjrsquared6),
AIC = c(aic5, aic6),
BIC = c(bic5, bic6),
Mallows_Cp = c(cp5, cp6)
)
print(model_comparison)## Model R_squared Adjusted_R_squared AIC BIC Mallows_Cp
## 1 Model 5 0.9343461 0.9239797 1721.470 1745.722 7.478926
## 2 Model 6 0.9330119 0.9237722 1720.818 1742.865 6.586070While 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+04Let’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