Komputasi Statistika

Tugas 6

Email             : calvin.riswandi@student.matanuniversity.ac.id
RPubs            : https://rpubs.com/calvinriswandy/
Jurusan          : Statistika Bisnis
Address         : ARA Center, Matana University Tower
                         Jl. CBD Barat Kav, RT.1, Curug Sangereng, Kelapa Dua, Tangerang, Banten 15810.

---

1 Library

setwd(getwd())
library(ggplot2)
library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library(broom)
library(ggpubr)

2 Simple Regression

2.1 import data

income.data = read.csv("income.data.csv")
summary(income.data)

##        X             income        happiness    
##  Min.   :  1.0   Min.   :1.506   Min.   :0.266  
##  1st Qu.:125.2   1st Qu.:3.006   1st Qu.:2.266  
##  Median :249.5   Median :4.424   Median :3.473  
##  Mean   :249.5   Mean   :4.467   Mean   :3.393  
##  3rd Qu.:373.8   3rd Qu.:5.992   3rd Qu.:4.503  
##  Max.   :498.0   Max.   :7.482   Max.   :6.863

2.1.1 linearitas

plot(happiness ~ income, data = income.data)

Karena pola grafik linear, maka kita akan menggunakan regresi linear.

2.1.2 independensi atau autokorelasi

karena uji regresi di atas linier maka sudah dipastikan kita mengetahui dependent dan independentnya. Lalu Maka kita tidak perlu melakukan cek lagi hubungan antar variabel karena variable dependentnya hanya ada satu.

2.1.3 normalitas

hist(income.data$happiness)

hasil dari histogram di atas kira-kira datanya bersebar normal atau berbentuk lonceng (bell-shaped), maka data kita merupakan normal.

2.2 Linear regression analisis

income.happiness.lm <- lm(happiness ~ income, data = income.data)

summary(income.happiness.lm)

## 
## Call:
## lm(formula = happiness ~ income, data = income.data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.02479 -0.48526  0.04078  0.45898  2.37805 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.20427    0.08884   2.299   0.0219 *  
## income       0.71383    0.01854  38.505   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.7181 on 496 degrees of freedom
## Multiple R-squared:  0.7493, Adjusted R-squared:  0.7488 
## F-statistic:  1483 on 1 and 496 DF,  p-value: < 2.2e-16

Dari hasil tersebut didapat bahwa p-value < dari alpha 5%, dan R-square sebesar 74,88%. Berarti terdapat 74,88% variabel yang sudah termsuk kedalam regresi. Regresi di atas sudah bisa kita gunakan walaupun belum terlalu baik. Kalau dilihat dari variabel dependen income memiliki pengaruh yang signifikan terhadap variabel independen happiness.
Jadi model refresinya adalah : happiness = 0.204 + 0.714*income

2.3 visualisasi Hasil dengan grafik

ggplot(income.data, aes(x=income, y=happiness))+
  geom_point()+ 
  geom_smooth(method="lm", col="black")+
  stat_regline_equation(label.x = 3, label.y = 7)+
  theme_bw() +
  labs(title = "Reported happiness as a function of income",
      x = "Income (x$10,000)",
      y = "Happiness score (0 to 10)")

## `geom_smooth()` using formula = 'y ~ x'

2.3.1 homosdesitas (kehomogenitas varians)

par(mfrow=c(1,1))
plot(income.happiness.lm)

Rdari hasil residunya, kita bisa bilang kalau model yang dibuat homogen (variansnya sama).

3 Multiple Regression

3.1 import data

heart.data = read.csv("heart.data.csv")
summary(heart.data)

##        X             biking          smoking        heart.disease    
##  Min.   :  1.0   Min.   : 1.119   Min.   : 0.5259   Min.   : 0.5519  
##  1st Qu.:125.2   1st Qu.:20.205   1st Qu.: 8.2798   1st Qu.: 6.5137  
##  Median :249.5   Median :35.824   Median :15.8146   Median :10.3853  
##  Mean   :249.5   Mean   :37.788   Mean   :15.4350   Mean   :10.1745  
##  3rd Qu.:373.8   3rd Qu.:57.853   3rd Qu.:22.5689   3rd Qu.:13.7240  
##  Max.   :498.0   Max.   :74.907   Max.   :29.9467   Max.   :20.4535

3.2 cek asumsi

Sama seperti simple regression, kita harus mengecek asumsi-asumsi yang ada.

3.2.1 linearitas

plot(heart.disease ~ biking, data=heart.data)

plot(heart.disease ~ smoking, data=heart.data)

meskipun hubungan antara smoking dan heart desease agak kurang jelas, namun tetap terlihat linier.

3.2.2 independensi atau autokorelasi

cor(heart.data$biking, heart.data$smoking)

## [1] 0.01513618

korelasi antara biking dan smoking hanya 0.015 atau 1.5%, maka kita simpulkan bahwa variable dependentnya tidak ber autokorelasi.

3.2.3 normalitas

hist(heart.data$heart.disease)

dari hasil histogram di atas histogram itu berbentuk lonceng maka bisa dinyatakan bahwa model di atas merupakan data yang normal.

3.3 Linear regression analisis

heart.disease.lm<-lm(heart.disease ~ biking + smoking, data = heart.data)

summary(heart.disease.lm)

## 
## Call:
## lm(formula = heart.disease ~ biking + smoking, data = heart.data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.1789 -0.4463  0.0362  0.4422  1.9331 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 14.984658   0.080137  186.99   <2e-16 ***
## biking      -0.200133   0.001366 -146.53   <2e-16 ***
## smoking      0.178334   0.003539   50.39   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.654 on 495 degrees of freedom
## Multiple R-squared:  0.9796, Adjusted R-squared:  0.9795 
## F-statistic: 1.19e+04 on 2 and 495 DF,  p-value: < 2.2e-16

Dari hasil tersebut didapat bahwa p-value < dari alpha 5%, dan R-square sebesar 97,96%. Berarti terdapat 97,96% variabel yang sudah termsuk kedalam regresi. Regresi di atas sudah sangat baik untuk digunakan.

3.4 visualisasi Hasil dengan grafik

plotting.data<-expand.grid(
  biking = seq(min(heart.data$biking), max(heart.data$biking), length.out=30),
    smoking=c(min(heart.data$smoking), mean(heart.data$smoking), max(heart.data$smoking)))

plotting.data$predicted.y <- predict.lm(heart.disease.lm, newdata=plotting.data)

plotting.data$smoking <- as.factor(round(plotting.data$smoking, digits = 2))


ggplot(heart.data, aes(x=biking, y=heart.disease)) +
  geom_point()+
  geom_line(data=plotting.data, aes(x=biking, y=predicted.y, color=smoking), size=1.25)+
  theme_bw() +
  labs(title = "Rates of heart disease (% of population) \n as a function of biking to work and smoking",
      x = "Biking to work (% of population)",
      y = "Heart disease (% of population)",
      color = "Smoking \n (% of population)")

## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

3.4.1 homosdesitas (kehomogenitas varians)

par(mfrow=c(1,1))
plot(heart.disease.lm)

Dalam survei kami terhadap 500 kota, kami menemukan hubungan yang signifikan antara frekuensi bersepeda ke tempat kerja dan frekuensi penyakit jantung serta frekuensi merokok dan frekuensi penyakit jantung (masing-masing p < 0 dan p < 0,001).

Secara khusus kami menemukan penurunan 0,2% (± 0,0014) pada frekuensi penyakit jantung untuk setiap peningkatan 1% dalam bersepeda, dan peningkatan 0,178% (± 0,0035) pada frekuensi penyakit jantung untuk setiap peningkatan 1% dalam merokok.