Data

df <- read_excel(path = "Data/data_02.xlsx", col_names = TRUE, sheet="sheet1")

head(df)

## # A tibble: 6 x 8
##   provinsi tahun ipm_p   idg keluhan_p ahh_p exp_p tidak_sekolah_p
##   <chr>    <dbl> <dbl> <dbl>     <dbl> <dbl> <dbl>           <dbl>
## 1 Aceh      2010  63.5  53.4      37.0  71.1  8.70            5.98
## 2 Aceh      2011  64.1  52.1      32.1  71.2  8.75            4.67
## 3 Aceh      2012  65.0  54.4      32.6  71.2  8.79            4.75
## 4 Aceh      2013  65.7  59.8      30.5  71.3  8.83            4.3 
## 5 Aceh      2014  66.9  65.1      32.0  71.3  8.89            3.76
## 6 Aceh      2015  67.5  65.6      30.3  71.5  8.91            4.19

dim(df)

## [1] 340   8

df$tahun <- as.factor(df$tahun )

str(df)

## tibble [340 x 8] (S3: tbl_df/tbl/data.frame)
##  $ provinsi       : chr [1:340] "Aceh" "Aceh" "Aceh" "Aceh" ...
##  $ tahun          : Factor w/ 10 levels "2010","2011",..: 1 2 3 4 5 6 7 8 9 10 ...
##  $ ipm_p          : num [1:340] 63.5 64.1 65 65.7 66.9 ...
##  $ idg            : num [1:340] 53.4 52.1 54.4 59.8 65.1 ...
##  $ keluhan_p      : num [1:340] 37 32.1 32.6 30.5 32 ...
##  $ ahh_p          : num [1:340] 71.1 71.2 71.2 71.3 71.3 ...
##  $ exp_p          : num [1:340] 8.7 8.75 8.79 8.83 8.89 ...
##  $ tidak_sekolah_p: num [1:340] 5.98 4.67 4.75 4.3 3.76 4.19 3.47 2.99 2.21 1.96 ...

skimr::skim(df)

Variable type: character

Variable type: factor

Variable type: numeric

df_numerik <- df[,3:8]
boxplot(df_numerik)

par(mfrow=c(2,3))
boxplot(df$ipm_p, xlab = "ipm_p")
boxplot(df$idg, xlab = "idg")
boxplot(df$keluhan_p, xlab = "keluhan_p")
boxplot(df$ahh_p, xlab = "ahh_p")
boxplot(df$exp_p, xlab = "exp_p")
boxplot(df$tidak_sekolah_p, xlab = "tidak_sekolah_p")

plot_histogram(data = df, nrow = 3, ncol = 2, geom_histogram_args = list (fill="#D27685"))

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

korelasi <- cor(df_numerik)
ggcorrplot(korelasi, type="lower", lab = TRUE)

Admin1 <- "Peta/Admin1Provinsi/idn_admbnda_adm1_bps_20200401.shp"
Admin1 <- st_read(Admin1)

## Reading layer `idn_admbnda_adm1_bps_20200401' from data source 
##   `D:\USK Model Ekonomi R\Praktik\Peta\Admin1Provinsi\idn_admbnda_adm1_bps_20200401.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 34 features and 12 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 95.01079 ymin: -11.00762 xmax: 141.0194 ymax: 6.07693
## Geodetic CRS:  WGS 84

data2010 <- df %>% filter(tahun==2010)
Indo_Geo_2010 <- geo_join(spatial_data=Admin1, 
                       data_frame=data2010, by_sp="ADM1_EN", 
                       by_df="provinsi", how = "inner")

## Warning: We recommend using the dplyr::*_join() family of functions instead.

Spasial2010<-ggplot()+
  geom_sf(data=Indo_Geo_2010,aes(fill=Indo_Geo_2010$ipm_p))+theme_bw()
Spasial2010

data2019 <- df %>% filter(tahun==2019)
Indo_Geo_2019 <- geo_join(spatial_data=Admin1, 
                       data_frame=data2019, by_sp="ADM1_EN", 
                       by_df="provinsi", how = "inner")

## Warning: We recommend using the dplyr::*_join() family of functions instead.

Spasial2019<-ggplot()+
  geom_sf(data=Indo_Geo_2019,aes(fill=Indo_Geo_2019$ipm_p))+theme_bw()
Spasial2019

Pengujian Multikolinearitas

Uji multikolinearitas adalah proses untuk memeriksa apakah terdapat ketergantungan yang tinggi antara dua atau lebih peubah independen dalam model regresi. Ketergantungan ini dapat menyebabkan masalah dalam estimasi koefisien regresi dan penafsiran hasil model.

Salah satu cara untuk menguji multikolinearitas adalah dengan menggunakan statistik VIF (Variance Inflation Factor). VIF mengukur tingkat keragaman dari koefisien regresi yang dipengaruhi oleh korelasi antar peubah independen dalam model.Nilai VIF yang tinggi (biasanya di atas 5 atau 10) menunjukkan adanya multikolinearitas yang signifikan dalam model

H0 : VIF < 10 (Tidak Multikolinearitas)

H1 : VIF > 10 (Terdapat Multikolinearitas)

library(car)

## Loading required package: carData

## 
## Attaching package: 'car'

## The following object is masked from 'package:dplyr':
## 
##     recode

model1 <- lm(ipm_p ~ idg + keluhan_p + ahh_p + exp_p + tidak_sekolah_p,df%>%filter(tahun==2010))
model2 <- lm(ipm_p ~ idg + keluhan_p + ahh_p + exp_p + tidak_sekolah_p,df%>%filter(tahun==2011))
model3 <- lm(ipm_p ~ idg + keluhan_p + ahh_p + exp_p + tidak_sekolah_p,df%>%filter(tahun==2012))
model4 <- lm(ipm_p ~ idg + keluhan_p + ahh_p + exp_p + tidak_sekolah_p,df%>%filter(tahun==2013))
model5 <- lm(ipm_p ~ idg + keluhan_p + ahh_p + exp_p + tidak_sekolah_p,df%>%filter(tahun==2014))
model6 <- lm(ipm_p ~ idg + keluhan_p + ahh_p + exp_p + tidak_sekolah_p,df%>%filter(tahun==2015))
model7 <- lm(ipm_p ~ idg + keluhan_p + ahh_p + exp_p + tidak_sekolah_p,df%>%filter(tahun==2016))
model8 <- lm(ipm_p ~ idg + keluhan_p + ahh_p + exp_p + tidak_sekolah_p,df%>%filter(tahun==2017))
model9 <- lm(ipm_p ~ idg + keluhan_p + ahh_p + exp_p + tidak_sekolah_p,df%>%filter(tahun==2018))
model10 <- lm(ipm_p ~ idg + keluhan_p + ahh_p + exp_p + tidak_sekolah_p,df%>%filter(tahun==2019))
model11 <- lm(ipm_p ~ idg + keluhan_p + ahh_p + exp_p + tidak_sekolah_p,df)


Multikol <- rbind(as.vector(vif(model1)),as.vector(vif(model2)),as.vector(vif(model3)),
           as.vector(vif(model4)),as.vector(vif(model5)),as.vector(vif(model6)),
           as.vector(vif(model7)),as.vector(vif(model8)),as.vector(vif(model9)),
           as.vector(vif(model10)),as.vector(vif(model11)))
           
rownames(Multikol) <- c("Tahun 2010","Tahun 2011","Tahun 2012",
                        "Tahun 2013","Tahun 2014","Tahun 2015","Tahun 2016",
                        "Tahun 2017","Tahun 2018","Tahun 2019","Tahun 2010-2019")

colnames(Multikol) <- c("idg","keluhan_p","ahh_p","exp_p","tidak_sekolah_p")

Multikol

##                      idg keluhan_p    ahh_p    exp_p tidak_sekolah_p
## Tahun 2010      1.459847  1.095796 1.350749 1.581056        1.181364
## Tahun 2011      1.377840  1.191890 1.427362 1.761132        1.218923
## Tahun 2012      1.309676  1.255514 1.367289 1.867328        1.221725
## Tahun 2013      1.460422  1.295906 1.274194 1.809225        1.207805
## Tahun 2014      1.072382  1.249225 1.277141 1.542944        1.133811
## Tahun 2015      1.087030  1.226466 1.274655 1.467806        1.153927
## Tahun 2016      1.091343  1.289929 1.373153 1.658655        1.100768
## Tahun 2017      1.158628  1.138590 1.332401 1.498504        1.201334
## Tahun 2018      1.062899  1.127595 1.345526 1.423362        1.161315
## Tahun 2019      1.091915  1.163698 1.349917 1.478440        1.213975
## Tahun 2010-2019 1.213838  1.099840 1.319259 1.575000        1.198021

Berdasarkan hasil VIF model OLS di atas , semua model pada tahun 2010-2019 memiliki VIF < 10 terima H0 artinya untuk semua peubah pada model tidak ada multikolinearitas.

library(plm)

## 
## Attaching package: 'plm'

## The following objects are masked from 'package:dplyr':
## 
##     between, lag, lead

Pemilihan Model Terbaik

Model Common Effect Model

Model Common Effect (MCE) adalah salah satu jenis model regresi panel data yang digunakan untuk memperkirakan efek dari peubah bebas terhadap peubah terikat pada setiap waktu atau periode, serta mengontrol efek konstan atau tetap pada setiap unit individu dalam panel.Dalam MCE, diasumsikan bahwa setiap unit individu memiliki intercept (efek tetap) yang sama di setiap waktu, dan hanya koefisien regresi peubah bebas yang berbeda-beda untuk setiap individu. Dengan kata lain, efek tetap atau konstan pada setiap individu dianggap homogen atau sama.

cem  <- plm(formula = ipm_p ~ idg + keluhan_p + ahh_p + exp_p + tidak_sekolah_p, 
            data = df, model = "pooling")
summary(cem )

## Pooling Model
## 
## Call:
## plm(formula = ipm_p ~ idg + keluhan_p + ahh_p + exp_p + tidak_sekolah_p, 
##     data = df, model = "pooling")
## 
## Balanced Panel: n = 34, T = 10, N = 340
## 
## Residuals:
##      Min.   1st Qu.    Median   3rd Qu.      Max. 
## -3.004731 -0.815794 -0.041955  0.794266  4.603826 
## 
## Coefficients:
##                   Estimate Std. Error  t-value  Pr(>|t|)    
## (Intercept)     -87.040063   2.543425 -34.2216 < 2.2e-16 ***
## idg               0.041261   0.011783   3.5016 0.0005253 ***
## keluhan_p        -0.057360   0.013026  -4.4036 1.436e-05 ***
## ahh_p             0.526113   0.030470  17.2668 < 2.2e-16 ***
## exp_p            12.861383   0.296227  43.4173 < 2.2e-16 ***
## tidak_sekolah_p  -0.198832   0.012978 -15.3212 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Total Sum of Squares:    10106
## Residual Sum of Squares: 573.11
## R-Squared:      0.94329
## Adj. R-Squared: 0.94244
## F-statistic: 1111.17 on 5 and 334 DF, p-value: < 2.22e-16

library(tseries)

## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo

res.cem <- residuals(cem)
(normal <- jarque.bera.test(res.cem))

## 
##  Jarque Bera Test
## 
## data:  res.cem
## X-squared = 4.227, df = 2, p-value = 0.1208

(homos <- bptest(cem))

## 
##  studentized Breusch-Pagan test
## 
## data:  cem
## BP = 46.04, df = 5, p-value = 8.915e-09

(autokol <- pbgtest(cem))

## 
##  Breusch-Godfrey/Wooldridge test for serial correlation in panel models
## 
## data:  ipm_p ~ idg + keluhan_p + ahh_p + exp_p + tidak_sekolah_p
## chisq = 216.59, df = 10, p-value < 2.2e-16
## alternative hypothesis: serial correlation in idiosyncratic errors

Model Fixed Effect Model Within Estimator

Model Fixed Effect (MFE) adalah salah satu jenis model regresi panel data yang digunakan untuk memperkirakan efek dari peubah bebas terhadap peubah terikat pada setiap waktu atau periode, serta mengontrol efek tetap atau konstan pada setiap unit individu dalam panel.Dalam MFE, diasumsikan bahwa efek tetap atau konstan pada setiap unit individu berbeda-beda di setiap waktu atau periode, dan hanya koefisien regresi peubah bebas yang sama untuk setiap individu. Dengan kata lain, efek tetap pada setiap individu dianggap heterogen atau berbeda-beda.

fem.twoway  <- plm(formula = ipm_p ~ idg + keluhan_p + ahh_p + exp_p + tidak_sekolah_p, 
             data = df, model = "within", effect= "twoways", index = c("provinsi","tahun"))
summary(fem.twoway)

## Twoways effects Within Model
## 
## Call:
## plm(formula = ipm_p ~ idg + keluhan_p + ahh_p + exp_p + tidak_sekolah_p, 
##     data = df, effect = "twoways", model = "within", index = c("provinsi", 
##         "tahun"))
## 
## Balanced Panel: n = 34, T = 10, N = 340
## 
## Residuals:
##      Min.   1st Qu.    Median   3rd Qu.      Max. 
## -2.477621 -0.127991  0.010706  0.158409  1.643066 
## 
## Coefficients:
##                   Estimate Std. Error t-value  Pr(>|t|)    
## idg              0.0222521  0.0060834  3.6579 0.0003016 ***
## keluhan_p       -0.0250657  0.0078165 -3.2068 0.0014913 ** 
## ahh_p            0.8545198  0.1037668  8.2350 6.073e-15 ***
## exp_p            6.9323436  0.4032464 17.1913 < 2.2e-16 ***
## tidak_sekolah_p -0.0719896  0.0182223 -3.9506 9.780e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Total Sum of Squares:    90.428
## Residual Sum of Squares: 30.069
## R-Squared:      0.66748
## Adj. R-Squared: 0.61396
## F-statistic: 117.229 on 5 and 292 DF, p-value: < 2.22e-16

res.fem.twoway<- residuals(fem.twoway)
(normal <- jarque.bera.test(res.fem.twoway))

## 
##  Jarque Bera Test
## 
## data:  res.fem.twoway
## X-squared = 3883.1, df = 2, p-value < 2.2e-16

(homos <- bptest(fem.twoway))

## 
##  studentized Breusch-Pagan test
## 
## data:  fem.twoway
## BP = 46.04, df = 5, p-value = 8.915e-09

(autokol <- pbgtest(fem.twoway))

## 
##  Breusch-Godfrey/Wooldridge test for serial correlation in panel models
## 
## data:  ipm_p ~ idg + keluhan_p + ahh_p + exp_p + tidak_sekolah_p
## chisq = 93.55, df = 10, p-value = 1.056e-15
## alternative hypothesis: serial correlation in idiosyncratic errors

Uji Chow

Uji Chow digunakan untuk menguji apakah koefisien regresi sama di seluruh kelompok (individual) atau tidak.

H0 : Model Common Effect Lebih Baik Digunakan

H1 : Model Fixed Effect Lebih Baik Digunakan

pooltest(cem,fem.twoway)

## 
##  F statistic
## 
## data:  ipm_p ~ idg + keluhan_p + ahh_p + exp_p + tidak_sekolah_p
## F = 125.56, df1 = 42, df2 = 292, p-value < 2.2e-16
## alternative hypothesis: unstability

Kesimpulan : Dengan tingkat keyakinan 95%, kita yakin bahwa model fixed effect merupakan model yang lebih baik untuk digunakan bila dibandingkan dengan model common effect.

Model Random Effect Model

Model Random Effect (MRE) adalah salah satu jenis model regresi panel data yang digunakan untuk memperkirakan efek dari peubah bebas terhadap peubah terikat pada setiap waktu atau periode, serta mengontrol efek tetap atau konstan pada setiap unit individu dalam panel.Dalam MRE, diasumsikan bahwa efek tetap atau konstan pada setiap unit individu dan efek regresi peubah bebas pada setiap individu merupakan peubah acak atau random. Dengan kata lain, efek tetap dan efek regresi dianggap sebagai gabungan antara faktor-faktor yang dapat diamati dan faktor-faktor yang tidak diamati.

rem.twoway  <- plm(formula = ipm_p ~ idg + keluhan_p + ahh_p + exp_p + tidak_sekolah_p, 
             data = df, model = "random", effect= "twoways", index = c("provinsi","tahun"))
summary(rem.twoway)

## Twoways effects Random Effect Model 
##    (Swamy-Arora's transformation)
## 
## Call:
## plm(formula = ipm_p ~ idg + keluhan_p + ahh_p + exp_p + tidak_sekolah_p, 
##     data = df, effect = "twoways", model = "random", index = c("provinsi", 
##         "tahun"))
## 
## Balanced Panel: n = 34, T = 10, N = 340
## 
## Effects:
##                   var std.dev share
## idiosyncratic 0.10298 0.32090 0.065
## individual    1.46700 1.21120 0.926
## time          0.01397 0.11818 0.009
## theta: 0.9165 (id) 0.5778 (time) 0.5765 (total)
## 
## Residuals:
##      Min.   1st Qu.    Median   3rd Qu.      Max. 
## -3.222964 -0.255579  0.026134  0.293941  2.621427 
## 
## Coefficients:
##                    Estimate  Std. Error  z-value  Pr(>|z|)    
## (Intercept)     -1.0234e+02  5.9807e+00 -17.1122 < 2.2e-16 ***
## idg              3.5989e-02  8.5043e-03   4.2319 2.318e-05 ***
## keluhan_p       -2.7621e-02  1.0419e-02  -2.6510  0.008025 ** 
## ahh_p            1.0669e+00  9.0440e-02  11.7966 < 2.2e-16 ***
## exp_p            1.0175e+01  4.8647e-01  20.9166 < 2.2e-16 ***
## tidak_sekolah_p -1.5528e-01  2.2480e-02  -6.9074 4.938e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Total Sum of Squares:    394.42
## Residual Sum of Squares: 70.108
## R-Squared:      0.82225
## Adj. R-Squared: 0.81959
## Chisq: 1545.06 on 5 DF, p-value: < 2.22e-16

res.rem.twoway<- residuals(rem.twoway)
(normal <- jarque.bera.test(res.rem.twoway))

## 
##  Jarque Bera Test
## 
## data:  res.rem.twoway
## X-squared = 1412.1, df = 2, p-value < 2.2e-16

(homos <- bptest(rem.twoway))

## 
##  studentized Breusch-Pagan test
## 
## data:  rem.twoway
## BP = 46.04, df = 5, p-value = 8.915e-09

(autokol <- pbgtest(rem.twoway))

## 
##  Breusch-Godfrey/Wooldridge test for serial correlation in panel models
## 
## data:  ipm_p ~ idg + keluhan_p + ahh_p + exp_p + tidak_sekolah_p
## chisq = 105.64, df = 10, p-value < 2.2e-16
## alternative hypothesis: serial correlation in idiosyncratic errors

Uji Haussman

H0 : Cov(Ui,x′it)=0 (Random Effect)

H1 : Cov(Ui,x′it)≠0 (Fixed Effect)

phtest(fem.twoway, rem.twoway)

## 
##  Hausman Test
## 
## data:  ipm_p ~ idg + keluhan_p + ahh_p + exp_p + tidak_sekolah_p
## chisq = 156.6, df = 5, p-value < 2.2e-16
## alternative hypothesis: one model is inconsistent

Kesimpulan : Dengan tingkat keyakinan 95%, kita yakin bahwa model fixed effect merupakan model yang lebih baik untuk digunakan bila dibandingkan dengan model random effect.

Fixed Effect Model

#efek individu dan waktu
plmtest(fem.twoway,type = "bp", effect="twoways" )

## 
##  Lagrange Multiplier Test - two-ways effects (Breusch-Pagan)
## 
## data:  ipm_p ~ idg + keluhan_p + ahh_p + exp_p + tidak_sekolah_p
## chisq = 849.51, df = 2, p-value < 2.2e-16
## alternative hypothesis: significant effects

#efek individual
plmtest(fem.twoway,type = "bp", effect="individual" )

## 
##  Lagrange Multiplier Test - (Breusch-Pagan)
## 
## data:  ipm_p ~ idg + keluhan_p + ahh_p + exp_p + tidak_sekolah_p
## chisq = 750.36, df = 1, p-value < 2.2e-16
## alternative hypothesis: significant effects

#efek waktu
plmtest(fem.twoway,type = "bp", effect="time" )

## 
##  Lagrange Multiplier Test - time effects (Breusch-Pagan)
## 
## data:  ipm_p ~ idg + keluhan_p + ahh_p + exp_p + tidak_sekolah_p
## chisq = 99.154, df = 1, p-value < 2.2e-16
## alternative hypothesis: significant effects

pcdtest(fem.twoway, test = c("lm"), index=NULL,w =NULL )

## 
##  Breusch-Pagan LM test for cross-sectional dependence in panels
## 
## data:  ipm_p ~ idg + keluhan_p + ahh_p + exp_p + tidak_sekolah_p
## chisq = 1364.5, df = 561, p-value < 2.2e-16
## alternative hypothesis: cross-sectional dependence

Fixed Effect Model dengan LSDV

Alasan menggunakan FEM LSDV dibandingkan dengan Within Estimator adalah sebagai berikut:

Konsistensi Estimator: FEM LSDV memberikan hasil yang lebih konsisten daripada Within Estimator ketika terdapat korelasi antara peubah bebas dan peubah tak terobservasi (unobserved variables). FEM LSDV memperhitungkan peubah tak terobservasi melalui dummy variables untuk setiap unit observasi, sehingga dapat menghilangkan efek dari peubah tak terobservasi tersebut pada estimasi koefisien peubah bebas.

Efisiensi Estimator: FEM LSDV juga lebih efisien dalam kasus panel data dengan jumlah unit waktu (time period) yang sedikit dan jumlah unit cross section yang banyak. Dalam kasus ini, Within Estimator cenderung menghasilkan varian yang besar dan bias yang signifikan.

Uji Hipotesis: FEM LSDV juga memungkinkan penggunaan uji hipotesis untuk perbedaan antar unit cross section (individual effects). Dalam FEM LSDV, uji hipotesis untuk perbedaan antar unit cross section dapat dilakukan melalui uji F.

Model FE yang akan dibentuk pada kasus ini adalah model FE dengan jumlah dummy k-1 dengan k adalah jumlah kategori. Hal ini dilakukan untuk mencegah adanya dummy trap yaitu suatu kondisi di mana dua atau lebih peubah dummy yang digunakan dalam model regresi saling berkaitan satu sama lain. Maka dari itu model FE LSDV yang digunakan adalah model dengan interesep.

fem.lsdv <- lm(ipm_p ~ idg + keluhan_p + ahh_p + exp_p + tidak_sekolah_p+factor(tahun)+factor(provinsi) , df)
summary(fem.lsdv)

## 
## Call:
## lm(formula = ipm_p ~ idg + keluhan_p + ahh_p + exp_p + tidak_sekolah_p + 
##     factor(tahun) + factor(provinsi), data = df)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.47762 -0.12799  0.01071  0.15841  1.64307 
## 
## Coefficients:
##                                      Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                        -57.666090   7.193626  -8.016 2.64e-14 ***
## idg                                  0.022252   0.006083   3.658 0.000302 ***
## keluhan_p                           -0.025066   0.007816  -3.207 0.001491 ** 
## ahh_p                                0.854520   0.103767   8.235 6.07e-15 ***
## exp_p                                6.932344   0.403246  17.191  < 2e-16 ***
## tidak_sekolah_p                     -0.071990   0.018222  -3.951 9.78e-05 ***
## factor(tahun)2011                    0.256692   0.086533   2.966 0.003262 ** 
## factor(tahun)2012                    0.679775   0.096654   7.033 1.44e-11 ***
## factor(tahun)2013                    1.124981   0.107769  10.439  < 2e-16 ***
## factor(tahun)2014                    1.610802   0.113956  14.135  < 2e-16 ***
## factor(tahun)2015                    1.984730   0.122508  16.201  < 2e-16 ***
## factor(tahun)2016                    2.029765   0.140904  14.405  < 2e-16 ***
## factor(tahun)2017                    2.515994   0.138561  18.158  < 2e-16 ***
## factor(tahun)2018                    2.828430   0.149897  18.869  < 2e-16 ***
## factor(tahun)2019                    3.041345   0.170313  17.857  < 2e-16 ***
## factor(provinsi)Bali                -0.767800   0.330597  -2.322 0.020895 *  
## factor(provinsi)Bangka Belitung     -3.275910   0.160549 -20.404  < 2e-16 ***
## factor(provinsi)Banten              -1.108178   0.196996  -5.625 4.33e-08 ***
## factor(provinsi)Bengkulu            -0.846599   0.188400  -4.494 1.01e-05 ***
## factor(provinsi)DI Yogyakarta        1.278246   0.548969   2.328 0.020571 *  
## factor(provinsi)DKI Jakarta          1.963663   0.392555   5.002 9.80e-07 ***
## factor(provinsi)Gorontalo           -3.431920   0.316377 -10.848  < 2e-16 ***
## factor(provinsi)Jambi               -3.425772   0.203266 -16.854  < 2e-16 ***
## factor(provinsi)Jawa Barat          -4.371011   0.319793 -13.668  < 2e-16 ***
## factor(provinsi)Jawa Tengah         -4.800365   0.471176 -10.188  < 2e-16 ***
## factor(provinsi)Jawa Timur          -2.966412   0.242573 -12.229  < 2e-16 ***
## factor(provinsi)Kalimantan Barat    -5.303226   0.225155 -23.554  < 2e-16 ***
## factor(provinsi)Kalimantan Selatan  -1.918472   0.264463  -7.254 3.65e-12 ***
## factor(provinsi)Kalimantan Tengah   -2.991689   0.167313 -17.881  < 2e-16 ***
## factor(provinsi)Kalimantan Timur    -2.724696   0.462453  -5.892 1.05e-08 ***
## factor(provinsi)Kalimantan Utara    -3.337941   0.403133  -8.280 4.47e-15 ***
## factor(provinsi)Kepulauan Riau       1.517991   0.253789   5.981 6.46e-09 ***
## factor(provinsi)Lampung             -3.374999   0.151585 -22.265  < 2e-16 ***
## factor(provinsi)Maluku               1.500949   0.487446   3.079 0.002273 ** 
## factor(provinsi)Maluku Utara        -1.939088   0.264370  -7.335 2.20e-12 ***
## factor(provinsi)NTB                 -0.863151   0.503078  -1.716 0.087270 .  
## factor(provinsi)NTT                 -1.830053   0.384843  -4.755 3.12e-06 ***
## factor(provinsi)Papua               -7.342913   0.679165 -10.812  < 2e-16 ***
## factor(provinsi)Papua Barat         -3.686533   0.451805  -8.160 1.01e-14 ***
## factor(provinsi)Riau                -2.404250   0.213717 -11.250  < 2e-16 ***
## factor(provinsi)Sulawesi Barat      -1.602653   0.588464  -2.723 0.006849 ** 
## factor(provinsi)Sulawesi Selatan    -1.165770   0.195872  -5.952 7.59e-09 ***
## factor(provinsi)Sulawesi Tengah     -0.807207   0.285871  -2.824 0.005074 ** 
## factor(provinsi)Sulawesi Tenggara   -2.071030   0.194920 -10.625  < 2e-16 ***
## factor(provinsi)Sulawesi Utara      -1.000168   0.227990  -4.387 1.61e-05 ***
## factor(provinsi)Sumatera Barat       1.025474   0.212869   4.817 2.34e-06 ***
## factor(provinsi)Sumatera Selatan    -2.599964   0.176657 -14.718  < 2e-16 ***
## factor(provinsi)Sumatera Utara      -0.244689   0.218320  -1.121 0.263303    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3209 on 292 degrees of freedom
## Multiple R-squared:  0.997,  Adjusted R-squared:  0.9965 
## F-statistic:  2082 on 47 and 292 DF,  p-value: < 2.2e-16

Interpretasi Hasil