3 R ile Basit Doğrusal Regresyon Uygulaması

Uzun uzun işlem yapmak istemeyen, sadece sonucu görmek isteyenler, aşağıdaki R’da Build-in Fonksiyonlar Yardımıyla Basit Doğrusal Regresyon kısmına geçebilir.

3.0.1 Veri Setine Genel Bakış

3.0.1.1 Veri setinin import edilmesi.

library(readr)
Advertising <- read_csv("Advertising.csv")
df <- Advertising <- Advertising[,2:5]
str(df)
## Classes 'tbl_df', 'tbl' and 'data.frame':    200 obs. of  4 variables:
##  $ TV       : num  230.1 44.5 17.2 151.5 180.8 ...
##  $ radio    : num  37.8 39.3 45.9 41.3 10.8 48.9 32.8 19.6 2.1 2.6 ...
##  $ newspaper: num  69.2 45.1 69.3 58.5 58.4 75 23.5 11.6 1 21.2 ...
##  $ sales    : num  22.1 10.4 9.3 18.5 12.9 7.2 11.8 13.2 4.8 10.6 ...
summary(df)
##        TV             radio          newspaper          sales      
##  Min.   :  0.70   Min.   : 0.000   Min.   :  0.30   Min.   : 1.60  
##  1st Qu.: 74.38   1st Qu.: 9.975   1st Qu.: 12.75   1st Qu.:10.38  
##  Median :149.75   Median :22.900   Median : 25.75   Median :12.90  
##  Mean   :147.04   Mean   :23.264   Mean   : 30.55   Mean   :14.02  
##  3rd Qu.:218.82   3rd Qu.:36.525   3rd Qu.: 45.10   3rd Qu.:17.40  
##  Max.   :296.40   Max.   :49.600   Max.   :114.00   Max.   :27.00

3.0.2 Korelasyon Matrisi

Değişkenlerin birbirleri üzerindeki etkisini incelemek için korelasyon matrisi kullanılabilir.

#library("PerformanceAnalytics")
chart.Correlation(df, histogram=TRUE, pch=19)
Tüm Değişkenler İçin Korelasyon Matrisi

Şekil 3.1: Tüm Değişkenler İçin Korelasyon Matrisi

Şekil 3.1’de TV ile sales değişkenleri arasında güçlü doğrusal bir ilişki gözlemlenmektedir.

Sadece TV ve sales değişkenlerinden oluşan yeni bir data frame oluşturalım.

df1 <- data.frame(sales = df$sales,TV = df$TV)
library(knitr)
kable(head(df1,6)) # Yeni oluşturduğumuz veri setinin ilk 6 satırı.
sales TV
22.1 230.1
10.4 44.5
9.3 17.2
18.5 151.5
12.9 180.8
7.2 8.7
chart.Correlation(df1, histogram=F, pch=10)
sales~TV Değişkenleri Arasındaki Korelasyon Matrisi

Şekil 3.2: sales~TV Değişkenleri Arasındaki Korelasyon Matrisi

Şekil 3.2’de görüldüğü gibi sales ve TV değişkenleri arasında doğrusal bir ilişki vardır.

Grafikte gözükse de cor() fonksiyonu yardımıyla sales ve TV değişkenleri arasındaki korelasyona ulaşılabilir.

cor(df1$sales,df1$TV)
## [1] 0.7822244

TV’ye verilen reklam ile Satışlar arasında güçlü bir ilişki olduğunu görüyoruz. Bu değişkenler ile regresyon modeli kurabiliriz. Reklamların satışlar üzerindeki etkisini araştırıdğımız için X = TV(bağımsız değişken), Y = Sales(bağımlı değişken) olarak belirleyebiliriz.

3.1 Model Katsayılarının Tahmini

ˆβ1=ni=1(xiˉx)(yiˉy)ni=1(xiˉx)2

x <- df1$TV
y <- df1$sales
b_1 <- sum((x-mean(x))*(y-mean(y)))/(sum((x-mean(x))^2))
b_1
## [1] 0.04753664

ˆβ0=ˉyˆβ1ˉx

b_0 <- mean(y) - b_1*mean(x)
b_0
## [1] 7.032594

Bulduğumuz katsayılar ile aşağıdaki eşitlik yazılabilir.

 sales 7.032594+0.04753664×TV

Yukarıdaki eşitliği yorumlayacak olursak. TV değişkenindeki 1 birimlik değişiklik sales değişkeninde 0.04753664 birimlik değişikliğe neden olur.

Bulduğumuz katsayılar ile regresyon doğrusunu çizdirecek olursa;

plot(x,y,pch=3,lwd = 1,xlab="x (TV)",ylab = "y (Sales)"); abline(b_0,b_1,col="red")
Regresyon Doğrusu

Şekil 3.3: Regresyon Doğrusu

3.2 Katsayıların ve Artıkların Varyansı

3.2.1 Artıkların Varyansı

RSE=RSS/(n2)

RSS

y_sapka <- b_0 + b_1*x # y şapka değerleri
e <- y - y_sapka # artıklar
rss <- sum(e^2)
rss
## [1] 2102.531

RSE

#residual standard error
rse <- sqrt(rss/(length(x)-2))
rse
## [1] 3.258656

3.2.2 Katsayıların Varyansları

SE(ˆβ0)2=σ2[1n+ˉx2ni=1(xiˉx)2],SE(ˆβ1)2=σ2ni=1(xiˉx)2

β0’ın standart hatası;

varbeta0 <- rse^2*(1/length(x) + (mean(x)^2)/ sum((x-mean(x))^2)) # beta 0'ın varyansı
se_beta0 <- sqrt(varbeta0)
se_beta0
## [1] 0.4578429

β1’in standart hatası;

varbeta1 <- rse^2/sum((x-mean(x))^2) # beta 1'in varyansı
se_beta1 <- sqrt(varbeta1)
se_beta1
## [1] 0.002690607

3.3 Katsayılar İçin Güven Aralığı

ˆβ1 için %95 güven aralığı; ˆβ1±2SE(ˆβ1)

beta1_alt_limit <- b_1 - 2 * se_beta1
beta1_ust_limit <- b_1 + 2 * se_beta1
beta1_ga <- c(beta1_alt_limit,beta1_ust_limit)
names(beta1_ga) <- c("Alt Limit","Üst Limit")
beta1_ga
##  Alt Limit  Üst Limit 
## 0.04215543 0.05291785

ˆβ0 için %95 güven aralığı; ˆβ0±2SE(ˆβ0)

beta0_alt_limit <- b_0 - 2 * se_beta0
beta0_ust_limit <- b_0 + 2 * se_beta0
beta0_ga <- c(beta0_alt_limit,beta0_ust_limit)
names(beta0_ga) <- c("Alt Limit","Üst Limit")
beta0_ga
## Alt Limit Üst Limit 
##  6.116908  7.948279

3.4 Katsayıların ve Modelin Anlamlılığı

3.4.1 Katsayıların Anlamlılığı

Katsayıların anlamlılığını ölçmek için önce katsayıların t değerlerini bulmalıyız.

t=ˆβ10SE(ˆβ1)

t_beta1 <- b_1 / se_beta1
t_beta1
## [1] 17.66763
t_beta0 <- b_0 / se_beta0
t_beta0
## [1] 15.36028

3.4.2 t tablo değerinin bulunması

qt(0.025,length(x)-2); qt(0.025,length(x)-2,lower.tail = FALSE)
## [1] -1.972017
## [1] 1.972017

3.4.3 Hipotezler

H0:β1=0 HA:β10

Çift taraflı yapılan hipotez testinde ˆβ1 ve ˆβ0 mutlak değerce 0.05 anlamlılık düzeyinde tablo değerimizden büyük olduğu için H0 hipotezi reddedilir. Katsayılarımız anlamlıdır.

Ayrıca hesapladığımız ˆβ0 ve ˆβ1 değelerleri, hesapladığımız %95 güven aralıklarının içinde olduğu için de anlamlı olduklarını söyleyebiliridk.

beta1_ga; paste("Beta1:", b_1); beta0_ga; paste("Beta0:", b_0)
##  Alt Limit  Üst Limit 
## 0.04215543 0.05291785
## [1] "Beta1: 0.0475366404330197"
## Alt Limit Üst Limit 
##  6.116908  7.948279
## [1] "Beta0: 7.03259354912769"

3.5 R2 Modelin Açıklayıcılığı

Modelin genel olarak anlamlığını incelemek içi R2 değerinden yararlanılabilir.

R2=TSSRSSTSS=1RSSTSS

tss <- sum((y-mean(y))^2) # Toplam Kareler Toplamı
r_2 <- 1 - rss/tss
r_2
## [1] 0.6118751

Aynı değeri X ve Y değişkenlerinin aralarındaki korelasyonun karesiyle de bulabiliriz.

Cor(X,Y)=ni=1(xiˉx)(yiˉy)ni=1(xiˉx)2ni=1(yiˉy)2

cor(x,y)^2
## [1] 0.6118751

Katsayılarımız anlamlı ancak modelimizin çok iyi sonuçlar vereceğini söyleyemeyiz.