Regres Linier

library(faraway)
data(stat500)
head(stat500, n = 2)

Dapatkah kita memprediksi kinerja mahasiswa pada ujian akhir berdasarkan skor tengah semester mereka? Apakah ada korelasi antara keduanya? Dapatkah kita memprediksi kinerja mahasiswa pada ujian akhir hanya berdasarkan skor ujian akhir mereka?

final <- stat500$final
midterm <- stat500$midterm

plot(final ~ midterm, xlim = c(0, 35), ylim = c(0,
35))
abline(0, 1)

x <- mean(midterm)
sdx <- sd(midterm)
y <- mean(final)
sdy <- sd(final)
plot(final ~ midterm)
arrows(x, min(final), x, max(final), code = 0)
arrows(min(midterm), y, max(midterm), y, code = 0)

Ketika bekerja dengan dua dimensi, yang (mungkin) menggunakan unit yang sangat berbeda dan memiliki rentang yang sangat berbeda. Oleh karena itu, dilakukan standarisasi nilai akhir akhir dan nilai tengah semester (yaitu, mengubah nilai menjadi nilai-z sehingga mereka memiliki rata-rata 0 dan standar deviasi 1)

scaledstat500 <- data.frame(scale(stat500))

final.original <- final
midterm.original <- midterm
final <- scaledstat500$final
midterm <- scaledstat500$midterm

plot nilai akhir yang terstandarisasi terhadap nilai nilai ujian tengah semester

plot(final ~ midterm)
arrows(mean(midterm), min(final), mean(midterm),
max(final), code = 0)
arrows(min(midterm), mean(final), max(midterm),
mean(final), code = 0)
text(1, 2, labels = expression(x[i] %*% y[i]),
cex = 1.2)
text(1.5, 2, labels = c(" = +ve"), cex = 1.2)
text(1, -2, labels = expression(x[i] %*% -y[i]),
cex = 1.2)
text(1.55, -2, labels = c(" = -ve"), cex = 1.2)
text(-1.1, -2, labels = expression(-x[i] %*%
-y[i]), cex = 1.2)
text(-0.5, -2, labels = c(" = +ve"), cex = 1.2)
text(-1, 2, labels = expression(-x[i] %*% y[i]),
cex = 1.2)
text(-0.5, 2, labels = c(" = -ve"), cex = 1.2)

Menghitung korelasi (seperti yang didefinisikan di atas) secara eksplisit, dan kemudian menggunakan fungsi bawaan R cor() untuk memverifikasi perhitungan kita.

sum(final * midterm)/(length(final) - 1)

## [1] 0.5452277

cor(midterm, final)

## [1] 0.5452277

Korelasi bernilai positif menunjukkan bahwa sebagian besar pasangan (x,y) terletak di kuadran pertama dan ketiga. Jadi kita dapat mengatakan, secara kasar,bahwa semakin tinggi nilai ujian tengah semester, semakin tinggi pula nilai akhirnya.

(one.SD.above <- subset(scaledstat500, 0.9 <
midterm & midterm < 1.1))

mean(one.SD.above$final)

## [1] 0.6578777

 (one.SD.below <- subset(scaledstat500, -0.9 >
midterm & midterm > -1.1))

mean(one.SD.below$final)

## [1] 0.304624

Jadi, ketika nilai ujian tengah semester 1 SD di atas rata-rata ujian tengah semester, nilai akhir hanya 0,65 SD di atas rata-rata nilai akhir. Ketika nilai ujian tengah semester 1 SD di bawah nilai rata-rata ujian tengah semester, maka nilai akhir hanya 0,30 SD di bawah nilai rata-rata ujian akhir.

di bawah rata-rata nilai akhir. Nilai akhir mengalami kemunduran terhadap nilai rata-rata. Nilai-nilai tersebut tidak akan mengalami kemunduran terhadap rata-rata jika rasio nilai akhir dan tengah semester adalah 1:1. Misalnya, Untuk membuat poin ini lebih konkret, mari kita pisahkan perbedaan antara kedua hasil ini dan gambarkan garis y = 0,475 × x, dan bandingkan dengan garis di mana perubahan 1 SD pada nilai ujian tengah semester menghasilkan perubahan 1 SD pada nilai akhir.

plot(final ~ midterm)
arrows(mean(midterm), min(final), mean(midterm),
max(final), code = 0)
arrows(min(midterm), mean(final), max(midterm),
mean(final), code = 0)
abline(0, 1, col = "black")
abline(0, 0.5452, col = "black", lty = 2)
text(1.5, 2, labels = c("1:1 ratio of change"),
col = "black")
text(1.45, 0.3, labels = c("0.5:1 ratio of change"),
col = "black")

summary(final~midterm)

##  Length   Class    Mode 
##       3 formula    call

Untuk mendapatkan kemiringan pada data yang tidak terstandardisasi, kita dapat memasukkan garis kuadrat terkecil pada data asli

finall <- final.original
midtermm <- midterm.original
lm.stat500 <- lm(finall ~ midtermm)
summary(lm.stat500)

## 
## Call:
## lm(formula = finall ~ midtermm)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -9.932 -2.657  0.527  2.984  9.286 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  15.0462     2.4822   6.062 1.44e-07 ***
## midtermm      0.5633     0.1190   4.735 1.67e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.192 on 53 degrees of freedom
## Multiple R-squared:  0.2973, Adjusted R-squared:  0.284 
## F-statistic: 22.42 on 1 and 53 DF,  p-value: 1.675e-05

plot(finall ~ midtermm)
abline(lm.stat500)
text(15, 24, labels = c("y = 15.0462 + 0.5633x"),
     cex = 1.2)

Regresi dari ujian tengah semester akhir akan memberikan kemiringan yang berbeda. Karena kita menggunakan regresi untuk memprediksi nilai y, maka yang menjadi perhatian adalah kesalahan dalam arah vertikal. Intersep adalah nilai akhir yang diprediksi ketika nilai ujian tengah semester adalah 0, dan kemiringan adalah peningkatan nilai akhir ketika nilai ujian tengah semester meningkat sebesar 1.

Uji Hipotesis Regresi

Dalam regresi, kita meminimalkan jumlah kuadrat residual ketika kita menyesuaikan garis regresi- itulah yang dilakukan oleh metode kuadrat terkecil! Residual ini dihitung oleh fungsi lm().

residuals <- residuals(lm(final ~ midterm))
tail(residuals, n = 5)

##          51          52          53          54          55 
## -0.14531326  0.60652517 -0.95275038  0.61103277 -0.05715639

sum(residuals^2)

## [1] 37.94724

ms.residuals <- sum(residuals^2)/53
ms.residuals

## [1] 0.7159857

Untuk setiap residual, kita memiliki nilai pada garis regresi yang sesuai dengan nilai rata-rata yang sesuai.

fitted.values <- fitted.values(lm(final ~ midterm))
tail(fitted.values)

##         50         51         52         53         54         55 
## -0.2635829  0.2480780 -0.6046901  0.2480780  0.7028876  0.3617804

mean(final)

## [1] -2.157444e-16

ms.fitted <- sum((fitted.values - mean(final))^2)/1

f.ratios <- rep(NA, 1000)
for (i in 1:1000) {
final.rand <- sample(final)
lm <- lm(final.rand ~ midterm)
residuals <- residuals(lm)
ms.residuals <- sum(residuals^2)/53
fitted.values <- fitted.values(lm)
ms.fitted <- sum((fitted.values - mean(final))^2)/1
f.ratios[i] <- ms.fitted/ms.residuals
}

main.title <- "Distribution of \n simulated F-ratios"
plot(density(f.ratios), xlim = range(0, 8), main = main.title,
xlab = "", ylab = "")

main.title <- "Distribution of \n theoretical F-ratios"
plot(function(x) df(x, 1, 53), -0.1, 8, main = main.title,
xlab = "", ylab = "")