Statistics Homework Solution

Import libraries

library(tidyverse)
library(GGally)
library(glmnet)
library(ggpubr)
library(MASS)

Soal 1

Untuk membandingkan hasil dari 2 program pelatihan pekerja di sebuah perusahaan, dilakukan eksperimen pada 20 orang pekerja. 10 orang secara acak dipilih untuk mengikuti program 1, sisanya dipilih untuk mengikuti program 2. Setelah pelatihan selesai, semua pekerja diuji kemampuannya dalam menyelesaikan sebuah pekerjaan. Data hasil pengujian diberikan sebagai berikut.

Waktu (menit)
Program 1	15	20	11	23	16	21	18	16	27	24
Program 2	23	31	13	19	23	17	28	26	25	28

(a) Tampilkan distribusi dari kedua data di atas dalam bentuk boxplot. Berikan analisis deskriptif dari output yang diperoleh.

Boxplot

program1 <- c(15,20,11,23,16,21,18,16,27,24)
program2 <- c(23,31,13,19,23,17,28,26,25,28) 
data <- data.frame( 
  group = rep(c("Program 1", "Program 2"), each = 10),
  time = c(program1, program2)
)

ggboxplot(data, x = "group", y = "time", color = "group",
          ylab = "Time (mins)", xlab = "")

Descriptive Statistics

summary(program1)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   11.0    16.0    19.0    19.1    22.5    27.0 

summary(program2)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   13.0    20.0    24.0    23.3    27.5    31.0 

Berdasarkan output statistika deskriptif di atas, diperoleh nilai rata-rata dari program 1 dan program 2 secara berturut-turut adalah 19.1 dan 23.3, dengan selisih rata-rata sebesar 4.2. Berdasarkan output boxplot, distribusi data sampel program 1 dan program 2 terlihat berbeda di mana waktu penyelesaian pekerjaan dari sampel program 1 cenderung lebih cepat dibandingkan program 2.

(b) Dari data hasil observasi di atas, dapatkah kita simpulkan bahwa rata-rata waktu penyelesaiaan pekerjaan hasil pelatihan program 1 lebih cepat dibandingkan hasil pelatihan program 2? (Lakukan pengujian dengan \(\alpha = 0.05\))

Sebelum dilakukan uji perbandingan rata-rata (t-test), akan dicek terlebih dahulu apakah varians kedua sampel berbeda menggunakan F-test.

var.test(program1, program2)

    F test to compare two variances

data:  program1 and program2
F = 0.75117, num df = 9, denom df = 9, p-value = 0.6769
alternative hypothesis: true ratio of variances is not equal to 1
95 percent confidence interval:
 0.1865797 3.0242006
sample estimates:
ratio of variances 
         0.7511686 

Berdasarkan output F-test di atas, diperoleh bahwa varians kedua sampel tidak berbeda secara signifikan dengan p-value sebesar 0.6869.

Selanjutnya, dilakukan two sample t-test untuk membandingkan apakah rata-rata waktu penyelesaian pekerjaan hasil pelatihan program 1 lebih cepat dibandingkan hasil pelatihan program 2. Misalkan rata-rata populasi dari waktu pekerjaan program 1 dan program 2 dinotasikan secara berturut-turut dengan \(\mu_1\) dan \(\mu_2\). Hipotesis dari persoalan ini adalah sebagai berikut.

Hipotesis statistik

\(H_0\): \(\mu_1 \geq \mu_2\)

\(H_1\): \(\mu_1 < \mu_2\)

Berdasarkan hasil F-test sebelumnya, maka digunakan asumsi kesamaan varians dalam t-test berikut dengan setting var.equal = TRUE. Hipotesis alternatif yang dipilih adalah \(H_1\): \(\mu_1 < \mu_2\) yang merupakan uji sisi kiri (left-tailed), maka digunakan setting alternative = "less". Berikut adalah hasil perhitungan uji hipotesis t-test dengan menggunakan function t.test().

t.test(program1, program2, alternative = "less", var.equal = TRUE)

    Two Sample t-test

data:  program1 and program2
t = -1.8055, df = 18, p-value = 0.04387
alternative hypothesis: true difference in means is less than 0
95 percent confidence interval:
       -Inf -0.1662568
sample estimates:
mean of x mean of y 
     19.1      23.3 

Berdasarkan output t-test di atas, diperoleh nilai p-value sebesar 0.04387, yang artinya secara statistik rata-rata waktu penyelesaian pekerjaan hasil pelatihan program 1 secara signifikan lebih cepat dibandingkan hasil pelatihan program 2 dengan level signifikansi \(\alpha = 0.05\).

Soal 2

Sebuah developer real estate sedang mempertimbangkan untuk berinvestasi dengan membangun mal di sebuah kota. Terdapat 3 kota (kota A, B, dan C) yang menjadi pertimbangan tempat mal tersebut dibangun. Salah satu pertimbangan pemilihan lokasi mal adalah besarnya pendapatan penduduk di wilayah sekitar mal. Oleh karena itu, dilakukan pengambilan sampel secara acak sebanyak 30 keluarga di tiap kota untuk ditanyakan besarnya pendapatan bulanan keluarga tersebut (dalam satuan juta rupiah). Berdasarkan hasil sampling yang diperoleh, lakukan analisis statistik yang tepat untuk memutuskan di kota mana sebaiknya dibangun mal. (Lakukan pengujian dengan \(\alpha = 0.05\). Data sampling tersedia pada Google Drive folder data dengan nama file mall_investment.csv)

Berdasarkan persoalan di atas, analisis yang sesuai untuk mengambil keputusan pemilihan lokasi pembangunan mal adalah Analysis of Variance (ANOVA). ANOVA digunakan untuk membandingkan rata-rata pendapatan dari kota A, B, dan C. Untuk mengimplementasikan ANOVA dengan R, terlebih dahulu dilakukan import data CSV ke dalam program R.

Load dataset

data <- read_csv("mall_investment.csv")

── Column specification ─────────────────────────────────────────────────────────────────────────────────────────────
cols(
  income = col_double(),
  city = col_character()
)

Berikut tampilan visualisasi data pendapatan dari kota A, B, dan C dalam boxplot.

Boxplot

ggboxplot(data, x = "city", y = "income", color = "city",
          ylab = "Income (IDR juta)", xlab = "")

Berdasarkan output boxplot di atas, dapat dilihat bahwa pendapatan penduduk di Kota C cenderung lebih tinggi dibandingkan Kota A dan B. Selanjutnya dilakukan ANOVA untuk menguji rata-rata pendapatan ketiga kota di atas. Misalkan rata-rata populasi dari pendapatan Kota A, B, dan C secara berturut-turut dinotasikan dengan \(\mu_A\), \(\mu_B\), dan \(\mu_C\). Hipotesis statistik dari persoalan di atas adalah sebagai berikut.

Hipotesis Statistik

\(H_0\): \(\mu_A = \mu_B = \mu_C\)

\(H_1\): Minimal ada \(\mu_i \neq \mu_j\) dengan \(i=A,B,C\), \(j=A,B,C\) dan \(i \neq j\).

Implementasi ANOVA di dalam R dapat menggunakan function aov().

ANOVA

fit <- aov(income ~ city, data=data)
summary(fit)

            Df Sum Sq Mean Sq F value Pr(>F)    
city         2  500.5  250.27   224.5 <2e-16 ***
Residuals   87   97.0    1.11                   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Berdasarkan output ANOVA di atas, diperoleh nilai p-value <2e-16 yang artinya terdapat nilai rata-rata yang berbeda secara signifikan di antara ketiga kota tersebut. Untuk mengecek kota mana yang berbeda, dilakukan uji lanjutan (post-hoc test) dengan menggunakan uji Tukey’s HSD (honestly significant difference). Implementasi Tukey’s HSD di dalam R dapat menggunakan fungsi TukeyHSD().

TukeyHSD(fit)

  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = income ~ city, data = data)

$city
          diff        lwr       upr     p adj
B-A -0.2877863 -0.9377745 0.3622019 0.5440096
C-A  4.8526334  4.2026452 5.5026216 0.0000000
C-B  5.1404196  4.4904314 5.7904078 0.0000000

Berdasarkan output dari uji Tukey’s HSD, diperoleh bahwa:

• Rata-rata kota A dan B tidak berbeda secara signifikan
• Rata-rata kota A dan C berbeda secara signifikan
• Rata-rata kota B dan C berbeda secara signifikan

data %>% group_by(city) %>% summarise(mean=mean(income))

Hasil di atas menunjukkan bahwa rata-rata populasi penduduk Kota C berbeda signifikan dan lebih tinggi dibandingkan Kota A dan Kota B, dengan rata-rata sampel sebesar IDR 15.002 juta. Maka, berdasarkan hasil ANOVA dan Tukey’s HSD, dapat direkomendasikan untuk membangun mal di wilayah Kota C karena penduduk Kota C memiliki rata-rata pendapatan yang lebih tinggi dibandingkan Kota A dan Kota B.

Soal 3

Lakukan analisis regresi pada dataset mtcars yang tersedia di R. Lakukan seleksi variabel untuk memperoleh model terbaik. Periksa juga asumsi residual model.

Data mtcars

data <- mtcars
head(data, 10)

data <- data %>% 
          mutate(cyl=as_factor(cyl), vs=as_factor(vs), am=as_factor(am), gear=as_factor(gear))
head(data,10)

mtcars data visualization

ggpairs(data, columns = 1:6, progress = FALSE)

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

ggpairs(data, columns = c(1,7:11), progress = FALSE)

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

Multiple regression model

m1 <- lm(mpg ~ ., data = data)
summary(m1)

Call:
lm(formula = mpg ~ ., data = data)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.2015 -1.2319  0.1033  1.1953  4.3085 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)  
(Intercept) 15.09262   17.13627   0.881   0.3895  
cyl6        -1.19940    2.38736  -0.502   0.6212  
cyl8         3.05492    4.82987   0.633   0.5346  
disp         0.01257    0.01774   0.708   0.4873  
hp          -0.05712    0.03175  -1.799   0.0879 .
drat         0.73577    1.98461   0.371   0.7149  
wt          -3.54512    1.90895  -1.857   0.0789 .
qsec         0.76801    0.75222   1.021   0.3201  
vs1          2.48849    2.54015   0.980   0.3396  
am1          3.34736    2.28948   1.462   0.1601  
gear4       -0.99922    2.94658  -0.339   0.7382  
gear5        1.06455    3.02730   0.352   0.7290  
carb         0.78703    1.03599   0.760   0.4568  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.616 on 19 degrees of freedom
Multiple R-squared:  0.8845,    Adjusted R-squared:  0.8116 
F-statistic: 12.13 on 12 and 19 DF,  p-value: 1.764e-06

Berdasarkan output model regresi di atas, diperoleh hasil bahwa hanya terdapat 2 variabel yang signifikan berpengaruh terhadap variabel respon mpg. Oleh karena itu, selanjutnya akan dilakukan seleksi variabel dengan menggunakan 2 metode, yaitu Stepwise Regression dan LASSO.

Stepwise Regression

m1.stepwise <- stepAIC(m1, direction = "both", trace = TRUE)

Start:  AIC=70.87
mpg ~ cyl + disp + hp + drat + wt + qsec + vs + am + gear + carb

       Df Sum of Sq    RSS    AIC
- gear  2    5.1061 135.16 68.103
- drat  1    0.9408 130.99 69.101
- disp  1    3.4354 133.49 69.705
- carb  1    3.9503 134.00 69.828
- vs    1    6.5693 136.62 70.447
- qsec  1    7.1353 137.19 70.579
- cyl   2   16.4500 146.50 70.682
<none>              130.05 70.870
- am    1   14.6316 144.68 72.282
- hp    1   22.1573 152.21 73.905
- wt    1   23.6065 153.66 74.208

Step:  AIC=68.1
mpg ~ cyl + disp + hp + drat + wt + qsec + vs + am + carb

       Df Sum of Sq    RSS    AIC
- drat  1     0.025 135.18 66.108
- carb  1     3.866 139.02 67.005
- vs    1     4.035 139.19 67.044
- disp  1     4.732 139.89 67.204
- qsec  1     4.941 140.10 67.251
- cyl   2    14.238 149.40 67.308
<none>              135.16 68.103
- am    1    15.929 151.09 69.668
- hp    1    18.284 153.44 70.163
+ gear  2     5.106 130.05 70.870
- wt    1    31.992 167.15 72.901

Step:  AIC=66.11
mpg ~ cyl + disp + hp + wt + qsec + vs + am + carb

       Df Sum of Sq    RSS    AIC
- vs    1     4.250 139.43 65.099
- carb  1     4.808 139.99 65.227
- disp  1     4.895 140.08 65.247
- qsec  1     4.918 140.10 65.252
- cyl   2    17.095 152.28 65.919
<none>              135.18 66.108
- am    1    16.829 152.01 67.863
+ drat  1     0.025 135.16 68.103
- hp    1    19.891 155.07 68.501
+ gear  2     4.191 130.99 69.101
- wt    1    33.543 168.73 71.201

Step:  AIC=65.1
mpg ~ cyl + disp + hp + wt + qsec + am + carb

       Df Sum of Sq    RSS    AIC
- carb  1     2.898 142.33 63.757
- disp  1     4.214 143.65 64.052
- cyl   2    13.993 153.43 64.160
<none>              139.43 65.099
- qsec  1    10.717 150.15 65.469
+ vs    1     4.250 135.18 66.108
- am    1    14.361 153.79 66.236
- hp    1    15.649 155.08 66.503
+ drat  1     0.241 139.19 67.044
+ gear  2     2.277 137.16 68.572
- wt    1    36.334 175.77 70.510

Step:  AIC=63.76
mpg ~ cyl + disp + hp + wt + qsec + am

       Df Sum of Sq    RSS    AIC
- disp  1     1.651 143.98 62.126
- cyl   2    11.107 153.44 62.162
- qsec  1     8.078 150.41 63.524
<none>              142.33 63.757
- hp    1    15.403 157.73 65.046
+ carb  1     2.898 139.43 65.099
+ vs    1     2.340 139.99 65.227
- am    1    17.424 159.75 65.453
+ drat  1     1.125 141.21 65.503
+ gear  2     3.935 138.40 66.860
- wt    1    40.707 183.04 69.807

Step:  AIC=62.13
mpg ~ cyl + hp + wt + qsec + am

       Df Sum of Sq    RSS    AIC
- cyl   2    16.085 160.07 61.515
- qsec  1     7.044 151.03 61.655
<none>              143.98 62.126
- hp    1    15.443 159.42 63.387
+ vs    1     2.744 141.24 63.511
- am    1    16.566 160.55 63.611
+ disp  1     1.651 142.33 63.757
+ drat  1     0.767 143.22 63.956
+ carb  1     0.334 143.65 64.052
+ gear  2     3.246 140.74 65.397
- wt    1    52.932 196.91 70.145

Step:  AIC=61.52
mpg ~ hp + wt + qsec + am

       Df Sum of Sq    RSS    AIC
- hp    1     9.219 169.29 61.307
<none>              160.07 61.515
+ cyl   2    16.085 143.98 62.126
+ disp  1     6.629 153.44 62.162
+ carb  1     3.227 156.84 62.864
+ drat  1     1.428 158.64 63.229
- qsec  1    20.225 180.29 63.323
+ vs    1     0.249 159.82 63.466
- am    1    25.993 186.06 64.331
+ gear  2     1.764 158.30 65.161
- wt    1    78.494 238.56 72.284

Step:  AIC=61.31
mpg ~ wt + qsec + am

       Df Sum of Sq    RSS    AIC
<none>              169.29 61.307
+ hp    1     9.219 160.07 61.515
+ carb  1     8.036 161.25 61.751
+ disp  1     3.276 166.01 62.682
+ drat  1     1.400 167.89 63.042
+ vs    1     0.000 169.29 63.307
+ cyl   2     9.862 159.42 63.387
- am    1    26.178 195.46 63.908
+ gear  2     0.185 169.10 65.272
- qsec  1   109.034 278.32 75.217
- wt    1   183.347 352.63 82.790

m1.stepwise

Call:
lm(formula = mpg ~ wt + qsec + am, data = data)

Coefficients:
(Intercept)           wt         qsec          am1  
      9.618       -3.917        1.226        2.936  

summary(m1.stepwise)

Call:
lm(formula = mpg ~ wt + qsec + am, data = data)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.4811 -1.5555 -0.7257  1.4110  4.6610 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   9.6178     6.9596   1.382 0.177915    
wt           -3.9165     0.7112  -5.507 6.95e-06 ***
qsec          1.2259     0.2887   4.247 0.000216 ***
am1           2.9358     1.4109   2.081 0.046716 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.459 on 28 degrees of freedom
Multiple R-squared:  0.8497,    Adjusted R-squared:  0.8336 
F-statistic: 52.75 on 3 and 28 DF,  p-value: 1.21e-11

LASSO

y <- as.matrix(mtcars %>% dplyr::select(mpg))
x <- as.matrix(mtcars %>% dplyr::select(-mpg))

cv.lasso <- cv.glmnet(x,y)
plot(cv.lasso)

cv.lasso$lambda.min

[1] 0.8007036

cv.lasso$lambda.1se

[1] 1.685404

coef(cv.lasso, s = "lambda.min")

11 x 1 sparse Matrix of class "dgCMatrix"
                      1
(Intercept) 36.00001676
cyl         -0.88608541
disp         .         
hp          -0.01168438
drat         .         
wt          -2.70814703
qsec         .         
vs           .         
am           .         
gear         .         
carb         .         

coef(cv.lasso, s = "lambda.1se")

11 x 1 sparse Matrix of class "dgCMatrix"
                       1
(Intercept) 32.956547393
cyl         -0.822871058
disp         .          
hp          -0.004701802
drat         .          
wt          -2.202101850
qsec         .          
vs           .          
am           .          
gear         .          
carb         .          

m1.lasso <- glmnet(x,y,lambda = cv.lasso$lambda.min)
m1.lasso$beta

10 x 1 sparse Matrix of class "dgCMatrix"
              s0
cyl  -0.88547684
disp  .         
hp   -0.01169485
drat  .         
wt   -2.70853300
qsec  .         
vs    .         
am    .         
gear  .         
carb  .         

# OLS estimate with selected variables from LASSO
m1.lasso2 <- lm(mpg ~ cyl+hp+wt, data = mtcars)
summary(m1.lasso2)

Call:
lm(formula = mpg ~ cyl + hp + wt, data = mtcars)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.9290 -1.5598 -0.5311  1.1850  5.8986 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 38.75179    1.78686  21.687  < 2e-16 ***
cyl         -0.94162    0.55092  -1.709 0.098480 .  
hp          -0.01804    0.01188  -1.519 0.140015    
wt          -3.16697    0.74058  -4.276 0.000199 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.512 on 28 degrees of freedom
Multiple R-squared:  0.8431,    Adjusted R-squared:  0.8263 
F-statistic: 50.17 on 3 and 28 DF,  p-value: 2.184e-11

Residual Diagnostics

Stepwise

par(mfrow=c(2,2))
plot(m1.stepwise,1)
plot(m1.stepwise,2)

plot(m1.stepwise,3)
plot(m1.stepwise,4)

shapiro.test(m1.stepwise$residuals)

    Shapiro-Wilk normality test

data:  m1.stepwise$residuals
W = 0.9411, p-value = 0.08043

hist(m1.stepwise$residuals, main = "", xlab = "residual")

LASSO

par(mfrow=c(2,2))
plot(m1.lasso2,1)
plot(m1.lasso2,2)

plot(m1.lasso2,3)
plot(m1.lasso2,4)

shapiro.test(m1.lasso2$residuals)

    Shapiro-Wilk normality test

data:  m1.lasso2$residuals
W = 0.93455, p-value = 0.05252

hist(m1.lasso2$residuals, main = "", xlab = "residual")