1. Mencari rata-rata

sum(data$Profit)/length(data$Profit)

## [1] 28.6569

mean(data$Profit)

## [1] 28.6569

2. Mencari median

length(data$Profit) 

## [1] 9994

Karena jumlah data kita genap, maka median adalah data ke 0.5*((N/2) + (N/2 + 1))

Untuk mencari median, data harus diurutkan dulu dari terkecil menjadi terbesar!

N = 9994
letakMedian_kiri = N/2
letakMedian_kanan = N/2+1
0.5*(sort(data$Profit)[letakMedian_kiri] + sort(data$Profit)[letakMedian_kanan])

## [1] 8.6665

median(data$Profit)

## [1] 8.6665

3. Mean vs Median

Dari profit yang kita hitung di atas, kita dapatkan dua jenis ukuran sentral, yaitu:

• 1. Mean = 28.6569
• 2. Median = 8.6665

Manakah ukuran yang lebih representatif?

Untuk membantu menjawab ini, mari kita pertimbangkan skenario lain yang lebih mudah dipahami.

Sekarang misalkan kita kumpulkan gaji 10 orang di sebuah perusahaan unicorn: (dalam juta rupiah)

gaji <- c(5, 6, 7, 5, 6, 7, 5, 6, 7, 7)
mean(gaji)

## [1] 6.1

median(gaji)

## [1] 6

Dari skenario di atas, besaran mean dan median dari sampel tersebut kurang lebih sama, sehingga apa yang kita pilih tidak terlalu masalah.

Sekarang misalkan data kita modifikasi sedikit, yaitu dengan mengikutsertakan gaji seorang CEO. Cukup satu orang saja!

gaji <- c(5, 6, 7, 5, 6, 7, 5, 6, 7, 7, 100) # misalkan gaji CEO 100 juta
mean(gaji)

## [1] 14.63636

median(gaji)

## [1] 6

Perhatikan bagaimana sekarang nilai mean melonjak tinggi, sedangkan nilai median TETAP SAMA.

Apa yang terjadi? Apakah kita adil dengan mengatakan bahwa gaji rata-rata di perusahaan unicorn ini 14.6 juta? Tentu tidak! Faktanya jika tidak menyertakan gaji CEO, gaji rata-rata hanyalah 6.1 juta-an!

Contoh di atas merupakan sebuah contoh bagaimana ukuran mean sangat sensitif terhadap outlier.

Sedangkan ukuran median tidak sensitif terhadap outlier. Untuk kasus di atas (seperti pengukuran gaji), maka nilai median lebih representatif. (Tentunya nanti bergantung kebutuhan kita apa)

Pesan moral: Berhati-hati ketika membaca berita yang melaporkan rata-rata, karena sangat berpotensi menimbulkan bias!

Mari kita tinjau kembali nilai mean dan median dari data kita.

mean(data$Profit)

## [1] 28.6569

median(data$Profit)

## [1] 8.6665

Di dalam statistika, ada modifikasi dari mean yang lebih tidak sensitif terhadap outlier, salah satu contohnya adalah Trimmed Mean. Dalam Trimmed Mean, kita akan membuang nilai ekstrim terbesar dan nilai ekstrim terkecil, atau dalam bahasa lain, kita hanya mengambil 95% yang berada di tengah-tengah data.

mean(data$Profit, trim = 0.05)

## [1] 19.16937

Angka yang dihasilkan Trimmed Mean tidak se-ekstrim mean biasa.

4. Modus (Mode)

Ukuran ini digunakan untuk variabel yang bersifat diskrit atau non-numerikal. Karena di R tidak ada fungsi bawaan untuk menghitung mode, maka kita perlu mendefinisikan fungsi kita sendiri.

getMode <- function(x) {
keys <- unique(x)
keys[which.max(tabulate(match(x, keys)))]
}

getMode(data$Month)

## [1] "November"

5. Menghitung simpanan baku dan variance

Simpangan baku dan variance menghitung sejauh mana nilai dalam distribusi berbeda satu sama lain. Dalam bahasa yang lebih sederhana: “seberapa menyimpang data ini dari nilai mean-nya”?

variance = sum((data$Profit - mean(data$Profit))^2) / (length(data$Profit))
sd = sqrt(variance)
variance

## [1] 54872.31

sd

## [1] 234.2484

variance = sum((data$Profit - mean(data$Profit))^2) / (length(data$Profit)-1)
sd = sqrt(variance)
variance

## [1] 54877.8

sd

## [1] 234.2601

var(data$Profit)

## [1] 54877.8

sd(data$Profit) # sqrt(var(data$Profit))

## [1] 234.2601

# Jika data adalah populasi
gaji <- c(5, 6, 7, 8, 9, 10, 5, 6, 8)
variance = sum((gaji - mean(gaji))^2) / (length(gaji))
sd = sqrt(variance)
variance

## [1] 2.765432

sd

## [1] 1.662959

# Jika data adalah sampel
gaji <- c(5, 6, 7, 8, 9, 10, 5, 6, 8)
variance = sum((gaji - mean(gaji))^2) / (length(gaji) - 1)
sd = sqrt(variance)
variance

## [1] 3.111111

sd

## [1] 1.763834

jumlahJam <- c(24, 24, 24, 24, 24, 24, 24)
var(jumlahJam)

## [1] 0

sd(jumlahJam)

## [1] 0

6. Menghitung range, IQR, quantile

Kita juga bisa menggunakan macam-macam ukuran statistik lainnya, misalnya: range, IQR, dan quantile. Ukuran-ukuran ini umumnya digunakan untuk menyusun boxplot.

min(data$Profit)

## [1] -6599.978

max(data$Profit)

## [1] 8399.976

range(data$Profit) # mengembalikkan min dan max sekaligus

## [1] -6599.978  8399.976

IQR (Interquantile Range) adalah selisih data Q3 dan Q1. Median adalah Q2.

• Q1: 25% data
• Q2: 50% data
• Q3: 75% data

median(data$Profit)

## [1] 8.6665

as.numeric(quantile(data$Profit, 0.50))

## [1] 8.6665

IQR(data$Profit)

## [1] 27.63525

as.numeric(quantile(data$Profit, 0.75) - quantile(data$Profit, 0.25))

## [1] 27.63525

Meskipun kita dapat menggunakan quantile() untuk mendapatkan Q3 dan Q1 satu per satu, kedua angka ini juga disajikan bersama dengan kuratil ke-0 (min()), Q2 (median()), dan ke-100 (max()).

Satu set angka ini dinamakan five-number summary, yang dapat dipanggil menggunakan fungsi fivenum() atau summary()

fivenum(data$Profit)

## [1] -6599.9780     1.7280     8.6665    29.3640  8399.9760

# Five-number summary + mean
summary(data$Profit)

##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
## -6599.978     1.729     8.666    28.657    29.364  8399.976

Kita coba menggunakan variabel lain:

# Five-number summary + mean
summary(data$Ship.Duration)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.000   3.000   4.000   3.958   5.000   7.000

Five-number summary ini juga dapat disajikan dalam bentuk boxplot.

boxplot(data$Ship.Duration)

Exercise 1: Saham manakah yang tingkat volatilitasnya paling tinggi?

saham.A <- c(9999, 10200, 10500, 9800, 9900, 10500, 10900, 11000, 11500, 10500)
saham.B <- c(10000, 9500, 9900, 11000, 12000, 9800, 11500, 10500, 11000, 11500)
saham.C <- c(10000, 10100, 10200, 10000, 9900, 10000, 9900, 10500, 10300, 10200)

plot(saham.B, type = "l", col = "green")
lines(saham.A, col = "red")
lines(saham.C, col = "blue")

7. Menghitung covariance dan correlation

Misalkan kita memilih dua sampel, X dan Y, yang berukuran sama, lalu kita bisa menghitung covariance-nya, yaitu bagaimana variasi dalam sampel X memiliki keterkaitan dengan variasi dalam sampel Y.

Metrik ini banyak digunakan dalam dunia keuangan, misalnya dalam memilih saham untuk sebuah portofolio.

sum((saham.A - mean(saham.A))*(saham.B - mean(saham.B)))/(length(saham.B)-1)

## [1] 26074.44

cov(saham.A, saham.B)

## [1] 26074.44

Untuk menghitung korelasi.

• Nilai korelasi memiliki rentang -1 <= Cor(X,Y) <= 1
• Cor(X,Y) = 0 artinya tidak ada korelasi sama sekali antara X dan Y
• Cor(X,Y) = 1 atau -1 artinya ada korelasi sempurna antara X dan Y

cor(saham.A, saham.B)

## [1] 0.05685123

cor(saham.A, saham.C)

## [1] 0.6081219

1. Menghitung nilai z

Nilai z adalah berapa standar deviasi skor tersebut dihitung dari mean. Distribusi data dianggap distribusi normal.

Contoh: Tinggi badan pria di Indonesia berdistribusi normal dengan rata-rata 160cm dan simpangan baku 7cm. Berapa peluang seorang pria yang dipilih secara acak lebih tinggi dari 175cm?

z <- (175-160)/7
pnorm(z, lower.tail = F)

## [1] 0.01606229

Jika rata-rata 160cm, maka berapa peluang seorang pria yang dipilih secara acak lebih tinggi dari 180cm?

z <- (180-160)/7
pnorm(z, lower.tail = F)

## [1] 0.002137367

Latihan: Misalkan kita memiliki rata-rata populasi 180cm dan simpangan baku 8cm, berapa probabilitas bahwa orang yang dipilih secara acak lebih pendek dari 174cm?

2. Menghitung interval kepercayaan

Pernah mendengar tingkat kepercayaan alfa = 0.05?

# Confidence interval 95%
qnorm(0.025)

## [1] -1.959964

# Confidence interval 99%
qnorm(0.005)

## [1] -2.575829

3. Uji Hipotesis dan p-value

Contoh: * a. Dari data sensus (populasi), misalkan pendapatan rata-rata per satu keluarga adalah 10 juta dengan simpanan baku 3 juta. * b. Sekarang kita melakukan sampling terhadap 30 keluarga, kita dapatkan rata-rata pendapatan per satu keluarga adalah 11 juta (sampel)

Pertanyaannya adalah: Apakah fluktuasi nilai pendapatan rata-rata yang ada di sampel kita ini hanya merupakan kebetulan (by chance karena randomness) atau memang keluarga ini berbeda dari keluarga pada umumnya?

Sekarang, kita asumsikan kita menggunakan tingkat kepercayaan 95%. Kita susun dulu hipotesis pengujian ini: * H0: mean sampel = mean sensus * H1: mean sampel TIDAK SAMA DENGAN mean sensus (two-sided test) * H1: mean sampel LEBIH KECIL atau LEBIH BESAR dari mean sensus (one-sided test)

z_stat <- (11 - 10) / (3 / sqrt(30))
z_stat

## [1] 1.825742

# P-value untuk uji two-sided
(1-pnorm(z_stat))*2

## [1] 0.06788915

Nilai z adalah 1.825. Nilai p-value 0.0679. Hasil TIDAK signifikan pada p < .05.

# P-value untuk uji one-sided
(1-pnorm(z_stat))

## [1] 0.03394458

Nilai z adalah 2.11. Nilai p-value 0.0339. Hasil signifikan pada p < .05.

4. Uji T-Test

Umumnya, uji-z digunakan ketika kita memiliki ukuran sampel yang cukup besar (aturan praktisnya adalah n >= 30) dan ketika simpangan baku populasi diketahui. Jika kondisi di atas tidak terpenuhi, kita dapat menggunakan pengukuran statistik yang dikenal sebagai uji-t Student.

Contoh: Diketahui 10 gaji karyawan dari perusahaan unicorn. Kita tahu dari data resmi bahwa rata-rata gaji adalah 10 juta, tapi kita tidak ketahui simpangan bakunya.

gaji <- c(9,8,8,9,10,11,12,8,7,9)
mean(gaji)

## [1] 9.1

• H0: mean gaji sampel SAMA dengan rata-rata data resmi
• H1: mean gaji sampel LEBIH KECIL dari rata-rata data resmi (one-sided test)

t.test(gaji, mu=10, alternative="less")

## 
##  One Sample t-test
## 
## data:  gaji
## t = -1.8676, df = 9, p-value = 0.04733
## alternative hypothesis: true mean is less than 10
## 95 percent confidence interval:
##      -Inf 9.983367
## sample estimates:
## mean of x 
##       9.1

• H0: mean gaji sampel SAMA dengan rata-rata data resmi
• H1: mean gaji sampel TIDAK SAMA dari rata-rata data resmi (two-sided test)

t.test(gaji, mu=10, alternative = "two.sided")

## 
##  One Sample t-test
## 
## data:  gaji
## t = -1.8676, df = 9, p-value = 0.09466
## alternative hypothesis: true mean is not equal to 10
## 95 percent confidence interval:
##   8.009879 10.190121
## sample estimates:
## mean of x 
##       9.1

Jika nilai p kurang dari 0,05, kita menolak H0, maka kita bisa menyimpulkan bahwa memang ada perbedaan yang signifikan. Jika nilai p lebih besar dari 0,05, kita tidak dapat menyimpulkan bahwa ada perbedaan yang signifikan.

1. skim()

Secara pribadi, saya menyukai fungsi skim() dari paket skimr daripada glimpse() atau summary() karena memberikan informasi yang lebih rinci tentang data, termasuk missing data, jumlah data unik, jenis kolom dan untuk jenis numerik, skim() juga mengembalikan ringkasan lima angka dengan histogram kecil.

skim(data)

Variable type: character

Variable type: Date

Variable type: numeric

2. ggpairs()

Dengan menggunakan ggpairs(), kita bisa membuat plot dari beberapa variabel sekaligus.

ggpairs(data[c('Sales','Profit','Quantity','Discount')])

3. Matriks Korelasi

Membuat matriks korelasi

cor_mat <- cor(data[c('Sales','Profit','Quantity','Discount')])
cor_mat

##                Sales      Profit   Quantity    Discount
## Sales     1.00000000  0.47906435 0.20079477 -0.02819012
## Profit    0.47906435  1.00000000 0.06625319 -0.21948746
## Quantity  0.20079477  0.06625319 1.00000000  0.00862297
## Discount -0.02819012 -0.21948746 0.00862297  1.00000000

Membuat plot dari matriks korelasi

corrplot(cor_mat, method="pie", type="lower", addCoef.col = "black")

4. Modelling

model_single <- lm(Profit ~ Sales, data=data)
get_regression_summaries(model_single)

## # A tibble: 1 x 9
##   r_squared adj_r_squared    mse  rmse sigma statistic p_value    df  nobs
##       <dbl>         <dbl>  <dbl> <dbl> <dbl>     <dbl>   <dbl> <dbl> <dbl>
## 1      0.23         0.229 42279.  206.  206.     2976.       0     1  9994

5. Menguji Asumsi Normalitas

Salah satu asumsi dalam regresi linier adalah bahwa residual berdistribusi normal. Sekarang kita akan memplot Studentized Residuals model kita.

srs <- studres(model_single)
hist(srs, freq=FALSE, 
     main="Distribution of Studentized Residuals",
     breaks=100)
xfit<-seq(min(srs),max(srs),length=40) 
yfit<-dnorm(xfit) 
lines(xfit, yfit)

Tampaknya data kita tidak berdistribusi normal dengan melihat grafik. Untuk memastikannya, kita harus menghitungnya secara statistik. Karena data kita melebihi 5000, kita harus menggunakan uji Kolmogorov-Smirnov.

ks.test(studres(model_single), data$Profit)

## Warning in ks.test(studres(model_single), data$Profit): p-value will be
## approximate in the presence of ties

## 
##  Two-sample Kolmogorov-Smirnov test
## 
## data:  studres(model_single) and data$Profit
## D = 0.76926, p-value < 2.2e-16
## alternative hypothesis: two-sided

• H0: Data mengikuti distribusi yang ditentukan
• HA: Data tidak mengikuti distribusi yang ditentukan

Karena p-value kita <0,5, maka H0 ditolak, sehingga data kita tidak terdistribusi normal.

6. Uji Heteroskedastisitas

Uji heteroskedastisitas bertujuan untuk menguji apakah dalam model regresi terjadi ketidaksamaan varian dari residual satu pengamatan ke pengamatan yang lain.

bptest(model_single)

## 
##  studentized Breusch-Pagan test
## 
## data:  model_single
## BP = 2723, df = 1, p-value < 2.2e-16

• H0 : Residual bersifat Homoskedastisitas
• HA : Residual bersifat Heteroskedastisitas

Karena p-value kita < 0,5, maka H0 ditolak, sehingga data error kita heteroskedastis, artinya data kita tidak terlalu cocok untuk menggunakan model regresi ini.

par(mfrow = c(2,2))
plot(model_single)

Ini adalah visualisasi dari residual model kita Dengan melihat Residuals vs Fitted, intinya terlihat seperti kipas, artinya semakin besar nilainya maka semakin besar pula residualnya. Kita harus mencoba untuk menormalkan mereka. Selain itu, dengan melihat Residuals vs Leverage, kita dapat melihat beberapa poin outlier yang akan kita lakukan analisis outlier.