data <- read_excel("Data Pengantar Ekonometrika Kel 10 (1).xlsx")
data <- data %>%
rename(
KabKota = `Kab/Kota`,
ketersediaan = `Ketersediaan Beras ton`,
produktivitas = `Produktivitas Tanaman Padi (Ku/Ha)`,
luas_panen = `Luas Panen (Ha)`
)
clean_numeric <- function(x){
x <- gsub(",", ".", x)
x <- gsub(" ", "", x)
x <- gsub("\\s+", "", x)
x <- gsub("[^0-9\\.]", "", x)
as.numeric(x)
}
data$ketersediaan <- clean_numeric(as.character(data$ketersediaan))
data$produktivitas <- clean_numeric(as.character(data$produktivitas))
data$luas_panen <- clean_numeric(as.character(data$luas_panen))
data$Tahun <- as.numeric(data$Tahun)
data$KabKota <- as.factor(data$KabKota)
str(data)
## tibble [190 × 5] (S3: tbl_df/tbl/data.frame)
## $ Tahun : num [1:190] 2020 2020 2020 2020 2020 2020 2020 2020 2020 2020 ...
## $ KabKota : Factor w/ 38 levels "Kabupaten Bangkalan",..: 18 21 27 29 3 9 14 11 7 2 ...
## $ ketersediaan : num [1:190] 48469 217880 62619 123798 113664 ...
## $ produktivitas: num [1:190] 43.5 58 55.2 58.4 56.6 ...
## $ luas_panen : num [1:190] 83941 377333 108446 214398 196848 ...
summary(data)
## Tahun KabKota ketersediaan produktivitas
## Min. :2020 Kabupaten Bangkalan : 5 Min. : 1901 Min. :43.51
## 1st Qu.:2021 Kabupaten Banyuwangi: 5 1st Qu.: 51214 1st Qu.:52.74
## Median :2022 Kabupaten Blitar : 5 Median :124233 Median :57.42
## Mean :2022 Kabupaten Bojonegoro: 5 Mean :146609 Mean :57.23
## 3rd Qu.:2023 Kabupaten Bondowoso : 5 3rd Qu.:214922 3rd Qu.:61.04
## Max. :2024 Kabupaten Gresik : 5 Max. :521920 Max. :72.72
## (Other) :160
## luas_panen
## Min. : 503.2
## 1st Qu.: 18002.4
## Median : 43661.6
## Mean : 87896.3
## 3rd Qu.: 84187.7
## Max. :886061.0
##
pdata <- pdata.frame(data, index = c("KabKota", "Tahun"))
library(car)
## Loading required package: carData
##
## Attaching package: 'car'
## The following object is masked from 'package:dplyr':
##
## recode
# Model OLS untuk cek VIF (karena plm tidak bisa langsung dihitung VIF)
model_ols <- lm(ketersediaan ~ produktivitas + luas_panen, data = data)
# Hitung nilai VIF
vif(model_ols)
## produktivitas luas_panen
## 1.000735 1.000735
model_pool <- plm(ketersediaan ~ produktivitas + luas_panen,
data = pdata, model = "pooling")
model_fe <- plm(ketersediaan ~ produktivitas + luas_panen,
data = pdata, model = "within")
model_re <- plm(ketersediaan ~ produktivitas + luas_panen,
data = pdata, model = "random")
summary(model_pool)
## Pooling Model
##
## Call:
## plm(formula = ketersediaan ~ produktivitas + luas_panen, data = pdata,
## model = "pooling")
##
## Balanced Panel: n = 38, T = 5, N = 190
##
## Residuals:
## Min. 1st Qu. Median 3rd Qu. Max.
## -102171 -92994 -17678 44014 343398
##
## Coefficients:
## Estimate Std. Error t-value Pr(>|t|)
## (Intercept) 1.2041e+05 7.9630e+04 1.5121 0.1322
## produktivitas -3.4234e+02 1.3801e+03 -0.2480 0.8044
## luas_panen 5.2099e-01 5.5078e-02 9.4591 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Total Sum of Squares: 3.0267e+12
## Residual Sum of Squares: 2.0453e+12
## R-Squared: 0.32425
## Adj. R-Squared: 0.31702
## F-statistic: 44.8645 on 2 and 187 DF, p-value: < 2.22e-16
summary(model_fe)
## Oneway (individual) effect Within Model
##
## Call:
## plm(formula = ketersediaan ~ produktivitas + luas_panen, data = pdata,
## model = "within")
##
## Balanced Panel: n = 38, T = 5, N = 190
##
## Residuals:
## Min. 1st Qu. Median 3rd Qu. Max.
## -34355.09 -4129.36 456.38 4207.10 44306.96
##
## Coefficients:
## Estimate Std. Error t-value Pr(>|t|)
## produktivitas 1.4478e+03 2.9632e+02 4.8859 2.615e-06 ***
## luas_panen 3.0631e-02 5.9780e-03 5.1241 9.071e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Total Sum of Squares: 1.8182e+10
## Residual Sum of Squares: 1.3537e+10
## R-Squared: 0.25548
## Adj. R-Squared: 0.061903
## F-statistic: 25.7358 on 2 and 150 DF, p-value: 2.4592e-10
summary(model_re)
## Oneway (individual) effect Random Effect Model
## (Swamy-Arora's transformation)
##
## Call:
## plm(formula = ketersediaan ~ produktivitas + luas_panen, data = pdata,
## model = "random")
##
## Balanced Panel: n = 38, T = 5, N = 190
##
## Effects:
## var std.dev share
## idiosyncratic 9.025e+07 9.500e+03 0.993
## individual 6.005e+05 7.749e+02 0.007
## theta: 0.01623
##
## Residuals:
## Min. 1st Qu. Median 3rd Qu. Max.
## -101476 -92214 -17862 43356 340345
##
## Coefficients:
## Estimate Std. Error z-value Pr(>|z|)
## (Intercept) 1.2142e+05 7.9779e+04 1.5220 0.1280
## produktivitas -3.4294e+02 1.3828e+03 -0.2480 0.8041
## luas_panen 5.0984e-01 5.4715e-02 9.3181 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Total Sum of Squares: 2.9298e+12
## Residual Sum of Squares: 1.999e+12
## R-Squared: 0.3177
## Adj. R-Squared: 0.3104
## Chisq: 87.0721 on 2 DF, p-value: < 2.22e-16
uji_chow <- pFtest(model_fe, model_pool)
uji_chow
##
## F test for individual effects
##
## data: ketersediaan ~ produktivitas + luas_panen
## F = 608.48, df1 = 37, df2 = 150, p-value < 2.2e-16
## alternative hypothesis: significant effects
uji_hausman <- phtest(model_fe, model_re)
uji_hausman
##
## Hausman Test
##
## data: ketersediaan ~ produktivitas + luas_panen
## chisq = 80.105, df = 2, p-value < 2.2e-16
## alternative hypothesis: one model is inconsistent
uji_lm <- plmtest(model_pool, type = "bp")
uji_lm
##
## Lagrange Multiplier Test - (Breusch-Pagan)
##
## data: ketersediaan ~ produktivitas + luas_panen
## chisq = 145.97, df = 1, p-value < 2.2e-16
## alternative hypothesis: significant effects
plmtest(model_re, effect = "individual")
##
## Lagrange Multiplier Test - (Honda)
##
## data: ketersediaan ~ produktivitas + luas_panen
## normal = 12.082, p-value < 2.2e-16
## alternative hypothesis: significant effects
plmtest(model_re, effect = "time")
##
## Lagrange Multiplier Test - time effects (Honda)
##
## data: ketersediaan ~ produktivitas + luas_panen
## normal = 8.877, p-value < 2.2e-16
## alternative hypothesis: significant effects
plmtest(model_re, effect = "twoways")
##
## Lagrange Multiplier Test - two-ways effects (Honda)
##
## data: ketersediaan ~ produktivitas + luas_panen
## normal = 14.82, p-value < 2.2e-16
## alternative hypothesis: significant effects
#FE
pFtest(model_fe, model_pool)
##
## F test for individual effects
##
## data: ketersediaan ~ produktivitas + luas_panen
## F = 608.48, df1 = 37, df2 = 150, p-value < 2.2e-16
## alternative hypothesis: significant effects
# RE (LM test)
plmtest(model_pool, type = "bp")
##
## Lagrange Multiplier Test - (Breusch-Pagan)
##
## data: ketersediaan ~ produktivitas + luas_panen
## chisq = 145.97, df = 1, p-value < 2.2e-16
## alternative hypothesis: significant effects
summary(model_re)
## Oneway (individual) effect Random Effect Model
## (Swamy-Arora's transformation)
##
## Call:
## plm(formula = ketersediaan ~ produktivitas + luas_panen, data = pdata,
## model = "random")
##
## Balanced Panel: n = 38, T = 5, N = 190
##
## Effects:
## var std.dev share
## idiosyncratic 9.025e+07 9.500e+03 0.993
## individual 6.005e+05 7.749e+02 0.007
## theta: 0.01623
##
## Residuals:
## Min. 1st Qu. Median 3rd Qu. Max.
## -101476 -92214 -17862 43356 340345
##
## Coefficients:
## Estimate Std. Error z-value Pr(>|z|)
## (Intercept) 1.2142e+05 7.9779e+04 1.5220 0.1280
## produktivitas -3.4294e+02 1.3828e+03 -0.2480 0.8041
## luas_panen 5.0984e-01 5.4715e-02 9.3181 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Total Sum of Squares: 2.9298e+12
## Residual Sum of Squares: 1.999e+12
## R-Squared: 0.3177
## Adj. R-Squared: 0.3104
## Chisq: 87.0721 on 2 DF, p-value: < 2.22e-16
summary(model_re)$r.squared
## rsq adjrsq
## 0.3176979 0.3104005
coeftest(model_re, vcov = vcovHC(model_re, type = "HC1"))
##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.2142e+05 1.1847e+05 1.0249 0.3067
## produktivitas -3.4294e+02 2.1064e+03 -0.1628 0.8708
## luas_panen 5.0984e-01 4.4982e-02 11.3345 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Hitung rata-rata per tahun seluruh Jatim
jatim_mean <- data %>%
group_by(Tahun) %>%
summarise(rata_ketersediaan = mean(ketersediaan))
# Plot line chart rata-rata
ggplot(jatim_mean, aes(x = Tahun, y = rata_ketersediaan)) +
geom_line(size = 1, color = "#2C75FF") +
geom_point(size = 3, color = "#2C75FF") +
labs(
title = "Rata-rata Ketersediaan Beras Jawa Timur per Tahun",
subtitle = "Data Gabungan Seluruh Kabupaten/Kota",
x = "Tahun",
y = "Rata-rata (Ton)"
) +
theme_minimal(base_size = 14)
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
jatim_total <- data %>%
group_by(Tahun) %>%
summarise(total_ketersediaan = sum(ketersediaan))
ggplot(jatim_total, aes(x = factor(Tahun), y = total_ketersediaan)) +
geom_line(aes(group = 1), size = 1.3, color = "#2CA02C") +
geom_point(size = 4, color = "#2CA02C") +
labs(
title = "Total Ketersediaan Beras Jawa Timur per Tahun",
subtitle = "Akumulasi Seluruh Kabupaten/Kota",
x = "Tahun",
y = "Total Ketersediaan (Ton)"
) +
theme_minimal(base_size = 14)
ggplot(jatim_total, aes(x = factor(Tahun), y = total_ketersediaan)) +
geom_bar(stat = "identity", fill = "#F0D34D") +
labs(
title = "Jumlah Ketersediaan Beras Jawa Timur Tahun 2020–2024",
subtitle = "Akumulasi Seluruh Kabupaten/Kota",
x = "Tahun",
y = "Ketersediaan (Ton)"
) +
theme_minimal(base_size = 14)
