Project ini merupakan Project Kelompok Mata Kuliah Kapita Selekta Statistika
Kehidupan manusia modern saat ini tidak terlepas dari berbagai jenis makanan yang salah satunya adalah cokelat. Cokelat dihasilkan dari biji buah kakao yang telah mengalami serangkaian proses pengolahan sehingga bentuk dan aromanya seperti yang terdapat di pasaran. Biji buah kakao (cokelat) yang telah difermentasi dijadikan serbuk yang disebut cokelat bubuk. Cokelat dalam bentuk bubuk ini banyak dipakai sebagai bahan untuk membuat berbagai macam produk makanan dan minuman, seperti susu, selai, roti, permen dan lain–lain. [1] Sebuah perusahaan yang memproduksi permen coklat terkenal di dunia memiliki banyak hektar pohon yang memproduksi buah kakao di Ghana dan Pantai Gading. Pantai Gading dan Ghana merupakan dua negara yang menyumbang lebih dari 60% produksi kakao global. [2]
Sejumlah penelitian telah dilakukan untuk mempelajari perkembangan kakao dari proses bunga mekar hingga dipanen dengan tujuan meningkatkan hasil panen buah pod kakao. Petani perlu memperbaiki hasil produksi kakao yang diperolehnya, dengan cara melihat faktor yang paling mempengaruhi produksi kakao secara signifikan. Berdasarkan penelitian sebelumya, ada tujuh faktor yang mempengaruhi produksi kakao. Faktor tersebut adalah curah hujan, luah lahan, pengendalian hama, pemangkasan, jumlah tenaga kerja, pengendalian gulma, dan jumlah pupuk. [3] Maka dari itu kami tertarik untuk membahas lebih lanjut faktor-faktor apa saja yang mempengaruhi produksi kakao secara signifikan dengan tujuan meningkatkan hasil panen buah pod kakao dan secara efisien menjelaskan panen masa depan agar berhasil mengantisipasi jumlah bahan yang diperlukan untuk menyiapkan produk perusahaan permen coklat tersebut. Salah satu analisis statistik yang dapat digunakan untuk permasalahan ini adalah analisisa regresi linear atau yang biasa disebut dengan istilah Ordinary Least Squares (OLS) regression.
Analisis regresi linier merupakan analisis yang digunakan untuk melihat hubungan variabel terikat dengan variabel bebas secara linier (garis lurus). Namun karena penelitian ini menggunakan satu variabel terikat dan lebih dari satu variabel bebas dengan jenis variabel bebas terdiri dari variabel kategorik dan variabel numerik, maka digunakanlah analisis regresi berganda dummy.
Pemilihan model dalam penelitian ini menggunakan metode backward analisis regresi dengan memilih model dengan nilai AIC terendah. Model regresi yang telah diperoleh dapat digunakan setelah melakukan uji asumsi terhadap model regresi linier dan deteksi terhadap pencilan. Asumsi klasik yang diuji yaitu normalitas, kehomogenan ragam, autokorelasi, dan multikolinieritas.
Data yang digunakan dalam penelitian ini adalah data kakao dalam file format csv. Dimana data berisi informasi tentang hasil pohon untuk 244 pohon berbeda yang diambil sampelnya di sebuah lokasi di Pantai Gading selama periode 18 tahun (1977-1994), selama bulan-bulan panen (antara September dan Februari). Variabel yang menjadi perhatian adalah pengukuran banyaknya pod buah yang dipanen (nama variabel “hv”) dan dicatat setiap bulan. Berikut merupakan gambaran umum dari variabel yang dicatat di lokasi bersama dengan hasil untuk analisis.
# Input packages yang digunakan dan input data
library(MASS)
library(tseries)
## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoo
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
library(dplyr)
##
## Attaching package: 'dplyr'
## The following object is masked from 'package:car':
##
## recode
## The following object is masked from 'package:MASS':
##
## select
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(rcompanion)
## TREE YEAR MONTH hv blk wilt loss can flu flw ty sml lrg pcp maxt mint
## 1 3 84 10 24 0 1 2 2 1 2 5 3 11 7.95 31.64 21.04
## 2 4 89 1 47 0 0 0 3 1 2 3 5 7 0.00 32.98 17.61
## 3 5 94 11 32 0 4 0 4 1 2 1 3 4 1.63 32.13 20.78
## 4 6 87 11 90 24 20 2 2 1 3 2 4 37 0.45 33.62 21.47
## 5 7 87 11 30 2 8 6 3 1 2 0 2 10 0.45 33.62 21.47
## 6 8 81 11 22 0 0 0 2 1 1 0 0 2 1.05 32.19 20.23
## midt smi lmi Monthv2
## 1 26.34 0.98 0.94 2-Oct
## 2 25.30 0.67 0.21 5-Jan
## 3 26.45 0.94 0.76 3-Nov
## 4 27.55 0.89 0.61 3-Nov
## 5 27.55 0.89 0.61 3-Nov
## 6 26.21 0.91 0.72 3-Nov
Metode penelitian yang digunakan adalah analisis deskriptif dan analisis regresi berganda dummy. Analisis deskriptif digunakan untuk memberikan gambaran karakteristik distribusi data dari variabel independen. Analisis regresi digunakan untuk memodelkan hubungan antara variabel independen dan variabel dependen. Karena variabel dalam komponen independen merupakan tipe data numerik dan kategorik, maka digunakan regresi berganda dummy. Adapun langkah-langkah yang dilakukan, yaitu:
Statistik deskriptif dilakukan untuk mengetahui nilai minimum, maksimum, mean, dan median dari setiap variabel yang ada dan akan di cek normalitas datanya.
#Cek Normalitas data
shapiro.test(kakao$hv)
##
## Shapiro-Wilk normality test
##
## data: kakao$hv
## W = 0.86769, p-value = 1.107e-13
plot(density(kakao$hv))
#Statistika Deskriptif
str(kakao)
## 'data.frame': 244 obs. of 20 variables:
## $ TREE : int 3 4 5 6 7 8 9 10 11 12 ...
## $ YEAR : int 84 89 94 87 87 81 83 86 85 94 ...
## $ MONTH : int 10 1 11 11 11 11 11 10 12 11 ...
## $ hv : int 24 47 32 90 30 22 72 28 19 56 ...
## $ blk : int 0 0 0 24 2 0 8 6 0 0 ...
## $ wilt : int 1 0 4 20 8 0 0 5 0 2 ...
## $ loss : int 2 0 0 2 6 0 2 0 0 0 ...
## $ can : int 2 3 4 2 3 2 2 3 3 4 ...
## $ flu : int 1 1 1 1 1 1 1 3 1 1 ...
## $ flw : int 2 2 2 3 2 1 3 3 2 2 ...
## $ ty : int 5 3 1 2 0 0 3 1 3 1 ...
## $ sml : int 3 5 3 4 2 0 0 1 3 0 ...
## $ lrg : int 11 7 4 37 10 2 1 12 4 3 ...
## $ pcp : num 7.95 0 1.63 0.45 0.45 1.05 1.43 4.24 0.29 1.63 ...
## $ maxt : num 31.6 33 32.1 33.6 33.6 ...
## $ mint : num 21 17.6 20.8 21.5 21.5 ...
## $ midt : num 26.3 25.3 26.4 27.6 27.6 ...
## $ smi : num 0.98 0.67 0.94 0.89 0.89 0.91 0.84 0.9 0.77 0.94 ...
## $ lmi : num 0.94 0.21 0.76 0.61 0.61 0.72 0.6 0.85 0.41 0.76 ...
## $ Monthv2: chr "2-Oct" "5-Jan" "3-Nov" "3-Nov" ...
summary(kakao)
## TREE YEAR MONTH hv
## Min. : 3.0 Min. :77.00 Min. : 1.000 Min. : 2.00
## 1st Qu.: 75.0 1st Qu.:80.00 1st Qu.:10.000 1st Qu.: 14.00
## Median :173.5 Median :84.00 Median :11.000 Median : 21.50
## Mean :186.8 Mean :84.78 Mean : 9.902 Mean : 25.59
## 3rd Qu.:299.2 3rd Qu.:89.00 3rd Qu.:11.000 3rd Qu.: 32.00
## Max. :395.0 Max. :94.00 Max. :12.000 Max. :102.00
## blk wilt loss can
## Min. : 0.000 Min. : 0.000 Min. : 0.0000 Min. :1.000
## 1st Qu.: 0.000 1st Qu.: 0.000 1st Qu.: 0.0000 1st Qu.:2.000
## Median : 0.000 Median : 0.000 Median : 0.0000 Median :3.000
## Mean : 2.307 Mean : 1.279 Mean : 0.6066 Mean :2.672
## 3rd Qu.: 1.000 3rd Qu.: 2.000 3rd Qu.: 1.0000 3rd Qu.:3.000
## Max. :58.000 Max. :20.000 Max. :12.0000 Max. :5.000
## flu flw ty sml
## Min. :1.000 Min. :1.000 Min. : 0.000 Min. : 0.00
## 1st Qu.:1.000 1st Qu.:2.000 1st Qu.: 0.000 1st Qu.: 0.00
## Median :1.000 Median :2.000 Median : 1.000 Median : 0.00
## Mean :1.492 Mean :2.033 Mean : 1.799 Mean : 1.27
## 3rd Qu.:2.000 3rd Qu.:2.000 3rd Qu.: 2.250 3rd Qu.: 2.00
## Max. :4.000 Max. :5.000 Max. :24.000 Max. :14.00
## lrg pcp maxt mint
## Min. : 0.000 Min. :0.000 Min. :30.20 Min. :17.61
## 1st Qu.: 0.000 1st Qu.:1.050 1st Qu.:31.92 1st Qu.:20.17
## Median : 1.500 Median :1.390 Median :32.13 Median :20.70
## Mean : 3.668 Mean :2.211 Mean :32.17 Mean :20.41
## 3rd Qu.: 4.000 3rd Qu.:1.960 3rd Qu.:32.92 3rd Qu.:21.10
## Max. :37.000 Max. :9.880 Max. :33.85 Max. :22.24
## midt smi lmi Monthv2
## Min. :24.91 Min. :0.5100 Min. :0.0800 Length:244
## 1st Qu.:25.95 1st Qu.:0.8225 1st Qu.:0.6000 Class :character
## Median :26.24 Median :0.9100 Median :0.7600 Mode :character
## Mean :26.29 Mean :0.8640 Mean :0.6985
## 3rd Qu.:26.70 3rd Qu.:0.9300 3rd Qu.:0.8200
## Max. :27.55 Max. :0.9900 Max. :0.9600
Dengan masing-masing tipe data dari 20 variabel di atas seperti berikut:
Terdapat 20 variabel yang terdiri dari variabel dependen dan variabel independen, dengan variabel dependen yaitu hv, dan variabel independen yaitu blk, wilt, loss, ty, sml, lrg, maxt, mint, midt, smi, lmi, can, flu, dan flw. Secara teori, variabel Tree, Year, Month, dan Monthv2 tidak digunakan karena variabel Tree, Year, dan Month hanya sebagai informasi umum dari urutan pohon dan waktu panen, sedangkan variabel Monthv2 berbentuk karakter maka tidak dapat digunakan. Untuk mengetahui variabel independen mana saja yang paling signifikan untuk digunakan di dalam model, kita dapat menggunakan uji multikolinearitas dan korelasi untuk mengeliminasi variabel-variabel yang tidak signifikan.
Pengujian Multikolinearitas dilakukan untuk mengetahui korelasi antara variabel independen dengan cara melihat nilai VIF. Hipotesis:
#Model Semua Variabel
models <- lm(hv ~ blk + wilt + loss + ty + sml + lrg + mint + midt + maxt + lmi + smi + pcp + as.factor(flu) + as.factor(flw) + as.factor(can), data = kakao)
#Cek Multikolinieritas
car::vif(models)
## GVIF Df GVIF^(1/(2*Df))
## blk 1.996335 1 1.412917
## wilt 1.802879 1 1.342713
## loss 1.769218 1 1.330119
## ty 2.446138 1 1.564014
## sml 2.806031 1 1.675121
## lrg 2.205639 1 1.485139
## mint 33100.617257 1 181.935750
## midt 47551.054788 1 218.062043
## maxt 18218.793662 1 134.977012
## lmi 11.406148 1 3.377299
## smi 6.141494 1 2.478204
## pcp 3.454294 1 1.858573
## as.factor(flu) 2.431701 3 1.159627
## as.factor(flw) 2.658288 4 1.129992
## as.factor(can) 2.822992 4 1.138515
Variabel yang memiliki nilai VIF>10 artinya ada korelasi antar variabel independen yang lain atau variabel-variabel tersebut memiliki informasi yang sama, maka variabel mint, midt, maxt, dan lmi akan dieliminasi dan tidak akan digunakan sebagai model analisis regresi.
#Buang Variabel yang mengandung multikolinieritas
kakao1 = select(kakao, -mint, -midt, -maxt)
Pengujian Korelasi digunakan untuk mengetahui ada atau tidaknya hubungan antara variabel independen dengan variabel dependen. Hipotesis:
Dengan melihat nilai p-value maka kita dapat menentukan variabel independen mana yang tidak memiliki hubungan signifikan dengan variabel dependen. Jika nilai p-value>0.05 maka variabel independen tidak memiliki hubungan signifikan dengan variabel dependen, sedangkan jika nilai p-value<0.05 maka variabel independen memiliki hubungan signifikan dengan variabel dependen.
#Uji Korelasi, karena datanya tidak normal maka korelasi spearman
library(GGally)
## Loading required package: ggplot2
## Registered S3 method overwritten by 'GGally':
## method from
## +.gg ggplot2
library(ggpubr)
cor.test(kakao$blk, kakao$hv, method = "spearman")
## Warning in cor.test.default(kakao$blk, kakao$hv, method = "spearman"): Cannot
## compute exact p-value with ties
##
## Spearman's rank correlation rho
##
## data: kakao$blk and kakao$hv
## S = 2136682, p-value = 0.06697
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.1174711
#tidak signifikan
cor.test(kakao$wilt, kakao$hv, method = "spearman")
## Warning in cor.test.default(kakao$wilt, kakao$hv, method = "spearman"): Cannot
## compute exact p-value with ties
##
## Spearman's rank correlation rho
##
## data: kakao$wilt and kakao$hv
## S = 1919273, p-value = 0.001128
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.207269
#signifikan
cor.test(kakao$loss, kakao$hv, method = "spearman")
## Warning in cor.test.default(kakao$loss, kakao$hv, method = "spearman"): Cannot
## compute exact p-value with ties
##
## Spearman's rank correlation rho
##
## data: kakao$loss and kakao$hv
## S = 2246635, p-value = 0.2622
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.0720562
#tidak signifikan
cor.test(kakao$ty, kakao$hv, method = "spearman")
## Warning in cor.test.default(kakao$ty, kakao$hv, method = "spearman"): Cannot
## compute exact p-value with ties
##
## Spearman's rank correlation rho
##
## data: kakao$ty and kakao$hv
## S = 1970465, p-value = 0.003523
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.186125
#signifikan
cor.test(kakao$sml, kakao$hv, method = "spearman")
## Warning in cor.test.default(kakao$sml, kakao$hv, method = "spearman"): Cannot
## compute exact p-value with ties
##
## Spearman's rank correlation rho
##
## data: kakao$sml and kakao$hv
## S = 1730136, p-value = 5.91e-06
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.2853898
#signifikan
cor.test(kakao$lrg, kakao$hv, method = "spearman")
## Warning in cor.test.default(kakao$lrg, kakao$hv, method = "spearman"): Cannot
## compute exact p-value with ties
##
## Spearman's rank correlation rho
##
## data: kakao$lrg and kakao$hv
## S = 1569579, p-value = 1.632e-08
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.3517055
#signifikan
cor.test(kakao$smi, kakao$hv, method = "spearman")
## Warning in cor.test.default(kakao$smi, kakao$hv, method = "spearman"): Cannot
## compute exact p-value with ties
##
## Spearman's rank correlation rho
##
## data: kakao$smi and kakao$hv
## S = 2733271, p-value = 0.0442
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## -0.1289424
#signifikan
cor.test(kakao$pcp, kakao$hv, method = "spearman")
## Warning in cor.test.default(kakao$pcp, kakao$hv, method = "spearman"): Cannot
## compute exact p-value with ties
##
## Spearman's rank correlation rho
##
## data: kakao$pcp and kakao$hv
## S = 2707518, p-value = 0.06504
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## -0.1183052
#tidak signifikan
Dari hasil didapat variabel-variabel yang memiliki nilai p-value>0.05 adalah variabel blk, loss, dan pcp sehingga ketiga variabel ini akan dieliminasi dan tidak akan digunakan untuk analisis regresi.
#Buang Variabel yang tidak signifikan
kakao2 = select(kakao1, -blk, -loss, -pcp)
Pemilihan model ini dilakukan untuk mendapatkan model terbaik yang akan digunakan untuk analisis regresi. Pemilihan model dilakukan dengan eliminasi backward yaitu mengeliminasi satu-satu variabel untuk mendapatkan model-model terbaik, kemudian memilih model terbaik dengan menggunakan perbandingan nilai AIC dengan melihat AIC terkecil dari setiap model yang didapatkan. Setelah dilakukan uji Multikolinearitas dan Korelasi maka didapatkan model yang akan digunakan untuk eliminasi backward:
hv = wilt + ty + sml + lrg + smi + as.factor(flu) + as.factor(flw) + as.faktor(can)
Dengan bantuan software R didapatkan hasil eliminasi backward untuk Model 1 sebagai berikut.
#eliminasi backward dimulai dengan model yang kompleks ke model yang lebih sederhana, dengan membuang variabel yang paling tidak signifikan.
model1 <- lm(hv ~ wilt + ty + sml + lrg + smi + as.factor(flu) + as.factor(flw) + as.factor(can), data = kakao2)
dropterm(model1, test = "F")
## Single term deletions
##
## Model:
## hv ~ wilt + ty + sml + lrg + smi + as.factor(flu) + as.factor(flw) +
## as.factor(can)
## Df Sum of Sq RSS AIC F Value Pr(F)
## <none> 51881 1341.7
## wilt 1 708.8 52590 1343.0 3.1013 0.0795754 .
## ty 1 283.1 52164 1341.1 1.2388 0.2668747
## sml 1 7.6 51889 1339.8 0.0332 0.8555561
## lrg 1 3077.7 54959 1353.8 13.4662 0.0003028 ***
## smi 1 1866.5 53748 1348.3 8.1666 0.0046635 **
## as.factor(flu) 3 1006.3 52888 1340.4 1.4677 0.2241301
## as.factor(flw) 4 3789.2 55670 1350.9 4.1448 0.0029243 **
## as.factor(can) 4 1624.9 53506 1341.2 1.7774 0.1342551
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Karena variabel sml memiliki nilai p-value paling tidak signifikan maka variable sml dieliminasi untuk model ke-2.
hv = wilt + ty + lrg + smi + as.factor(flu) + as.factor(flw) + as.faktor(can)
Didapatkan hasil eliminasi backward untuk model 2 sebagai berikut:
#Membuang variabel yang paling tidak signifikan yaitu sml
model2 <- lm(hv ~ wilt + ty + lrg + smi + as.factor(flu) + as.factor(flw) + as.factor(can), data = kakao2)
dropterm(model2, test = "F")
## Single term deletions
##
## Model:
## hv ~ wilt + ty + lrg + smi + as.factor(flu) + as.factor(flw) +
## as.factor(can)
## Df Sum of Sq RSS AIC F Value Pr(F)
## <none> 51889 1339.8
## wilt 1 708.5 52597 1341.1 3.1130 0.0790094 .
## ty 1 472.2 52361 1340.0 2.0748 0.1511250
## lrg 1 3240.5 55129 1352.5 14.2386 0.0002053 ***
## smi 1 1863.1 53752 1346.4 8.1866 0.0046124 **
## as.factor(flu) 3 1001.2 52890 1338.4 1.4664 0.2244739
## as.factor(flw) 4 3834.8 55724 1349.2 4.2125 0.0026102 **
## as.factor(can) 4 1634.1 53523 1339.3 1.7950 0.1306819
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Karena variabel flu memiliki nilai p-value paling tidak signifikan maka variable flu dieliminasi untuk model ke-3.
hv = wilt + ty + lrg + smi + as.factor(flw) + as.faktor(can)
Didapatkan hasil eliminasi backward untuk model 3 sebagai berikut:
#Membuang variabel paling tidak signifikan yaitu flu
model3 <- lm(hv ~ wilt + ty + lrg + smi + as.factor(flw) + as.factor(can), data = kakao2)
dropterm(model3, test = "F")
## Single term deletions
##
## Model:
## hv ~ wilt + ty + lrg + smi + as.factor(flw) + as.factor(can)
## Df Sum of Sq RSS AIC F Value Pr(F)
## <none> 52890 1338.4
## wilt 1 737.7 53628 1339.8 3.2221 0.0739588 .
## ty 1 376.5 53266 1338.2 1.6443 0.2010278
## lrg 1 3655.9 56546 1352.7 15.9672 8.67e-05 ***
## smi 1 2953.6 55844 1349.7 12.8999 0.0004014 ***
## as.factor(flw) 4 3778.9 56669 1347.3 4.1261 0.0030057 **
## as.factor(can) 4 1533.1 54423 1337.4 1.6740 0.1568341
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Karena variabel ty memiliki nilai p-value paling tidak signifikan maka variable ty dieliminasi untuk model ke-4.
hv = wilt + lrg + smi + as.factor(flw) + as.faktor(can)
Didapatkan hasil eliminasi backward untuk model 4 sebagai berikut:
#Membuang variabel paling tidak signifikan yaitu ty
model4 <- lm(hv ~ wilt + lrg + smi + as.factor(flw) + as.factor(can), data = kakao2)
dropterm(model4, test = "F")
## Single term deletions
##
## Model:
## hv ~ wilt + lrg + smi + as.factor(flw) + as.factor(can)
## Df Sum of Sq RSS AIC F Value Pr(F)
## <none> 53266 1338.2
## wilt 1 819.5 54086 1339.9 3.5691 0.0601105 .
## lrg 1 3427.5 56694 1351.4 14.9282 0.0001450 ***
## smi 1 2912.8 56179 1349.2 12.6868 0.0004468 ***
## as.factor(flw) 4 3422.9 56689 1345.3 3.7271 0.0058315 **
## as.factor(can) 4 1367.9 54634 1336.3 1.4894 0.2061419
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Karena variabel can memiliki nilai p-value paling tidak signifikan maka variabel can dieliminasi untuk model ke-5.
hv = wilt + lrg + smi + as.factor(flw)
Didapatkan hasil eliminasi backward untuk model 5 sebagai berikut:
#Membuang variabel paling tidak signifikan yaitu can
model5 <- lm(hv ~ wilt + lrg + smi + as.factor(flw), data = kakao2)
dropterm(model5, test = "F")
## Single term deletions
##
## Model:
## hv ~ wilt + lrg + smi + as.factor(flw)
## Df Sum of Sq RSS AIC F Value Pr(F)
## <none> 54634 1336.3
## wilt 1 627.6 55262 1337.1 2.7109 0.100997
## lrg 1 3874.6 58509 1351.1 16.7367 5.896e-05 ***
## smi 1 2422.7 57057 1344.9 10.4651 0.001391 **
## as.factor(flw) 4 3557.4 58192 1343.7 3.8416 0.004807 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Karena variabel wilt memiliki nilai p-value paling tidak signifikan maka variabel wilt dieliminasi untuk model ke-6.
hv = lrg + smi + as.factor(flw)
Didapatkan hasil eliminasi backward untuk model 6 sebagai berikut:
#Membuang variabel paling tidak signifikan yaitu wilt
model6 <- lm(hv ~ lrg + smi + as.factor(flw), data = kakao2)
dropterm(model6, test = "F")
## Single term deletions
##
## Model:
## hv ~ lrg + smi + as.factor(flw)
## Df Sum of Sq RSS AIC F Value Pr(F)
## <none> 55262 1337.1
## lrg 1 8762.0 64024 1371.0 37.577 3.64e-09 ***
## smi 1 2653.1 57915 1346.6 11.378 0.0008678 ***
## as.factor(flw) 4 3565.2 58827 1344.4 3.822 0.0049591 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Karena nilai p-value semua variabel sudah signifikan maka variabel-variabel sudah signifikan dan tidak dieliminasi lagi.
Dari hasil eliminasi backward telah didapatkan 6 model yang selanjutnya akan dibandingkan nilai AIC masing-masing model.
#Kriteria kebaikan model dengan AIC
AIC=c(AIC(model1),AIC(model2),AIC(model3),AIC(model4),AIC(model5),AIC(model6))
Model=c("Model 1","Model 2","Model 3","Model 4","Model 5","Model 6")
Kriteria=data.frame(Model,AIC)
Kriteria
## Model AIC
## 1 Model 1 2036.170
## 2 Model 2 2034.206
## 3 Model 3 2032.869
## 4 Model 4 2032.600
## 5 Model 5 2030.787
## 6 Model 6 2031.574
Dari perbandingan nilai AIC ke-6 model di tabel maka didapatkan model terbaik adalah model 5 yang memiliki nilai AIC terkecil, maka model 5 yang akan digunakan untuk analisis regresi selanjutnya.
#maka didapat model terbaik adalah:
modelfit <- lm(hv ~ wilt + lrg + smi + as.factor(flw), data = kakao2)
Berikut merupakan model 5 atau model terbaik yang akan dipakai:
hv = wilt + lrg + smi + as.factor(flw)
Terlihat pada model bahwa terdapat satu variabel dengan tipe data skala kategorik yaitu variabel flw, maka selanjutnya model yang akan dibentuk pada saat estimasi parameter adalah model Regresi Linear Berganda Dummy.
Setelah mendapatkan model terbaik dengan melakukan pemilihan model menggunakan eliminasi backward, selanjutnya akan dilakukan uji asumsi klasik untuk melihat apakah model regresi yang dipilih memenuhi asumsi-asumsi dari regresi linear atau tidak.
Uji normalitas bertujuan untuk melihat apakah residual model regresi berdistribusi normal atau tidak. Hipotesis:
Dengan menggunakan bantuan software R dilakukan uji normalitas dengan melihat visualisasi QQ-Plot, histogram, dan metode Kolmogorov-Smirnov.
#Normalitas
residual = resid(modelfit)
ks.test(residual, y= pnorm)
## Warning in ks.test(residual, y = pnorm): ties should not be present for the
## Kolmogorov-Smirnov test
##
## One-sample Kolmogorov-Smirnov test
##
## data: residual
## D = 0.55229, p-value < 2.2e-16
## alternative hypothesis: two-sided
qqPlot ( residual, distribution = "norm", main = "Normal QQ Plot" )
## [1] 60 35
hist(residual)
shapiro.test(residual)
##
## Shapiro-Wilk normality test
##
## data: residual
## W = 0.88111, p-value = 6.866e-13
#Residual model tidak normal
Berdasarkan hasil uji QQ-Plot, histogram, dan Kolmogorov-Smirnov untuk menguji normalitas residual terlihat bahwa pada visualisasi QQ-Plot sebaran residual cenderung tidak mengikuti garis lurus yang ada, sama halnya dengan histogram yang menunjukkan bentuk skewness positif yang mana kedua hasil visualisasi ini menunjukkan residual tidak berdistribusi normal. Selanjutnya hasil dari metode Kolmogorov-Smirnov dimana nilai p-value<2.2e-16 atau kurang dari α=0.05 maka artinya tolak H_0 yang berarti residual tidak berdistribusi normal.
Karena residual model tidak berdistribusi normal maka selanjutnya akan dilakukan transformasi terhadap variabel dependen atau hv.
#Transformasi variabel hv
boxcox(kakao2$hv~1)
p <- powerTransform(kakao2$hv)
hv_trans <- bcPower(kakao2$hv, p$lambda)
#model dengan data hasil transformasi
modelfit2 <- lm(hv_trans ~ wilt + lrg + smi + as.factor(flw), data = kakao2)
summary(modelfit2)
##
## Call:
## lm(formula = hv_trans ~ wilt + lrg + smi + as.factor(flw), data = kakao2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.76480 -0.57931 -0.05234 0.64310 2.56492
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.00885 0.50276 9.963 < 2e-16 ***
## wilt 0.01974 0.03171 0.623 0.534191
## lrg 0.05804 0.01337 4.341 2.11e-05 ***
## smi -1.99403 0.58645 -3.400 0.000791 ***
## as.factor(flw)2 0.64060 0.15482 4.138 4.88e-05 ***
## as.factor(flw)3 0.52305 0.20030 2.611 0.009598 **
## as.factor(flw)4 1.06973 0.34323 3.117 0.002056 **
## as.factor(flw)5 0.28570 0.98515 0.290 0.772066
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.9751 on 236 degrees of freedom
## Multiple R-squared: 0.2269, Adjusted R-squared: 0.2039
## F-statistic: 9.893 on 7 and 236 DF, p-value: 8.107e-11
Dengan menggunakan bantuan software R dilakukan transformasi menggunakan transformasi BoxCox dengan parameter lambda sebesar λ=0.1654258, variabel hv yang telah ditransformasi diberi nama hv_trans.
Selanjutnya diuji ulang normalitas model dengan variabel dependen yang telah ditransformasi, berikut merupakan hasilnya:
#Normalitas
residual = resid(modelfit2)
ks.test(residual, y=pnorm)
## Warning in ks.test(residual, y = pnorm): ties should not be present for the
## Kolmogorov-Smirnov test
##
## One-sample Kolmogorov-Smirnov test
##
## data: residual
## D = 0.060881, p-value = 0.3263
## alternative hypothesis: two-sided
qqPlot (residual, distribution = "norm", main = "Normal QQ Plot" )
## [1] 159 60
hist(residual)
shapiro.test(residual)
##
## Shapiro-Wilk normality test
##
## data: residual
## W = 0.99424, p-value = 0.48
Berdasarkan hasil uji QQ-Plot, histogram, dan Kolmogorov-Smirnov setelah variabel hv ditransformasi terlihat bahwa pada visualisasi QQ-Plot sebaran residual mengikuti garis lurus yang ada, begitu pula dengan histogram yang menunjukkan bentuk lonceng yang mana kedua hasil visualisasi ini menunjukkan residual berdistribusi normal. Selanjutnya hasil dari metode Kolmogorov-Smirnov dimana nilai p-value=0.3263 atau lebih dari α=0.05 maka artinya terima H0 yang berarti residual model berdistribusi normal.
Uji Heteroskedastisitas bertujuan menguji apakah dalam model regresi terjadi ketidaksamaan variansi dari residual satu pengamatan ke pengamatan yang lain. Hipotesis:
#Heteroskedastisitas
bptest(modelfit2,studentize = F, data=kakao2)
##
## Breusch-Pagan test
##
## data: modelfit2
## BP = 7.2199, df = 7, p-value = 0.4063
ncvTest(modelfit2)
## Non-constant Variance Score Test
## Variance formula: ~ fitted.values
## Chisquare = 2.277345, Df = 1, p = 0.13128
Berdasarkan uji Breusch-Pagan dihasilkan nilai p-value=0.4063 dimana nilai ini lebih besar dari α=0.05 maka artinya terima H0 yang berarti tidak terjadi heteroskedastisitas pada residual model.
Uji autokorelasi bertujuan untuk menguji apakah dalam model regresi linear ada korelasi antara residual pada periode sebelumnya. Hipotesis:
#Autokorelasi
durbinWatsonTest(modelfit2)
## lag Autocorrelation D-W Statistic p-value
## 1 0.1073432 1.785158 0.118
## Alternative hypothesis: rho != 0
Berdasarkan uji Durbin-Watson dihasilkan nilai p-value=0.102 dimana nilai ini lebih besar dari α=0.05 maka artinya terima H0 yang berarti tidak terjadi autokorelasi pada residual model.
Uji multikolinearitas bertujuan untuk mengetahui ada tidaknya hubungan antara variabel independen atau bebas. Hipotesis:
#Multikolinearitas
vif(modelfit2)
## GVIF Df GVIF^(1/(2*Df))
## wilt 1.542958 1 1.242159
## lrg 1.622320 1 1.273703
## smi 1.064681 1 1.031834
## as.factor(flw) 1.127056 4 1.015063
Dari hasil uji multikolinearitas di atas, dengan melihat nilai VIF semua variabel berada di bawah nilai 10 maka dapat dikatakan bahwa tidak terjadi multikolinearitas di antara variabel independen.
Berdasarkan uji asumsi klasik model regresi diperoleh bahwa semua asumsi telah terpenuhi. Selanjutnya model regresi linear yang sudah diuji dapat digunakan dan diestimasi parameternya, seperti pada penjelasan sebelumnya model yang akan dibentuk adalah model Regresi Linear Berganda Dummy.
summary(modelfit2)
##
## Call:
## lm(formula = hv_trans ~ wilt + lrg + smi + as.factor(flw), data = kakao2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.76480 -0.57931 -0.05234 0.64310 2.56492
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.00885 0.50276 9.963 < 2e-16 ***
## wilt 0.01974 0.03171 0.623 0.534191
## lrg 0.05804 0.01337 4.341 2.11e-05 ***
## smi -1.99403 0.58645 -3.400 0.000791 ***
## as.factor(flw)2 0.64060 0.15482 4.138 4.88e-05 ***
## as.factor(flw)3 0.52305 0.20030 2.611 0.009598 **
## as.factor(flw)4 1.06973 0.34323 3.117 0.002056 **
## as.factor(flw)5 0.28570 0.98515 0.290 0.772066
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.9751 on 236 degrees of freedom
## Multiple R-squared: 0.2269, Adjusted R-squared: 0.2039
## F-statistic: 9.893 on 7 and 236 DF, p-value: 8.107e-11
Berdasarkan hasil estimasi parameter model regresi linear di atas, karena model yang akan dibentuk adalah regresi linear berganda dengan variabel dummy maka terlihat bahwa variabel flw terbagi menjadi beberapa variabel dengan keterangan tambahan di belakang variabelnya berupa angka 2 sampai 5 yang berarti variabel flw dengan angka 1 sudah otomatis oleh R dijadikan kategori baseline pada model ini.
Dari hasil estimasi parameter dengan R di atas juga terlihat informasi mengenai R-Squared atau Koefisien Determinasi yang bermakna sebagai nilai informasi untuk memprediksi seberapa besar kontribusi pengaruh yang diberikan variabel X secara simultan (bersama-sama) terhadap variabel Y. Nilai R-Squared yang dihasilkan pada model ini adalah sebesar 0.2039 atau sekitar 20% pengaruh yang diberikan variabel X secara simultan terhadap variabel Y. Karena nilai R-Squared ini dianggap masih terlalu kecil, maka penguji akan mencoba melihat apakah terdapat outlier di dalam data dan membuang outlier tersebut.
outlier <- rstandard(modelfit2)
outlier[which(abs(outlier)>2)]
## 7 35 36 60 70 77 115 153
## 2.395678 2.376564 2.408535 2.644249 -2.287392 2.527110 2.185291 -2.611770
## 159 171 190 206 235
## -2.898529 2.096124 -2.408565 -2.410418 -2.135084
Berdasarkan hasil identifikasi outlier, terlihat bahwa ada beberapa outlier pada data ke-7, 35, 36, 60, 70, dst. Maka dari itu data-data yang teridentifikasi sebagai outlier ini selanjutnya akan dibuang dan model akan diuji ulang dimulai dari uji asumsi klasik sampai estimasi parameter model.
# Model ulang
library(tidyverse)
## -- Attaching packages --------------------------------------- tidyverse 1.3.1 --
## v tibble 3.1.4 v purrr 0.3.4
## v tidyr 1.1.3 v stringr 1.4.0
## v readr 2.0.1 v forcats 0.5.1
## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()
## x dplyr::recode() masks car::recode()
## x dplyr::select() masks MASS::select()
## x purrr::some() masks car::some()
kakao1 <- kakao1 %>%
mutate(hv_trans)
kakao3 <- kakao1[-c(7, 35, 36, 60, 70, 77, 115, 153, 159, 171, 190, 206, 235), ]
kakao3<- data.frame(kakao3)
modelfit3 <- lm(hv_trans ~ wilt + lrg + smi + as.factor(flw), data = kakao3)
summary(modelfit3)
##
## Call:
## lm(formula = hv_trans ~ wilt + lrg + smi + as.factor(flw), data = kakao3)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.84063 -0.53566 -0.04674 0.64525 1.92071
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.021142 0.439882 11.415 < 2e-16 ***
## wilt 0.008429 0.027472 0.307 0.759278
## lrg 0.055249 0.011474 4.815 2.72e-06 ***
## smi -1.930625 0.512188 -3.769 0.000210 ***
## as.factor(flw)2 0.576107 0.136793 4.212 3.68e-05 ***
## as.factor(flw)3 0.442843 0.175043 2.530 0.012100 *
## as.factor(flw)4 1.044565 0.293706 3.557 0.000459 ***
## as.factor(flw)5 0.228547 0.836772 0.273 0.785006
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.8276 on 223 degrees of freedom
## Multiple R-squared: 0.2613, Adjusted R-squared: 0.2381
## F-statistic: 11.27 on 7 and 223 DF, p-value: 3.316e-12
dropterm(modelfit3, test = "F")
## Single term deletions
##
## Model:
## hv_trans ~ wilt + lrg + smi + as.factor(flw)
## Df Sum of Sq RSS AIC F Value Pr(F)
## <none> 152.73 -79.582
## wilt 1 0.0645 152.79 -81.485 0.0941 0.7592776
## lrg 1 15.8781 168.60 -58.734 23.1842 2.716e-06 ***
## smi 1 9.7307 162.46 -67.314 14.2081 0.0002096 ***
## as.factor(flw) 4 15.5490 168.28 -65.186 5.6759 0.0002282 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Hipotesis:
#Normalitas
residual = resid(modelfit3)
ks.test(residual, y=pnorm)
## Warning in ks.test(residual, y = pnorm): ties should not be present for the
## Kolmogorov-Smirnov test
##
## One-sample Kolmogorov-Smirnov test
##
## data: residual
## D = 0.064494, p-value = 0.2918
## alternative hypothesis: two-sided
qqPlot (residual, distribution = "norm", main = "Normal QQ Plot" )
## 137 18
## 130 17
hist(residual)
shapiro.test(residual)
##
## Shapiro-Wilk normality test
##
## data: residual
## W = 0.98963, p-value = 0.09626
Berdasarkan hasil uji QQ-Plot, histogram, dan Kolmogorov-Smirnov setelah outlier pada data dibuang terlihat bahwa pada visualisasi QQ-Plot sebaran residual mengikuti garis lurus yang ada, begitu pula dengan histogram yang menunjukkan bentuk lonceng yang mana kedua hasil visualisasi ini menunjukkan residual berdistribusi normal. Selanjutnya hasil dari metode Kolmogorov-Smirnov dimana nilai p-value=0.2918 atau lebih dari α=0.05 maka artinya terima H0 yang berarti residual model berdistribusi normal.
Hipotesis:
#Heteroskedastisitas
bptest(modelfit3,studentize = F, data=kakao3)
##
## Breusch-Pagan test
##
## data: modelfit3
## BP = 3.8028, df = 7, p-value = 0.8022
ncvTest(modelfit3)
## Non-constant Variance Score Test
## Variance formula: ~ fitted.values
## Chisquare = 1.32298, Df = 1, p = 0.25006
Berdasarkan uji Breusch-Pagan dihasilkan nilai p-value=0.8022 dimana nilai ini lebih besar dari α=0.05 maka artinya terima H0 yang berarti tidak terjadi heteroskedastisitas pada residual model
Hipotesis:
#Autokorelasi
durbinWatsonTest(modelfit3)
## lag Autocorrelation D-W Statistic p-value
## 1 0.07599639 1.847851 0.274
## Alternative hypothesis: rho != 0
Berdasarkan uji Durbin-Watson dihasilkan nilai p-value=0.27 dimana nilai ini lebih besar dari α=0.05 maka artinya terima H0 yang berarti tidak terjadi autokorelasi pada residual model.
Hipotesis:
#Multikolinearitas
vif(modelfit3)
## GVIF Df GVIF^(1/(2*Df))
## wilt 1.551559 1 1.245616
## lrg 1.628858 1 1.276267
## smi 1.060936 1 1.030018
## as.factor(flw) 1.135322 4 1.015991
Dari hasil uji multikolinearitas di atas, dengan melihat nilai VIF semua variabel berada di bawah nilai 10 maka dapat dikatakan bahwa tidak terjadi multikolinearitas di antara variabel independen.
Berdasarkan uji ulang asumsi klasik model regresi yang datanya telah dibuang outliernya diperoleh bahwa semua asumsi telah terpenuhi.
Selanjutnya model akan diestimasi ulang parameternya dengan bentuk model regresi berganda variabel dummy yang sama seperti sebelumnya.
summary(modelfit3)
##
## Call:
## lm(formula = hv_trans ~ wilt + lrg + smi + as.factor(flw), data = kakao3)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.84063 -0.53566 -0.04674 0.64525 1.92071
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.021142 0.439882 11.415 < 2e-16 ***
## wilt 0.008429 0.027472 0.307 0.759278
## lrg 0.055249 0.011474 4.815 2.72e-06 ***
## smi -1.930625 0.512188 -3.769 0.000210 ***
## as.factor(flw)2 0.576107 0.136793 4.212 3.68e-05 ***
## as.factor(flw)3 0.442843 0.175043 2.530 0.012100 *
## as.factor(flw)4 1.044565 0.293706 3.557 0.000459 ***
## as.factor(flw)5 0.228547 0.836772 0.273 0.785006
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.8276 on 223 degrees of freedom
## Multiple R-squared: 0.2613, Adjusted R-squared: 0.2381
## F-statistic: 11.27 on 7 and 223 DF, p-value: 3.316e-12
Dari hasil estimasi parameter dengan model ini dihasilkan nilai R-Squared atau Koefisien Determinasi adalah sebesar 0.2381 atau sekitar 23% pengaruh yang diberikan variabel X secara simultan terhadap variabel Y. Nilai R-Squared pada model ini terlihat lebih baik daripada model sebelumnya karena menjelaskan pengaruh variabel X terhadap variabel Y lebih besar. Selanjutnya melihat nilai statistik F atau uji Anova, hasil menunjukkan bahwa nilai p-value=3.316e-12 dimana nilai ini lebih kecil dari α=0.05 maka artinya tolak H0 pada hipotesis Anova yang berarti minimal ada satu variabel X nyata pengaruhnya terhadap variabel Y.
Setelah dilakukan pembentukan model dan estimasi parameternya, selanjutnya akan dilakukan uji Goodness of Fit untuk mengetahui apakah model yang sudah dipilih ini layak untuk digunakan atau tidak. Hipotesis:
#Kriteria Kebaikan Model dengan hoslem.test
library(ResourceSelection)
## ResourceSelection 0.3-5 2019-07-22
modeldata <- model.frame(modelfit3)
hoslem.test(modeldata$hv,fitted(modelfit3))
##
## Hosmer and Lemeshow goodness of fit (GOF) test
##
## data: modeldata$hv, fitted(modelfit3)
## X-squared = -1.047, df = 8, p-value = 1
Berdasarkan hasil di atas dapat dilihat bahwa nilai p-value=1 dimana nilai ini lebih besar dari α=0.05 maka artinya terima H0 yang berarti model dapat diterima
Setelah terpilihnya model regresi linear berganda dummy, selanjutnya akan diinterpretasi hasil estimasi parameter-parameter yang telah diketahui. Berikut penjelasan interpretasi parameter dari masing-masing variabel:
hv_trans = 5.021142 + 0.008429(wilt) + 0.055249(lrg) - 1.930625(smi) + 0.576107(〖flw〗2) + 0.442843(〖flw〗3) + 1.044565(〖flw〗4) + 0.228547(〖flw〗5)
Berdasarkan hasil analisis yang telah dilakukan pada data hasil panen kakao di Ghana dan Pantai Gading selama periode pengamatan 18 tahun (1977 – 1994) maka dapat disimpulkan bahwa:
hv_trans = 5.021142 + 0.008429(wilt) + 0.055249(lrg) - 1.930625(smi) + 0.576107(〖flw〗2) + 0.442843(〖flw〗3) + 1.044565(〖flw〗4) + 0.228547(〖flw〗5)