Input Data

Berikut merupakan data Ketahanan Pangan provinsi Sulawesi Selatan Tahun 2020-2023 dengan variabel rasio ketersediaan beras (Y),konsumsi beras(Ton)(X1), Produktivitas Padi (GKG)(Kw/Ha)(X2),luas lahan panen (ha)(X3), dan jumlah penduduk (ratusan ribu).

library(readxl)
data <- read_excel("D:/Semester 5/Analisis Data Panel/Data Projek/data_tambah.xlsx")
head(data)

Eksplorasi Data

library(psych)
## Warning: package 'psych' was built under R version 4.3.3
describeBy(data$Y,group = data$Tahun) #Rasio Ketersediaan Beras
## 
##  Descriptive statistics by group 
## group: 2020
##    vars  n mean   sd median trimmed  mad  min max range skew kurtosis   se
## X1    1 24 2.91 2.37   2.04    2.72 1.82 0.04 7.8  7.75 0.78    -0.68 0.48
## ------------------------------------------------------------ 
## group: 2021
##    vars  n mean   sd median trimmed mad  min  max range skew kurtosis   se
## X1    1 24 3.11 2.53   2.28    2.87   2 0.04 8.67  8.63 0.87     -0.4 0.52
## ------------------------------------------------------------ 
## group: 2022
##    vars  n mean   sd median trimmed  mad  min   max range skew kurtosis   se
## X1    1 24 3.26 2.84   2.23    2.92 2.04 0.04 10.33 10.28 1.05     0.04 0.58
## ------------------------------------------------------------ 
## group: 2023
##    vars  n mean   sd median trimmed  mad  min  max range skew kurtosis   se
## X1    1 24 2.78 2.37   2.13     2.5 1.75 0.05 8.77  8.72 0.97    -0.01 0.48
describeBy(data$X1,group = data$Tahun) #Konsumsi Beras
## 
##  Descriptive statistics by group 
## group: 2020
##    vars  n     mean       sd   median  trimmed      mad      min      max
## X1    1 24 40198.49 32052.03 32550.57 34182.19 12056.92 14676.24 166953.1
##       range skew kurtosis      se
## X1 152276.9 2.68     7.52 6542.59
## ------------------------------------------------------------ 
## group: 2021
##    vars  n     mean       sd   median  trimmed      mad      min      max
## X1    1 24 40727.11 32574.69 33094.86 34582.12 12264.86 14854.93 169741.7
##       range skew kurtosis      se
## X1 154886.8 2.69     7.57 6649.28
## ------------------------------------------------------------ 
## group: 2022
##    vars  n     mean      sd   median  trimmed      mad      min      max
## X1    1 24 40870.31 32787.5 33323.71 34655.98 12354.97 14893.35 170909.2
##       range skew kurtosis      se
## X1 156015.8  2.7     7.63 6692.72
## ------------------------------------------------------------ 
## group: 2023
##    vars  n     mean       sd   median  trimmed      mad      min      max
## X1    1 24 42705.26 30665.92 35288.86 37197.05 14445.99 15442.72 158080.2
##       range skew kurtosis      se
## X1 142637.5 2.38     5.86 6259.66
describeBy(data$X2,group = data$Tahun) #Produktivitas Padi
## 
##  Descriptive statistics by group 
## group: 2020
##    vars  n  mean   sd median trimmed  mad   min max range skew kurtosis   se
## X1    1 24 48.04 5.77     47   47.75 4.54 38.61  61 22.39 0.61    -0.44 1.18
## ------------------------------------------------------------ 
## group: 2021
##    vars  n  mean   sd median trimmed  mad   min   max range skew kurtosis   se
## X1    1 24 51.44 6.99   50.7   51.21 6.02 36.84 66.46 29.62 0.26    -0.41 1.43
## ------------------------------------------------------------ 
## group: 2022
##    vars  n  mean   sd median trimmed  mad   min   max range skew kurtosis  se
## X1    1 24 50.16 5.88  48.92   49.84 5.23 41.39 61.88 20.49 0.54    -0.83 1.2
## ------------------------------------------------------------ 
## group: 2023
##    vars  n  mean   sd median trimmed  mad   min   max range skew kurtosis   se
## X1    1 24 49.09 5.15  48.55   49.06 5.62 40.08 58.38  18.3 0.12    -1.01 1.05
describeBy(data$X3,group = data$Tahun) #Luas Panen
## 
##  Descriptive statistics by group 
## group: 2020
##    vars  n     mean    sd median trimmed      mad min    max  range skew
## X1    1 24 40677.46 41336  26937 33962.4 25686.04 998 164096 163098 1.48
##    kurtosis      se
## X1     1.59 8437.68
## ------------------------------------------------------------ 
## group: 2021
##    vars  n     mean       sd  median trimmed      mad min    max  range skew
## X1    1 24 41048.33 41505.99 28283.5 34193.6 28103.42 968 165260 164292 1.52
##    kurtosis      se
## X1     1.71 8472.37
## ------------------------------------------------------------ 
## group: 2022
##    vars  n     mean       sd   median  trimmed      mad    min      max
## X1    1 24 43253.51 46598.92 27266.89 34548.32 26453.59 963.45 186094.8
##       range skew kurtosis      se
## X1 185131.3 1.74     2.45 9511.96
## ------------------------------------------------------------ 
## group: 2023
##    vars  n     mean       sd   median  trimmed      mad    min      max
## X1    1 24 40324.59 42346.77 27062.44 32715.26 27849.45 962.16 170329.8
##       range skew kurtosis   se
## X1 169367.6 1.67     2.26 8644
describeBy(data$X4,group = data$Tahun) #Jumlah Penduduk
## 
##  Descriptive statistics by group 
## group: 2020
##    vars  n   mean    sd median trimmed   mad    min     max   range skew
## X1    1 24 378.06 275.6 310.42  327.97 123.5 137.07 1423.88 1286.81 2.43
##    kurtosis    se
## X1     6.12 56.26
## ------------------------------------------------------------ 
## group: 2021
##    vars  n   mean     sd median trimmed    mad    min     max   range skew
## X1    1 24 381.54 277.84 313.16   331.1 124.98 138.16 1434.44 1296.28 2.42
##    kurtosis    se
## X1     6.08 56.71
## ------------------------------------------------------------ 
## group: 2022
##    vars  n   mean    sd median trimmed    mad    min     max   range skew
## X1    1 24 385.84 280.1 316.79  335.18 126.78 139.66 1444.88 1305.22 2.41
##    kurtosis    se
## X1     6.01 57.18
## ------------------------------------------------------------ 
## group: 2023
##    vars  n  mean     sd median trimmed   mad    min     max   range skew
## X1    1 24 390.1 282.31 320.36  339.23 128.5 141.18 1454.96 1313.78  2.4
##    kurtosis    se
## X1     5.94 57.63

Plot Eksplorasi

library(plotly)
## Warning: package 'plotly' was built under R version 4.3.3
## Loading required package: ggplot2
## 
## Attaching package: 'ggplot2'
## The following objects are masked from 'package:psych':
## 
##     %+%, alpha
## 
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
## 
##     last_plot
## The following object is masked from 'package:stats':
## 
##     filter
## The following object is masked from 'package:graphics':
## 
##     layout
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
#Rasio Ketersediaan Beras
plot_ly(data, 
        x = ~factor(Tahun), 
        y = ~Y, 
        text = ~paste('Kabupaten/Kota: ', `Kab/Kota`, '<br>Value: ', Y),  # Menampilkan informasi saat hover
        hoverinfo = 'text',  # Tampilkan info saat hover
        type = "scatter",
        mode = "marker",
        color = ~factor(`Kab/Kota`))%>%
  layout(
    title = "Ketersediaan Beras",  # Judul grafik
    xaxis = list(title = "Tahun"),  # Label sumbu X
    yaxis = list(title = "Ketersediaan Beras"),  # Label sumbu Y
    legend = list(title = list(text = 'Kab/Kota')))  # Judul legend
## Warning in RColorBrewer::brewer.pal(N, "Set2"): n too large, allowed maximum for palette Set2 is 8
## Returning the palette you asked for with that many colors
## Warning in RColorBrewer::brewer.pal(N, "Set2"): n too large, allowed maximum for palette Set2 is 8
## Returning the palette you asked for with that many colors
#Konsumsi Beras
plot_ly(data, 
        x = ~factor(Tahun), 
        y = ~X1, 
        type = 'scatter', 
        text = ~paste('Kabupaten/Kota: ', `Kab/Kota`, '<br>Value: ', X1),  # Menampilkan informasi saat hover
        hoverinfo = 'text',  # Tampilkan info saat hover
        mode = "marker",
        color = ~factor(`Kab/Kota`))%>%  # Gunakan 'Kab/Kota' sebagai nama untuk legend
  layout(title = "Konsumsi Beras",
         xaxis = list(title = "Tahun"),  # Setel ticks x-axis
         yaxis = list(title = "Konsumsi Beras (Ton)"),
         legend = list(title = list(text = 'Kabupaten/Kota')))  # Judul
## Warning in RColorBrewer::brewer.pal(N, "Set2"): n too large, allowed maximum for palette Set2 is 8
## Returning the palette you asked for with that many colors
## Warning in RColorBrewer::brewer.pal(N, "Set2"): n too large, allowed maximum for palette Set2 is 8
## Returning the palette you asked for with that many colors
#Produktivitas Padi
plot_ly(data, 
        x = ~factor(Tahun), 
        y = ~X2, 
        type = 'scatter', 
        text = ~paste('Kabupaten/Kota: ', `Kab/Kota`, '<br>Value: ', X2),  # Menampilkan informasi saat hover
        hoverinfo = 'text',  # Tampilkan info saat hover
        mode = "marker",
        color = ~factor(`Kab/Kota`))%>%   # Gunakan 'Kab/Kota' sebagai nama untuk legend
  layout(title = "Produktivitas Padi (GKG)(Kw/Ha)",
         xaxis = list(title = "Tahun"),  # Setel ticks x-axis
         yaxis = list(title = "Produktivitas Padi"),
         legend = list(title = list(text = 'Kabupaten/Kota')))  # Judul
## Warning in RColorBrewer::brewer.pal(N, "Set2"): n too large, allowed maximum for palette Set2 is 8
## Returning the palette you asked for with that many colors
## Warning in RColorBrewer::brewer.pal(N, "Set2"): n too large, allowed maximum for palette Set2 is 8
## Returning the palette you asked for with that many colors
#Luas Lahan Panen
plot_ly(data, 
        x = ~factor(Tahun), 
        y = ~X3, 
        type = 'scatter', 
        text = ~paste('Kabupaten/Kota: ', `Kab/Kota`, '<br>Value: ', X3),  # Menampilkan informasi saat hover
        hoverinfo = 'text',  # Tampilkan info saat hover
        mode = "marker",
        color = ~factor(`Kab/Kota`))%>%
  layout(title = "Luas Lahan Panen",
         xaxis = list(title = "Tahun"),  # Setel ticks x-axis
         yaxis = list(title = "Luas Lahan Panen (ha)"),
         legend = list(title = list(text = 'Kabupaten/Kota')))  # Judul
## Warning in RColorBrewer::brewer.pal(N, "Set2"): n too large, allowed maximum for palette Set2 is 8
## Returning the palette you asked for with that many colors
## Warning in RColorBrewer::brewer.pal(N, "Set2"): n too large, allowed maximum for palette Set2 is 8
## Returning the palette you asked for with that many colors
#Jumlah Penduduk
plot_ly(data, 
        x = ~factor(Tahun), 
        y = ~X4, 
        type = 'scatter', 
        text = ~paste('Kabupaten/Kota: ', `Kab/Kota`, '<br>Value: ', X4),  # Menampilkan informasi saat hover
        hoverinfo = 'text',  # Tampilkan info saat hover
        mode = "marker",
        color = ~factor(`Kab/Kota`))%>%
  layout(title = "Jumlah Penduduk",
         xaxis = list(title = "Tahun"),  # Setel ticks x-axis
         yaxis = list(title = "Jumlah Penduduk (Ribuan)"),
         legend = list(title = list(text = 'Kabupaten/Kota')))  # Judul
## Warning in RColorBrewer::brewer.pal(N, "Set2"): n too large, allowed maximum for palette Set2 is 8
## Returning the palette you asked for with that many colors
## Warning in RColorBrewer::brewer.pal(N, "Set2"): n too large, allowed maximum for palette Set2 is 8
## Returning the palette you asked for with that many colors

Model Regresi Data Panel

library(plm)
## Warning: package 'plm' was built under R version 4.3.3
## 
## Attaching package: 'plm'
## The following objects are masked from 'package:dplyr':
## 
##     between, lag, lead
com =  plm(Y~X1+X2+X3+X4, data=data,
              model = "pooling")
summary(com)
## Pooling Model
## 
## Call:
## plm(formula = Y ~ X1 + X2 + X3 + X4, data = data, model = "pooling")
## 
## Balanced Panel: n = 24, T = 4, N = 96
## 
## Residuals:
##       Min.    1st Qu.     Median    3rd Qu.       Max. 
## -3.4244664 -0.5229990  0.0031608  0.5684066  2.7513192 
## 
## Coefficients:
##                Estimate  Std. Error t-value  Pr(>|t|)    
## (Intercept) -3.0947e+00  9.9449e-01 -3.1119  0.002484 ** 
## X1           1.3662e-04  4.5376e-05  3.0108  0.003372 ** 
## X2           1.0181e-01  1.9175e-02  5.3093 7.728e-07 ***
## X3           5.2840e-05  2.9808e-06 17.7270 < 2.2e-16 ***
## X4          -1.7591e-02  5.2532e-03 -3.3486  0.001183 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Total Sum of Squares:    594.8
## Residual Sum of Squares: 109.14
## R-Squared:      0.81651
## Adj. R-Squared: 0.80844
## F-statistic: 101.233 on 4 and 91 DF, p-value: < 2.22e-16
fix = plm(Y~X1+X2+X3+X4, data=data,
              model = "within",effect = "individual")
summary(fix)
## Oneway (individual) effect Within Model
## 
## Call:
## plm(formula = Y ~ X1 + X2 + X3 + X4, data = data, effect = "individual", 
##     model = "within")
## 
## Balanced Panel: n = 24, T = 4, N = 96
## 
## Residuals:
##       Min.    1st Qu.     Median    3rd Qu.       Max. 
## -0.6091676 -0.0575424  0.0066621  0.0957061  0.7251311 
## 
## Coefficients:
##       Estimate  Std. Error t-value  Pr(>|t|)    
## X1 -2.7168e-05  1.4176e-05 -1.9165    0.0595 .  
## X2  4.0377e-02  6.8021e-03  5.9359 1.100e-07 ***
## X3  7.6676e-05  7.1239e-06 10.7632 2.481e-16 ***
## X4 -5.8216e-03  4.8224e-03 -1.2072    0.2315    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Total Sum of Squares:    13.159
## Residual Sum of Squares: 3.5556
## R-Squared:      0.7298
## Adj. R-Squared: 0.62251
## F-statistic: 45.9161 on 4 and 68 DF, p-value: < 2.22e-16
fixef(fix)
##                 Bantaeng                    Barru                     Bone 
##                  0.29546                  1.38423                 -2.39993 
##                Bulukumba                 Enrekang                     Gowa 
##                  1.15619                  0.39325                  2.71763 
##                Jeneponto        Kepulauan Selayar                     Luwu 
##                  1.31011                 -0.58476                  1.02511 
##               Luwu Timur               Luwu Utara                 Makassar 
##                  1.27591                  1.01096                 10.85962 
##                    Maros                   Palopo Pangkajene dan Kepulauan 
##                  1.15148                 -0.41876                  1.03918 
##                 Parepare                  Pinrang        Sindereng Rappang 
##                 -0.35133                  1.46087                  2.00090 
##                   Sinjai                  Soppeng                  Takalar 
##                  0.86150                  2.28612                  0.87871 
##              Tana Toraja             Toraja Utara                     Wajo 
##                  0.85664                  0.84371                 -0.48860
rand =  plm(Y~X1+X2+X3+X4,
                data = data,
                model = "random")
summary(rand)
## Oneway (individual) effect Random Effect Model 
##    (Swamy-Arora's transformation)
## 
## Call:
## plm(formula = Y ~ X1 + X2 + X3 + X4, data = data, model = "random")
## 
## Balanced Panel: n = 24, T = 4, N = 96
## 
## Effects:
##                   var std.dev share
## idiosyncratic 0.05229 0.22867 0.043
## individual    1.15250 1.07354 0.957
## theta: 0.8941
## 
## Residuals:
##       Min.    1st Qu.     Median    3rd Qu.       Max. 
## -0.7218123 -0.0825882 -0.0079109  0.0854220  0.9375377 
## 
## Coefficients:
##                Estimate  Std. Error z-value  Pr(>|z|)    
## (Intercept) -6.1615e-01  5.5114e-01 -1.1180   0.26358    
## X1          -2.6620e-05  1.5235e-05 -1.7473   0.08059 .  
## X2           4.1758e-02  7.2185e-03  5.7848 7.259e-09 ***
## X3           6.0402e-05  4.7081e-06 12.8293 < 2.2e-16 ***
## X4           3.9820e-04  1.9777e-03  0.2013   0.84043    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Total Sum of Squares:    19.682
## Residual Sum of Squares: 5.6354
## R-Squared:      0.71369
## Adj. R-Squared: 0.7011
## Chisq: 226.832 on 4 DF, p-value: < 2.22e-16
ranef(rand)
##                 Bantaeng                    Barru                     Bone 
##              -0.24889692               1.10281005              -4.11394648 
##                Bulukumba                 Enrekang                     Gowa 
##              -0.37949623              -0.35765526              -0.82118648 
##                Jeneponto        Kepulauan Selayar                     Luwu 
##              -0.25942333              -0.87338544               0.08487769 
##               Luwu Timur               Luwu Utara                 Makassar 
##               0.58775182               0.12943480               2.38432966 
##                    Maros                   Palopo Pangkajene dan Kepulauan 
##              -0.14300949              -1.00679584              -0.17637736 
##                 Parepare                  Pinrang        Sindereng Rappang 
##              -0.74325005               0.83116149               1.94456530 
##                   Sinjai                  Soppeng                  Takalar 
##               0.12265066               2.11350988              -0.03784358 
##              Tana Toraja             Toraja Utara                     Wajo 
##              -0.12292215               0.03453653              -0.05143927
resid(rand,effect = "individu")
##             1             2             3             4             5 
## -0.2036442056 -0.2992225769 -0.2271010028 -0.2769123170  1.2582882709 
##             6             7             8             9            10 
##  1.0716820756  1.4081521255  0.7231520657 -3.9977932753 -3.9195671250 
##            11            12            13            14            15 
## -4.4400173201 -4.2850574079 -0.4112615072 -0.3494980161 -0.3502424743 
##            16            17            18            19            20 
## -0.4242006019 -0.2956418185 -0.3516922900 -0.3727508661 -0.4267628329 
##            21            22            23            24            25 
## -0.9009827260 -0.8208376206 -0.8191719749 -0.7810107278 -0.2474496926 
##            26            27            28            29            30 
## -0.2426212314 -0.2534480828 -0.3059443028 -0.7914328956 -1.5115616657 
##            31            32            33            34            35 
## -0.6446559904 -0.5855165979 -0.0005862334  0.0822101874  0.0933892521 
##            36            37            38            39            40 
##  0.1683484252  0.7193755726  0.7891215748  0.4280729631  0.4411033929 
##            41            42            43            44            45 
## -0.0157556191  0.1913247668  0.2575067053  0.0905357956  2.4881301303 
##            46            47            48            49            50 
##  2.6322437541  2.5605545904  1.9645668840 -0.1187380312 -0.0392642019 
##            51            52            53            54            55 
## -0.0939085928 -0.3266154752 -0.8919683640 -1.2305175239 -1.0199131056 
##            56            57            58            59            60 
## -0.9304625408 -0.1630803431 -0.1730179617 -0.1699867553 -0.2074266109 
##            61            62            63            64            65 
## -0.7430386690 -0.7481296721 -0.9245335809 -0.5910194518  0.8354478959 
##            66            67            68            69            70 
##  1.1404361487  1.1664353438  0.2200362733  1.7001602361  2.0567890238 
##            71            72            73            74            75 
##  2.6958897495  1.4136468655  0.1243521114  0.1667772158  0.1637529221 
##            76            77            78            79            80 
##  0.0412850270  2.2370838777  2.1209058625  2.4233896887  1.7685497479 
##            81            82            83            84            85 
## -0.0271672799 -0.0528385349 -0.0310845689 -0.0420008784 -0.2684102871 
##            86            87            88            89            90 
## -0.0129634132  0.0167385905 -0.2326304326 -0.0578826350  0.0453166482 
##            91            92            93            94            95 
##  0.2082666826 -0.0559876523 -0.6110948446  0.1530732693  0.2969141204 
##            96 
## -0.0469834273

Model Terbaik

Pada pemilihan model terbaik, terpilih model tetap yang paling baik digunakan berdasarkan uji Chow dan uji Hausman

pFtest(fix,com)
## 
##  F test for individual effects
## 
## data:  Y ~ X1 + X2 + X3 + X4
## F = 87.795, df1 = 23, df2 = 68, p-value < 2.2e-16
## alternative hypothesis: significant effects
phtest(fix,rand)
## 
##  Hausman Test
## 
## data:  Y ~ X1 + X2 + X3 + X4
## chisq = 10.622, df = 4, p-value = 0.03116
## alternative hypothesis: one model is inconsistent

Uji Asumsi

library(pcse)
library(performance)
## Warning: package 'performance' was built under R version 4.3.3
library(moments)
library(lmtest)
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
library(car)
## Loading required package: carData
## 
## Attaching package: 'car'
## The following object is masked from 'package:dplyr':
## 
##     recode
## The following object is masked from 'package:psych':
## 
##     logit
demean <- model.matrix(fix)#Ubah ke matrix
vif(lm(demean[, 1] ~ ., data = as.data.frame(demean)))#Multikol
## Warning in summary.lm(object, ...): essentially perfect fit: summary may be
## unreliable
##       X1       X2       X3       X4 
## 1.232470 1.067082 1.026075 1.220629
check_collinearity(fix)#Multikol (Bagus)
pwartest(fix)#Autokorelasi(Bagus)
## 
##  Wooldridge's test for serial correlation in FE panels
## 
## data:  fix
## F = 1.238, df1 = 1, df2 = 70, p-value = 0.2697
## alternative hypothesis: serial correlation
bptest(formula = fix,varformula = ~X1+X2+X3+X4,data = data)#heteroskedastisitas (bagus)
## 
##  studentized Breusch-Pagan test
## 
## data:  fix
## BP = 4.35, df = 4, p-value = 0.3607
#Normalitas
re = as.numeric(fix$residuals)
jarque.test(re)#Normalitas (Tidak Normal)
## 
##  Jarque-Bera Normality Test
## 
## data:  re
## JB = 28.941, p-value = 5.195e-07
## alternative hypothesis: greater
shapiro.test(residuals(fix))#Normalitas (Tidak Normal)
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(fix)
## W = 0.92075, p-value = 2.23e-05
hist(fix$residuals,xlab = "Residual",col = "blue")