Adityo Sanjaya Profile
Pada September s/d Oktober 2015, muncul kebakaran hutan di beberapa provinsi di Indonesia. Hal ini menyebabkan munculnya kabut asap di beberapa provinsi tersebut dan beberapa negara tetangga. Dengan menggunakan dua variables yang paling menggambarkan dugaan tersebut, luas kebakaran hutan dan luas perkebunan kelapa sawit, saya menemukan terdapat kemungkinan bahwa kebakaran hutan ditujukan untuk pembesaran kelapa sawit. Saya menemukan bahwa besarnya luas kebakaran hutan di tahun sebelumnya di setiap provinsi, secara rata-rata memprediksi besarnya lahan perkebunan pada masa sekarang.
Berikut adalah visualisasi kebakaran hutan dan persebaaran asap di Indonesia contoh visualisasi
Analisis ini dibagi menjadi empat bagian:
Saya menggunakan data dari 2 sumber:
Kedua data tersebut merupakan data longitudinal. Saya memilih subset dari sample yang memiliki nilai individu x time terbanyak untuk mendapatkan maximal data, dan meminimalisir pengaruh unbalanced data. Saya mendapatkan data sebagai berikut:
In [136]:
DataSawitBersih <- read.csv(file="Data Sawit Bersih.csv",head=TRUE,sep=",")
DataSawitBersih
Out[136]:
Tabel 1. Data luas kebakaran hutan dan luas perkebunan kelapa sawit dalam satu provinsi
Kemudian saya menggunakan angka sebagai penanda provinsi sesuai gambar 1, e.g, 1 adalah Jambi, 2 adalah Kalimantan Barat, dst, sebagai berikut:
In [137]:
DataSawitUnbalanced <- read.csv(file="DataSawitUnbalanced.csv",head=TRUE,sep=",")
DataSawitUnbalanced
Out[137]:
Tabel 2. Data luas kebakaran hutan pada periode sekarang dan luas lahan sawit pada tahun berikutnya.
Pada EDA, saya mencoba untuk menjelaskan karakteristik pada data menggunakan data visualization.
In [138]:
summary(DataSawitUnbalanced$Api)
Out[138]:
In [139]:
summary(DataSawitUnbalanced$Sawit)
Out[139]:
In [140]:
Variance.Api <- var(DataSawitUnbalanced$Api)
Variance.Sawit <- var(DataSawitUnbalanced$Sawit)
Covariance <- cov(DataSawitUnbalanced$Api, DataSawitUnbalanced$Sawit)
plotvar <- cbind(Variance.Api, Variance.Sawit, Covariance)
plotvar
Out[140]:
Kemudian saya mengubah skala dari variabel luas lahan kebakaran dan luas lahan kelapa sawit mengikuti z-distribution, seperti yang disarankan oleh Gelman Gelman (2008) pdf, namun menggunakan hanya satu standard deviation
\begin{equation} X_i standardized = \left(\frac{X_i-\mu}{\sigma}\right). \end{equation}
Sehingga saya mendapatkan variabel yang memiliki rata-rata yang sama, yang distandarisasi dengan standard deviation. Sehingga mendapatkan data dengan nilai mean yang sama sebagai berikut:
In [141]:
Luas_Api <- scale(DataSawitUnbalanced$Api, center = T, scale = T)
Luas_Sawit_t1 <- scale(DataSawitUnbalanced$Sawit, center = T, scale = T)
In [72]:
summary(Luas_Api)
Out[72]:
In [68]:
summary(Luas_Sawit_t1)
Out[68]:
In [142]:
Variance.Api1 <- var(Luas_Api)
Variance.Sawit1 <- var(Luas_Sawit_t1)
Covariance1 <- cov(Luas_Api, Luas_Sawit_t1)
Correlation1 <- cor(Luas_Api, Luas_Sawit_t1)
plotvar1 <- cbind(Variance.Api1, Variance.Sawit1, Covariance1, Correlation1 )
plotvar1
Out[142]:
In [166]:
Api_ts <- ts(Luas_Api)
Sawit_ts <- ts(Luas_Sawit_t1)
Pseudo_time_series <- cbind( Sawit_ts, Api_ts)
# plotting data, dianggap time series
#plot(Pseudo_time_series, col = "red")
Gambar 1. Pola variabel kebakaran hutan dan luas lahan kelapa sawit.
Penemuan:
Berikut adalah plotting dari scatterplot:
In [165]:
#qplot(Luas_Api, Luas_Sawit_t1)
gambar 2. Scatterplot Luas Kebakaran Hutan dan Luas Lahan Kelapa Sawit
Kesimpulan:
Pada bagian ini saya menggunakan metode estimasi yang sederhana namun salah. Hal ini dikarenakan metode yang saya gunakan dalam menganggap data adalah cross section. Namun hal ini memberikan gambaran bahwa terdapat pola yang sama, yaitu Luas kebakaran berpengaruh pada besarnya luas lahan kelapa sawit di periode mendatang. Lagi, saya mengigatkan bahwa ini adalah metode yang salah.
Pada bagian ini saya menggunakan metode OLS cross section untuk mengestimasi pengaruh luas kebakaran hutan pada luas lahan kelapa sawit.
metode OLS yang digunakan:
\begin{equation} y=\beta_0 +\beta_1 x +\varepsilon,\quad i=1,\dots,n.\! \end{equation}yang mana:
y = Luas Lahan kelapa sawit (periode berikutnya)
x = Luas Lahan hutan terbakar (periode sekarang)
In [167]:
OLS <- lm(Luas_Sawit_t1 ~ Luas_Api )
summary(OLS)
Out[167]:
Ups, Luas_Api signifikan mempengaruhi Luas_sawit
In [106]:
require(foreign)
require(MASS)
Pada bagian ini saya menggunakan metode Robust Regression dengan iterated re-weighted least squares cross section untuk mengestimasi pengaruh luas kebakaran hutan pada luas lahan kelapa sawit. Metode ini digunakan karena adanya kemungkinan outlier mengakibatkan beta Api biased up.
metode Robust Regression yang digunakan:
\begin{equation} y=\beta_0 +\beta_1 x +\varepsilon,\quad i=1,\dots,n.\! \end{equation}yang mana:
y = Luas Lahan kelapa sawit (periode berikutnya)
x = Luas Lahan hutan terbakar (periode sekarang)
In [109]:
Robust_Regression <- rlm(Luas_Sawit_t1 ~ Luas_Api )
summary(Robust_Regression)
Out[109]:
In [ ]:
In [131]:
library(plm)
panel_plm <- plm(Sawit ~ Api , model = "between", data = DataSawitUnbalanced)
summary(panel_plm)
Out[131]:
In [134]:
panel_plmb <- plm(Sawit ~ Api , model = "between", data = DataSawitUnbalanced)
panel_plmw <- plm(Sawit ~ Api , model = "within", data = DataSawitUnbalanced)
summary(panel_plmb)
summary(panel_plmw)
Out[134]:
Out[134]: