Gerekli Kütüphaneler

# Paketlerin yüklü olduğundan emin olunuz, kurulu değilse console'dan install.packages("paketadi") ile kurunuz.
library(ggplot2)
library(dplyr)
library(broom)
library(ggpubr)
library(ISLR)
library(caret)
library(mice)
library(PerformanceAnalytics)

Basit Doğrusal Regresyon - İlk Örnek

Aylık Gelir ve Kültürel Harcama

x <- c(850,1000,1400,1750,2000,2200,2500,3400)
y <- c(80,90,175,225,250,275,300,345)
dataz <- data.frame(x,y)
dataz
##      x   y
## 1  850  80
## 2 1000  90
## 3 1400 175
## 4 1750 225
## 5 2000 250
## 6 2200 275
## 7 2500 300
## 8 3400 345
plot(x,y,xlab = "aylık gelir",ylab = "Kulturel harcama")
model <- lm(y~x)
summary(model)
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -39.55 -24.55  13.10  20.73  22.98 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  9.03129   25.99569   0.347 0.740142    
## x            0.11045    0.01272   8.681 0.000129 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 28.15 on 6 degrees of freedom
## Multiple R-squared:  0.9262, Adjusted R-squared:  0.914 
## F-statistic: 75.35 on 1 and 6 DF,  p-value: 0.000129
abline(model,col = "red")

predictd <- predict(model)
datazz <- data.frame(dataz,"predictions" = predictd)
datazz
##      x   y predictions
## 1  850  80    102.9112
## 2 1000  90    119.4783
## 3 1400 175    163.6571
## 4 1750 225    202.3135
## 5 2000 250    229.9253
## 6 2200 275    252.0147
## 7 2500 300    285.1488
## 8 3400 345    384.5511
error <- y-predictd
error
##         1         2         3         4         5         6         7         8 
## -22.91124 -29.47829  11.34291  22.68646  20.07471  22.98531  14.85121 -39.55109
datazzz <- data.frame(datazz,"error" = error)
datazzz
##      x   y predictions     error
## 1  850  80    102.9112 -22.91124
## 2 1000  90    119.4783 -29.47829
## 3 1400 175    163.6571  11.34291
## 4 1750 225    202.3135  22.68646
## 5 2000 250    229.9253  20.07471
## 6 2200 275    252.0147  22.98531
## 7 2500 300    285.1488  14.85121
## 8 3400 345    384.5511 -39.55109
sum(error)
## [1] 8.526513e-14

İkinci Örnek - İticinin Yaşı ve Kesme Değeri

#defterdeki ikinci örnek soru;
x <- c(15.5,23.75,8,17,5.5,19,24,2.5,7.5,11,13,3.75,25,9.75,22,18,6,12.5,2,21.5)
y <- c(2158.70,1678.15,2316,2061.30,2207.5,1708.3,1784.7,2575,2357.9,2256.7,2165.2,2399.55,1779.8,2336.75,1765.3,2053.5,2414.4,2200.5,2654.2,1753.7)
plot(x,y,xlab = "İticinin Yaşı", ylab = "Kesme Değeri")
model2 <- lm(y~x)
summary(model2)
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -215.98  -50.68   28.74   66.61  106.76 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 2627.822     44.184   59.48  < 2e-16 ***
## x            -37.154      2.889  -12.86 1.64e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 96.11 on 18 degrees of freedom
## Multiple R-squared:  0.9018, Adjusted R-squared:  0.8964 
## F-statistic: 165.4 on 1 and 18 DF,  p-value: 1.643e-10
abline(model2,col = "red")

Manuel Regresyon Hesaplamaları

Teslim Hacmi ve Teslim Süresi Örneği

#Örnek;
# x teslim hacmi, y teslim süresi
x <- c(7,3,3,4,6,7,2,7,30,5,16,10,4,6,9,10,6,7,3,17,10,26,9,8,4)
y <- c(16.68,11.5,12.03,14.88,13.75,18.11,8,17.83,79.24,21.5,40.33,21,13.5,19.75,24,29,15.35,19,9.5,35.10,17.9,52.3,18.75,19.83,10.75)
plot(x,y)

# B1 şapka ve B0 şapka bulalım Sxy ve Sxx bulunmalı. Formüller defterde var.

sxx <- sum((x-mean(x))^2)
sxx
## [1] 1136.56
sxy <- sum((x-mean(x))*y)
sxy
## [1] 2472.999
beta1sapka <- sxy/sxx
# B1sapka = sum(i = 1' den n'e kadar, ci*yi)    ci = (xi - xbar) / sxx
beta1sapka
## [1] 2.175863
beta0sapka <- mean(y) - beta1sapka*mean(x)
# B0sapka = sum(i = 1' den n'e kadar, di*yi)    di = (1/n) - ci*xbar

beta0sapka
## [1] 3.322637
# ysapka = 3.32 + 2.1761x çıkmış oldu.
# Bunun hazırı da var ve onu kullanacağız normalde fakat bu da önemli

model5 <- lm(y~x)
summary(model5)
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -7.5951 -1.8737 -0.3502  2.1798 10.6415 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   3.3226     1.3716   2.422   0.0237 *  
## x             2.1759     0.1241  17.536 8.31e-15 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.183 on 23 degrees of freedom
## Multiple R-squared:  0.9304, Adjusted R-squared:  0.9274 
## F-statistic: 307.5 on 1 and 23 DF,  p-value: 8.312e-15
#Hacmi sıfırken teslim süresi 3.32 imiş. Bu ikisi arası ilişki pozitif yönlü. Malın 
#hacmindeki 1 birimlik artış 2.17 birimlik saat bazı bir artışa sebep olur.

# çıkan tabloda Pr(>|t|) değeri ile alpha değeri kıyaslanır. 0.05 alpha değeri.
# 0.05' den kçüük olduğu için H0 reddedilir. Yani B0 anlamlı. B0 modelimde olmalı.

#B1 ise intercept değerinin altındaki x kısmından bakılır.
# onda ise 0' çok yakın bir değer o yüzden H0 yine reddedilir. B1 anlamlı yani. 
# Teslim hacmi ve teslim süresi anlamlı ilişkiye sahip. Bu model anlamlı bir model çkıtı.


#4.183 olarak çıkan sigma şapka. Yani standart hata denen şey sigma şapkadır.

tahmin <- predict(model5)
xx <- data.frame("Gercek" = y, "Tahmin" = tahmin)
xx
##    Gercek    Tahmin
## 1   16.68 18.553681
## 2   11.50  9.850227
## 3   12.03  9.850227
## 4   14.88 12.026091
## 5   13.75 16.377817
## 6   18.11 18.553681
## 7    8.00  7.674364
## 8   17.83 18.553681
## 9   79.24 68.598537
## 10  21.50 14.201954
## 11  40.33 38.136450
## 12  21.00 25.081271
## 13  13.50 12.026091
## 14  19.75 16.377817
## 15  24.00 22.905407
## 16  29.00 25.081271
## 17  15.35 16.377817
## 18  19.00 18.553681
## 19   9.50  9.850227
## 20  35.10 40.312314
## 21  17.90 25.081271
## 22  52.30 59.895083
## 23  18.75 22.905407
## 24  19.83 20.729544
## 25  10.75 12.026091
#Şimdi aradaki farklar ile hatalrı bulalım

error <- y-tahmin
data.frame(xx,"Artik" = error)
##    Gercek    Tahmin      Artik
## 1   16.68 18.553681 -1.8736806
## 2   11.50  9.850227  1.6497726
## 3   12.03  9.850227  2.1797726
## 4   14.88 12.026091  2.8539093
## 5   13.75 16.377817 -2.6278173
## 6   18.11 18.553681 -0.4436806
## 7    8.00  7.674364  0.3256360
## 8   17.83 18.553681 -0.7236806
## 9   79.24 68.598537 10.6414634
## 10  21.50 14.201954  7.2980460
## 11  40.33 38.136450  2.1935497
## 12  21.00 25.081271 -4.0812705
## 13  13.50 12.026091  1.4739093
## 14  19.75 16.377817  3.3721827
## 15  24.00 22.905407  1.0945928
## 16  29.00 25.081271  3.9187295
## 17  15.35 16.377817 -1.0278173
## 18  19.00 18.553681  0.4463194
## 19   9.50  9.850227 -0.3502274
## 20  35.10 40.312314 -5.2123136
## 21  17.90 25.081271 -7.1812705
## 22  52.30 59.895083 -7.5950834
## 23  18.75 22.905407 -4.1554072
## 24  19.83 20.729544 -0.8995439
## 25  10.75 12.026091 -1.2760907
sum(error)
## [1] 1.882938e-13
#Sıfır kabul edilebilecek bir değer geldi.
sum(y)
## [1] 559.58
sum(tahmin)
## [1] 559.58
#Aynı çıktılar.


#SSE hesaplarsak (Formül defterde var.)

sse <- sum((y-tahmin)^2)
sse
## [1] 402.4373
# sgima^2 şapka yani hata kareler ortalaması için ise

hko <- sse/(length(x)-2)
hko
## [1] 17.49728
# Bunun karekökü standart error olur. Tabloda olan.
sigmasapka <- sqrt(hko)
sigmasapka
## [1] 4.182974
#Ayrı ayrı betaların standart hataları ise;
sebeta1 <- sqrt(hko/sxx)
sebeta0 <- sqrt(hko * ((1/length(y)) + (mean(x)^2) / sum(x - mean(x))^2))

#B1 t değeri ise;
thesapbeta1 <- beta1sapka/ sebeta1
thesapbeta1
## [1] 17.53649
#Yani kısaca bu değerler hep veriliyor tabloda fakat biz de kendimiz bulabiliyoruz.

#F değeri ile ise model anlamlılığı test edilir. B1 t değerinin karesi alınırsa 
#F değeri çıkar. 
#(Alpha değeri ile de zaten anlamsızlığa bakabiliyorduk.)

Güven Aralıkları

# Teslim süresi teslim hacmi örneği için; sxx sxy beta1sapka beta0sapka bulmuştuk.
# B0 ve B1 için güven aralıklarını bulacağız. (Formülleri defterde var.)
# B0 şapka ve B1 şapka için ilk önce onların standart hatalarını bulamlıyız.

dataz7 <- data.frame("actual" = y, "tahmin" = tahmin)
dataz7
##    actual    tahmin
## 1   16.68 18.553681
## 2   11.50  9.850227
## 3   12.03  9.850227
## 4   14.88 12.026091
## 5   13.75 16.377817
## 6   18.11 18.553681
## 7    8.00  7.674364
## 8   17.83 18.553681
## 9   79.24 68.598537
## 10  21.50 14.201954
## 11  40.33 38.136450
## 12  21.00 25.081271
## 13  13.50 12.026091
## 14  19.75 16.377817
## 15  24.00 22.905407
## 16  29.00 25.081271
## 17  15.35 16.377817
## 18  19.00 18.553681
## 19   9.50  9.850227
## 20  35.10 40.312314
## 21  17.90 25.081271
## 22  52.30 59.895083
## 23  18.75 22.905407
## 24  19.83 20.729544
## 25  10.75 12.026091
error <- y-tahmin
data.frame(dataz7,"Artik" = error)
##    actual    tahmin      Artik
## 1   16.68 18.553681 -1.8736806
## 2   11.50  9.850227  1.6497726
## 3   12.03  9.850227  2.1797726
## 4   14.88 12.026091  2.8539093
## 5   13.75 16.377817 -2.6278173
## 6   18.11 18.553681 -0.4436806
## 7    8.00  7.674364  0.3256360
## 8   17.83 18.553681 -0.7236806
## 9   79.24 68.598537 10.6414634
## 10  21.50 14.201954  7.2980460
## 11  40.33 38.136450  2.1935497
## 12  21.00 25.081271 -4.0812705
## 13  13.50 12.026091  1.4739093
## 14  19.75 16.377817  3.3721827
## 15  24.00 22.905407  1.0945928
## 16  29.00 25.081271  3.9187295
## 17  15.35 16.377817 -1.0278173
## 18  19.00 18.553681  0.4463194
## 19   9.50  9.850227 -0.3502274
## 20  35.10 40.312314 -5.2123136
## 21  17.90 25.081271 -7.1812705
## 22  52.30 59.895083 -7.5950834
## 23  18.75 22.905407 -4.1554072
## 24  19.83 20.729544 -0.8995439
## 25  10.75 12.026091 -1.2760907
# tahmin ve artıkların olduğu dataframi tekrar yazdık.
#SSE hata kareler toplamını bulalım.
SSE <- sum((error)^2)
SSE
## [1] 402.4373
sigmakaresapka <- SSE/(length(y) - 2)
sigmakaresapka
## [1] 17.49728
sigmasapka <- sqrt(sigmakaresapka)
sigmasapka
## [1] 4.182974
#B1 standart hatası
sebeta1 <- sqrt(sigmakaresapka/sxx)
#B0 standart hatası
sebeta0 <- sqrt((sigmakaresapka*((1/length(y))+(mean(x)^2/sxx))))

#B0 için güven aralığı
beta0interval <- c(beta0sapka - (qt(1-0.025,23)*sebeta0), beta0sapka + (qt(1-0.025,23)*sebeta0))

#B1 için güven aralığı                    
beta1interval <- c(beta1sapka - (qt(1-0.025,23)*sebeta1), beta1sapka + (qt(1-0.025,23)*sebeta1))

beta0interval; beta1interval                    
## [1] 0.4852853 6.1599896
## [1] 1.919192 2.432535
#El ile yaptığımız bu her şeyi paket ile yapabiliriz;
confint(model5, level = 0.95)
##                 2.5 %   97.5 %
## (Intercept) 0.4852853 6.159990
## x           1.9191918 2.432535

Varyans Analizi Tablosu

# Varyans analizi tablosu
# SST bulalım
SST <- sum((y-mean(y))^2)
#SSR bullaım
SSR <- beta1sapka*sxy
SSR <- sum((tahmin - mean(y))^2)# Bu şekilde de bulunabilir.
# SSR + SSE = SST

Anova_table <- rbind(c(1,SSR,SSR/1,(SSR/1)/(SSE/(length(y)-2))),
                     c((length(y)-2),SSE,SSE/(length(y)-2),(SSR/1)/(SSE/(length(y)-2))),
                     c((length(y)-1),SST,SST/(length(y)-1),(SSR/1)/(SSE/(length(y)-2))))
colnames(Anova_table) <- c("SD","Kareler Toplamı","Kareler ORT","Fhesap")
rownames(Anova_table) <- c("Regresyon","Hata","Genel")
Anova_table
##           SD Kareler Toplamı Kareler ORT   Fhesap
## Regresyon  1       5380.9082  5380.90822 307.5284
## Hata      23        402.4373    17.49728 307.5284
## Genel     24       5783.3455   240.97273 307.5284
#El işle yaptığımız tabloyu direk böyle de yapabiliyoruz.
anova(model5)
## Analysis of Variance Table
## 
## Response: y
##           Df Sum Sq Mean Sq F value    Pr(>F)    
## x          1 5380.9  5380.9  307.53 8.312e-15 ***
## Residuals 23  402.4    17.5                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Rkare <- SSR/SST

summary(model5)
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -7.5951 -1.8737 -0.3502  2.1798 10.6415 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   3.3226     1.3716   2.422   0.0237 *  
## x             2.1759     0.1241  17.536 8.31e-15 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.183 on 23 degrees of freedom
## Multiple R-squared:  0.9304, Adjusted R-squared:  0.9274 
## F-statistic: 307.5 on 1 and 23 DF,  p-value: 8.312e-15

Çoklu Doğrusal Regresyon

Matrislerle Çözüm

# Matrislerle Çözüm

y <- c(6,6,11,9,16,16,4,8,11,13,13,9,17,17,12,6,5,12,8,9)
x1 <- c(1,2,2,1,3,1,1,3,2,3,1,1,3,2,4,1,1,3,1,2)
x2 <- c(1,1,2,3,3,5,1,1,2,2,4,2,3,4,1,1,1,2,2,2)

model1 <- lm(y~x1+x2)
summary(model1)
## 
## Call:
## lm(formula = y ~ x1 + x2)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.7443 -0.7314  0.1787  0.8277  1.4511 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -0.5040     0.7419  -0.679    0.506    
## x1            2.4023     0.2558   9.393 3.84e-08 ***
## x2            2.9486     0.2094  14.081 8.41e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.071 on 17 degrees of freedom
## Multiple R-squared:  0.9381, Adjusted R-squared:  0.9308 
## F-statistic: 128.8 on 2 and 17 DF,  p-value: 5.365e-11
zzzz <- data.frame(y,x1,x2)
model2 <- lm(y~.,data = zzzz)

# Model analmlılığına bakalım.
# p değeri çok çok küçük alpha > p yani B0 reddedilir.
# Pr(>|t|) değerlerine bakarsak B0 reddedilemiyor. Yani modelden çıkarmalıyız.
# anova tablosu yaparsak

anova(model1)
## Analysis of Variance Table
## 
## Response: y
##           Df  Sum Sq Mean Sq F value    Pr(>F)    
## x1         1  68.036  68.036  59.352 6.075e-07 ***
## x2         1 227.277 227.277 198.269 8.409e-11 ***
## Residuals 17  19.487   1.146                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#B0'ı çıkaralım
model3 <- lm(y~x1+x2-1)
summary(model3)
## 
## Call:
## lm(formula = y ~ x1 + x2 - 1)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.82770 -0.79039  0.03524  0.77915  1.62139 
## 
## Coefficients:
##    Estimate Std. Error t value Pr(>|t|)    
## x1   2.2755     0.1721   13.22 1.04e-10 ***
## x2   2.8507     0.1496   19.05 2.23e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.055 on 18 degrees of freedom
## Multiple R-squared:  0.9919, Adjusted R-squared:  0.991 
## F-statistic:  1105 on 2 and 18 DF,  p-value: < 2.2e-16
# !!Anlamlı olan modeli yorumlamalıyız.!!
# Eğer b2 anlamsız olsaydı direkt x2 yi almazdık.
# Artık anlamlı oldu. Adjsuted R squared yükseldi ve hata azaldı.

#Güven aralıkları
confint(model3)
##       2.5 %   97.5 %
## x1 1.913934 2.636973
## x2 2.536352 3.165147

Gerçek Veri Seti Uygulaması - KC House Data

Veri Hazırlama ve Eğitim-Test Ayrımı

#cArtık büyük data çekip onun üzerinde çalışabiliriz.
# Sağ üstteki importtan da yapabiliriz.
#kc house data verisetini yükledik
names(kc_house_data)

# Bu iki değişken üzerine model kuracağız
modeldata <- kc_house_data[c("price","sqft_living")]
head(modeldata)

#Eğitim test parçalaması yapmalyız.

# Datanın istediğimiz yüzdelikte kısmını seçtik
set.seed((145))
sampleind <- sample(1:nrow(modeldata),size = 0.8*nrow(modeldata))

trainset <- modeldata[sampleind,]
testset <- modeldata[-sampleind,]
nrow(trainset);nrow(testset)

cor(trainset)
# x ve y arası ilişki %70 çıktı.
plot(trainset$price,trainset$sqft_living)
# Burada aykırı değerler de olabilir ama etkili değer de olabilir.
# Aykırı var mı anlamalıyız.
# Eğer aykırılar atıldıktan sonra model performansı düşerse bu değerler
# aykırı değildir etkili edğerlerdir.

Aykırı Değer Tespiti

distance <- mahalanobis(trainset,center = colMeans(trainset),cov = cov(trainset))
cutoff <- qchisq(p = 0.95,df = 2)
ids <- which(distance > cutoff)

trainsetrem <- trainset[-ids,]
nrow(trainset);nrow(trainsetrem)
plot(trainsetrem$price,trainsetrem$sqft_living) # Veri daha iyi hale geldi
sum(is.na(trainsetrem))#Missing değer yok
cor(trainsetrem)
#Eksik veri varsa onlar ıdoldurmak için bu paket kullanılıyor.(indirilmedi indir)
# install.packages("mice")
library(mice)
md.pattern(trainset)
# Zaten temiz bu.

Model Kurma ve Karşılaştırma

# Artık modeli kurabiliriz.
model1 <- lm(price~sqft_living, data = trainset)
model2 <- lm(price~sqft_living, data = trainsetrem)
summary(model1)
# Model 1 aslında anlamlı gözüküyor ama standart hatamız çok yüksek ve 
# adjusted E squared düşük bir değer. Bu yeterli bir model değil.
summary(model2)
# Standart hata azaldı fakat r squared de azaldı. Etkili değerleri de çıkarmış
# olabiliriz. Burada yeni değişkenler eklememiz lazım. (Burayı bu ders yapmıyoruz.)
#Güven aralıkları
confint(model1) # MOdel1 güven aralığına baktık. Aldığımız atl ve üst değerlere
# bakarsak sıfırı içine almıyor. Bu nedenle anlamlı olur. B0 ve B1.
confint(model2) # Bu da anlamlı.

AIC(model1,k = 3)# Normalde 2 parametremiz var ama varyansı da alıyor. O yüzden
# 3 parametre giriyoruz.
AIC(model2,k=3) # Bu daha düşük çıktı yani daha iyi.
BIC(model1)
BIC(model2) # 2. model burada da daha iyi bir performansa sahip olduğunu görüyoruz.
# Fakat sadece bu ikisi yeterli diil.

Model Performans Değerlendirmesi

# Tahminler ( Testteki performansına da bakalım.)
model1pred <- predict(model1,testset)
model2pred <- predict(model2,testset)
model1preddata <- data.frame("actual" = testset$price, "tahmin" = model1pred)
model2preddata <- data.frame("actual" = testset$price, "tahmin" = model2pred)
model1hata <- model1preddata$actual - model1preddata$tahmin
model2hata <- model2preddata$actual - model2preddata$tahmin
mse1 <- sum(model1hata^2) / (nrow(model1preddata))
mse2 <- sum(model2hata^2) / (nrow(model2preddata))
mse1;mse2 # hata kareler ortalaması bunlar.
# RMSE <- Hata kareler ortalamasının kökü
# MAE <- Hata ortalamasının murlak değeri
#Burada gördük ki ilk model daha iyi gelmiş oldu.

rmse1 <- sqrt(mse1)
rmse2 <- sqrt(mse2)
rmse1;rmse2
# Burada da ilk model daha anlamlı çıktı yani biz aykırı değerleri çıkarırken
# çıkarmamalıydık

# install.packages("caret")
library(caret)
#Bu yaptığımız şeyler bu pakaetteki fonksiyonlarla da oluyo
RMSE(model1preddata$tahmin, model1preddata$actual)
RMSE(model2preddata$tahmin, model2preddata$actual)

MAE(model1pred,testset$price)
MAE(model2pred,testset$price) # bu ise aykırıdan arındırdığımız durum daha iyi çıktı

Çok Değişkenli Model

# Bu model değişken yetersizliğinden dolayı pek de yeterli olmuyor.
# Başka değişkenler de ekleyerek yeni model yapalım

modeldata <- kc_house_data[c("price","sqft_living")]

new_data <- kc_house_data[c("price","sqft_living","grade","bathrooms","sqft_living15")]
cor(new_data)


set.seed((145))
sampleind <- sample(1:nrow(new_data),size = 0.8*nrow(new_data))
nrow(new_data)
trainset <- new_data[sampleind,]
testset <- new_data[-sampleind,]
nrow(trainset);nrow(testset)
cor(trainset)

# install.packages("PerformanceAnalytics")
library(PerformanceAnalytics)

chart.Correlation(trainset,histogram = TRUE, pch = 19)

#Aykırı kontrolü

distance <- mahalanobis(trainset,center = colMeans(trainset),cov = cov(trainset))
cutoff <- qchisq(p = 0.95,df = 5)
ids <- which(distance > cutoff)

trainsetrem <- trainset[-ids,]
nrow(trainset);nrow(trainsetrem)

#Missing Kontrolü

sum(is.na(trainsetrem)) # Missing değer yok.
md.pattern(trainsetrem)

# Model Kurma
model11 <- lm(price ~., data = trainset)
model22 <- lm(price~., data = trainsetrem)

summary(model11)
summary(model22) # yine çok yeterli bir model çıkmadı. Yeni değişkenler eklemek de 
# bizi iyi bir modele çıkarmadı
# Ayrıca modelde anlamsız bir değer çıkmadı. Çıksaydı eğer modelden tek tek
# çıkarmak ve ayrı ayrı model kurmak gerekirdi. P değeri daha büyük olan
# modelden çıkarıp sonra diğeri hala anlamsız oluyorsa onu da çıkarmak lazım.
# Diğer kısımları kendin devam et.

Görselleştirme

# Burada ilk hali için grafikler yapıyoruz. Yaptığımız şeyin devamı diil.
par(mfrow = c(1,2))
hist(trainset$price)
hist(trainset$sqft_living)
dev.off # par komutunu ortadan kaldırdı, 
plot(trainset$price,trainset$sqft_living) # iki sayısal değişken arası scatter plot yaptık


# install.packages("PerformanceAnalytics")
library(PerformanceAnalytics)

chart.Correlation(trainset,histogram = TRUE, pch = 19)

MTCARS Veri Seti Analizi

#--------------------------------------------------------------------------------
# Başka data çekelim
df <- mtcars
head(df)
##                    mpg cyl disp  hp drat    wt  qsec vs am gear carb
## Mazda RX4         21.0   6  160 110 3.90 2.620 16.46  0  1    4    4
## Mazda RX4 Wag     21.0   6  160 110 3.90 2.875 17.02  0  1    4    4
## Datsun 710        22.8   4  108  93 3.85 2.320 18.61  1  1    4    1
## Hornet 4 Drive    21.4   6  258 110 3.08 3.215 19.44  1  0    3    1
## Hornet Sportabout 18.7   8  360 175 3.15 3.440 17.02  0  0    3    2
## Valiant           18.1   6  225 105 2.76 3.460 20.22  1  0    3    1
set.seed(145)
sampleindx <- sample(1:nrow(df),size = 0.75 * nrow(df))
trainset <- df[sampleindx,]
testset <- df[-sampleindx,]
nrow(trainset);nrow(testset)
## [1] 24
## [1] 8
names(df)
##  [1] "mpg"  "cyl"  "disp" "hp"   "drat" "wt"   "qsec" "vs"   "am"   "gear"
## [11] "carb"
chart.Correlation(trainset,histogram = TRUE, pch = 19)

md.pattern(df)# mice paketi eksik değer var mı diye kontrol.
##  /\     /\
## {  `---'  }
## {  O   O  }
## ==>  V <==  No need for mice. This data set is completely observed.
##  \  \|/  /
##   `-----'

##    mpg cyl disp hp drat wt qsec vs am gear carb  
## 32   1   1    1  1    1  1    1  1  1    1    1 0
##      0   0    0  0    0  0    0  0  0    0    0 0
model_mtcars <- lm(mpg~hp+wt+qsec+disp, data = trainset)
summary(model_mtcars) # mpg deki değişimin %81'i kullanılan bağımsız 
## 
## Call:
## lm(formula = mpg ~ hp + wt + qsec + disp, data = trainset)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.5872 -1.2896 -0.2797  0.9926  5.6381 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)   
## (Intercept) 26.293804  10.004311   2.628  0.01655 * 
## hp          -0.016949   0.024026  -0.705  0.48910   
## wt          -4.382860   1.387300  -3.159  0.00516 **
## qsec         0.553454   0.523101   1.058  0.30331   
## disp         0.001387   0.013345   0.104  0.91833   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.618 on 19 degrees of freedom
## Multiple R-squared:  0.8459, Adjusted R-squared:  0.8134 
## F-statistic: 26.07 on 4 and 19 DF,  p-value: 1.741e-07
# değişkenler ile açıklanabiliyor. İlk olarak modelin anlamlılığına bakarız.
# p değeri alphadan küçükse model anlamlı.
# Anlamsız değişkenler varken tahmin modelini yazmayız. İlk olarak anlamsızları çıkarmak
# lazım. P değeri 0.05 den en uzak olan değişkeni ilk çıkarırız.
model2 <- lm(mpg~hp+wt+qsec, data = trainset)
summary(model2)
## 
## Call:
## lm(formula = mpg ~ hp + wt + qsec, data = trainset)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.5904 -1.3105 -0.3799  0.9497  5.6648 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 26.27797    9.75264   2.694 0.013946 *  
## hp          -0.01593    0.02138  -0.745 0.464893    
## wt          -4.27776    0.92575  -4.621 0.000165 ***
## qsec         0.54484    0.50355   1.082 0.292131    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.553 on 20 degrees of freedom
## Multiple R-squared:  0.8458, Adjusted R-squared:  0.8226 
## F-statistic: 36.56 on 3 and 20 DF,  p-value: 2.611e-08
model3 <- lm(mpg~wt+qsec, data = trainset)
summary(model3) # model tamamen anlamlı oldu
## 
## Call:
## lm(formula = mpg ~ wt + qsec, data = trainset)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.9572 -1.6304 -0.0224  1.0025  5.4781 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  20.4383     5.7420   3.559  0.00185 ** 
## wt           -4.8539     0.5035  -9.640 3.67e-09 ***
## qsec          0.8488     0.2921   2.906  0.00845 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.526 on 21 degrees of freedom
## Multiple R-squared:  0.8415, Adjusted R-squared:  0.8264 
## F-statistic: 55.75 on 2 and 21 DF,  p-value: 3.984e-09
# mpg^ = 20,4 - 4.85wt + 0.84qsec
# Şimdi predictionları bulalım
model3pred <- predict(model3,testset)
RMSE(model3pred,testset$mpg) #caret paketi indirmen lazım.
## [1] 2.828962
MAE(model3pred,testset$mpg) # Eğer 2 tane model yapsaydık ikisinin bu değerlerine
## [1] 2.250896
# bakardık hangisi daha düşükse aslında o model daha iyi.
f_tablo <- qf(0.95,1,14)

Regresyon Varsayımları

Sabit Varyans (Homoscedasticity)

#Vize sonrası

#Sabit Varyans
x <- seq(0,1,length.out = 50)
y <- 2 * x + rnorm(50,sd = 0.5)
fit <- lm(y ~ x)
summary(fit)
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.76604 -0.38468  0.00023  0.32107  1.13397 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.01468    0.13152  -0.112    0.912    
## x            2.15138    0.22664   9.493 1.35e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.472 on 48 degrees of freedom
## Multiple R-squared:  0.6524, Adjusted R-squared:  0.6452 
## F-statistic: 90.11 on 1 and 48 DF,  p-value: 1.354e-12
# Sabit varyans durumu(homoscedasticity)

# artıklar y eksenine, yi tahmin değerleri x eksenine yazılarak 
# grafik yapılır.

par(mar = c(2, 2, 2, 2)) # Alt, Sol, Üst, Sağ kenar boşluklarını küçültür
plot(fit$fitted.values, residuals(fit), main = "Sabit Varyans", xlab = "ysapka", ylab = "Hatalar")

Değişen Varyans (Heteroscedasticity)

# Değişen varyans durumu(heterocedasticity)

x <- seq(0,1,length.out = 50)
y_het <- 2*x + rnorm(50,sd = 0.5*x)
fit2 <- lm(y_het ~ x)
plot(fit2$fitted.values,residuals(fit2), main = "Değişken Varyans",xlab = "ysapka",
     ylab = "Hatalar")
abline(h = 0, col ="red")

Doğrusal Olmayan Yapı

# Doğrusal olmayan yapı değişen varyans

x <- seq(0,1,length.out = 50)
y_nl <- x^2 + rnorm(50,sd = 0.05)
fit3 <- lm(y_nl ~ x)
plot(fit3$fitted.values,residuals(fit3), main = "Non lineer yapı",xlab = "ysapka",
     ylab = "Hatalar")
abline(h = 0, col ="red")

Savings Veri Seti - Varsayım Kontrolleri

# install.packages("remotes")
# remotes::install_version("faraway", version = "1.0.7")
library(faraway)
data(savings, package = "faraway")
names(savings)
model1 <- lm(sr ~pop15 + pop75 + dpi + ddpi, data = savings)
summary(model1)
# Model anlamlı yani lineer bir ilişki var fakat R^2 düşük o yüzden model açkıklama
# konusunda iyi bir model değil.

plot(model1$fitted.values,residuals(model1), main = "Sabit varyans kontrolu",xlab = "ysapka",
     ylab = "Hatalar")
abline(h = 0, col ="red")

# Breusch-Pagan testi
# Yapılarak tam olarak belli olmayan grafiği anlayacağız.
# H0: varyans sabittir.
# H1: Varyans değişkendir. (sorun var)
# install.packages("lmtest")
library(lmtest)
bptest(model1)
# BP değeri p değerinden büyük olduğundan H0 reddedilemez. Yani sabit varyanslı
# Çıktı ve iyi oldu.

# Değişken varyans olsaydı ilk yapacağımız şey dönüşüm yapmak olacaktır.
# Bağımlı değişkene ilk akla gelen logaritmik dönüşüm yapılır.
# y nin logaritması alınır ve yeni model oluşturulur.(zor değil klasik)
# Bundan sonra karekök döüşümü vs. yapılır hangisi olursa 
# Uyguladğımız dönüşümün ters dönüşümünü de yapmalıyız.
# Bu işe yaramazsa ağırlıklı en küçük kareler yöntemi ile model değiştirilmeye çalışılır.
# Bu konu şimdinin konusu değil.

Normallik Varsayımı

# NOT: Üstteki chunk eval=FALSE olduğu için (yani model1 yaratılmadığı için) 
# buranın çalışması hata verir. Bu yüzden burası da eval=FALSE yapılmıştır.

#Normallik varsayımı:
# Modelin artıklarının nırmal dağılıma uyup uymadığı ile ilgilidir.
# Bunu test etmek için hem grafiksel hem de istatistiksel yöntemler kullarınlıt.

# Grafik olarak qqplot sıklıkla tercih edilir.
# Test olarak da Shapiro testi kullanılabilir.

qqnorm(residuals(model1))
qqline(residuals(model1), col = "red", lwd = 2)
# Gözlemler çizgi üzerinde çok yakınlar yani normallikten bir sapma yok.

#Shapirotest
shapiro.test(residuals(model1))
# H0: normal dağılır
# H1: normallikten sapar.
# W değeri p değerinden büyük çıktığından H0 reddedilemez. Yani normaldir.

# Eğer normallikten sapsaydı ilk olarak gözlem sayısı artırılabiliyorsa artırılır.
# Yoksa y üzerinde dönüşüm yapılır bazen de x üzerinde de dönüşüm yapılabilir.

Otokorelasyon (İlişkili Hatalar)

# NOT: Üstteki chunk (model1) çalışmadığı için burası da çalıştırılmamalıdır.

# İlişkili hatalar varsayımı
# Bir gözlemdeki hatanın kendinden önceki veya sonraki hatayla ilişkili olması durumudur.
# H0: otokorelasyon yoktur. (ro = 0)
# H1: otokorelasyon vardır. (ro != 0)
# Durbin-Watson testi yapılır. 
# 0 - 1,5 arası çıkarsa pozitif otokorelasyon
# 1,5 - 2,5 arası ise otokorelasyon yok.
# 2,5 - 4 arası ise negatif otokorelasyon
# ilk grafikle incelersek
n <- length(residuals(model1))
plot(tail(residuals(model1), n-1) ~ head(residuals(model1),n-1)
     ,xlab = expression(hat(epsilon)[i]),
     ylab = expression(hat(epsilon)[i+1]))
abline(h = 0, v = 0, col = grey(0.75))
# Grafikteki saçılım rastgele olduğundan otokorelasyon yok deriz.

# test ile de bakalım
dwtest(model1)
# 1.5 ve 2.5 arası çıktı. Yani H0 reddedilemez. Otokorelasyon yok.

Stepwise Regresyon (Adımsal Regresyon)

# Stepwise Regression(Adımsal Regresyon)

# rEGRESYON MODELİ OLUŞTURULURKEN AMACIMIZ BAĞIMSIZ DEĞİŞKENLERİ EN İYİ ŞEKİLDE SEÇMEKTİR
# Stepwise yönteminde değişkenler belirli bir sistematikle çıkartılarak ya da 
# eklenerek farklı regresyon modelleri oluşturulup o modeller arasından en iyisini
# seçmemizi sağlayan yöntemdir.

# 1-) Forward
# 2-) Backward
# 3-) Both

# Forward - Tek tek modele değişkenlerin eklenmesiyle elde edilir.
# Başta sadece B0 ile modelimiz var, sonra bütün değişkenleri tek tek ekliyor.
# en düşük AIC değerine sahip değişkenleri sırayla tek tek ekleriz.

# Backward - Full modelden tek tek değişken çıkarılarak işlem yapılır.
# AIC değeri iyice azalana kadar tam modelden bağımsız değişkenler çıkarılır.

# Both - değişkenler duruma göre hem eklemip hem çıkarılabilir.

#Forward
step(lm(sr~1,data = savings),direction = "forward",
     scope = ~pop15 + pop75 + dpi + ddpi)

# çıktı olarak ilk olarak hangilerini eklediğini ve neye sebep olduğunu görüyoruz.
#AIC değerini en çok düşüren pop15 ilk eklenmiş.
# Sonra bu öodelin üstüne diğer değişkenleri ayrı ayrı eklediğinde hangisi AIC en çok
# düşürürse onu da ekler.
# en son AIC değeri düşmediğinde modeli bitirir.

#Backward
step(lm(sr ~ pop15 + pop75 + dpi + ddpi,data = savings))
# ekstra belirtmeye gerek yok direction diye bunu
# Tek tek değişkenleri AIC değerini, modelden çıktığı takdirde düşürüyorsa onu çıkarırız.
# son model en son oluşur.

#Both
step(lm(sr~1,data = savings),direction = "both",
     scope = ~pop15 + pop75 + dpi + ddpi)

# İlk olarak başlangıçç modeline en düşük AIC değerli değişkeni ekledi.
# Artık başta eklediği değişkeni çıkarırsak ne olur onu da gösteriyor.
# en son kararı sonda yine verir.

Dummy Değişkenler

Teori ve Örnekler

# Dummy değişken
# Bağımsaz değişken nitel olduğu zaman kullanılır.
# Regresyon analiiznde sıklıkça kullanılan değişkenler niceldir.
# sıcaklık,uzaklık,basınç vs.
# bazı durumlarda regresyonda bağımsız değişkenler olarak nitel (kategorik)
# değişkenler kullanmak gerekir. istihdam durumu, vardiya(gündüz,akşam,gece)
# cinsiyet.
# Nitel bir değişkeni yanıt(y) üzerinde sahip olabileceği etkiyi dikkate almak için
# dummy değişkenler kullanmak gerekir. 
# Eğer cinsiyetteki gibi erkek ve kadın olarak 2 tane kategori varsa 2-1 tane dummy değişken
# oluştururuz. Yeni bir X3 dummy değişkeni oluşturduk ve kadına 0 erkeğe ise 1 sayısını verdik.
# Örnek; X = eğitim düzeyi, alt kategorileri ise (ilköğretim,lise,lisans,lisansüstü)
# 4-1 = 3 tane dummy değişken oluşturmalıyız.

x <- data.frame(x1 = c(0,1,0,0), x2 = c(0,0,1,0), x3 = c(0,0,0,1), EgitimDuzeyi = c("ilkogretim","lise","lisans","lisansustu"))
x
##   x1 x2 x3 EgitimDuzeyi
## 1  0  0  0   ilkogretim
## 2  1  0  0         lise
## 3  0  1  0       lisans
## 4  0  0  1   lisansustu
# Bu şekilde her alt kategori için 3 tane kukla değişken oluşturmuş olduk.

#Örnek:Bir makine mühendisinin tonra üzerinde kullanılan kesici bir aletin ömrünü
#dk başına torna hızına(X1) ve kullanılan kesici alet türüne göre incelemek istediği
#varsayılsın.Alet türü A ve B olmak üzere ikiye ayrılmaktadır.Buna göre uygun modeli
#belirleyiniz.
Y  <- c(18.73,14.52,17.43,14.54,13.44,24.39,13.34,22.71,12.68,19.32,30.16,27.09,25.40,
        26.05,33.49,35.62,26.07,36.78,34.95,43.67)
X1  <- c(610,950,720,840,980,530,680,540,890,730,670,770,880,1000,760,590,910,650,810,
         500)
ALETturu <- c("A","A","A","A","A","A","A","A","A","A","B","B","B","B","B","B","B","B","B","B")

# Bu önekte aletturu deeğişkeni 2 farklı kategoriye sahip yani bir tane kukla değişken oluştururuz.
# X2 = 0 <- A alet turu
#      1 <- B alet turu

# Normal model = Y = B0 + B1X1 + B2X2 + Ebsulon
# A alet turu olan model = Y = B0 +B1X1 + Ebsulon
# B alet turu olan model = Y = B0 + B2 + B1X1 + Ebsulon

# modeli oluşturusak 
# ysapka <- 36,986 - 0,027X1 + 15,004X2 modeli ortaya çıkar.
# Bu modelde B alet türünün A alet türü üzerinde 15 birimlik daha fazla etkiye sahip ol-
# duğunu görüyoruz. Çünkü A alet türünün modelinde  X2 değişkeni yok, fakat B alet türünde var.
# ana modelde olan X2 bu dediğimiz farkı anlatır bu yüzden.



# örnek
# yaş ile şiddetli depresyon tedavisinde kullanılan 3 farklı yöntemin (A,B,C)
# Tedavi etkinliği üzerindeki etkisi inceleniyor. Bu amaçla hastalıpın tanısına ve şiddetine göre 
# 3 farklı yöntem uygulanan 36 hasta rastgele seçiliyor.
# Buna göre uygun modeli belirle ve yorumla. Y tedavi etkinliği, X1 yaş
Y <- (c(56,40,41,45,60))
X1 <- c(21, 23, 30, 19, 18)
tedaviyontemi <- c("A","B","B","C","A")
a <- data.frame(X2 = c(0,1,0), x3 = c(0,0,1), Yontem = c("C","A","B")) # 2 kukla oluşturduk ve değerleri tek tek yazdık.
# ysapka <- 22,29059 + 0,66446X1 + 10,25276X2 + 0,44518X3
# C yontemi için regresyon denkleml <- x2 = x3 = 0    ysapka <- 22,29059 + 0,66446X1
# A yöntemi için regresyon denklemş <- x2 = 1, x3 = 0    ysapka <- 22,29059 + 0,66446X1 + 10,2527
#                                                               = 32,54335 + 0,66446X1
# B yontemi için regresyon denklemi <- x2 = 0, x3 = 1    ysapka <- 22,29059 + 0,66446X1 + 0,44518
#                                                               = 22,73577 + 0,66446X1
# Buradan A yönteminin en etkili yöntem olduğunu görebiliyoruz.

R Uygulaması

# R orneği

data <- data.frame(
  y = c(18.73,14.52,17.43,14.54,13.44,24.39,13.34,22.71,12.68,19.32,30.16,27.09,25.40,
        26.05,33.49,35.62,26.07,36.78,34.95,43.67),
  x1= c(610,950,720,840,980,530,680,540,890,730,670,770,880,1000,760,590,910,850,810,
        500),
  Alet_turu = c("A","A","A","A","A","A","A","A","A","A","B","B","B","B","B","B","B",
                "B","B","B")
)

data
##        y   x1 Alet_turu
## 1  18.73  610         A
## 2  14.52  950         A
## 3  17.43  720         A
## 4  14.54  840         A
## 5  13.44  980         A
## 6  24.39  530         A
## 7  13.34  680         A
## 8  22.71  540         A
## 9  12.68  890         A
## 10 19.32  730         A
## 11 30.16  670         B
## 12 27.09  770         B
## 13 25.40  880         B
## 14 26.05 1000         B
## 15 33.49  760         B
## 16 35.62  590         B
## 17 26.07  910         B
## 18 36.78  850         B
## 19 34.95  810         B
## 20 43.67  500         B
class(data$Alet_turu)
## [1] "character"
data$Alet_turu <- as.factor(data$Alet_turu)# faktor ataması yaptık dummy atamak için bu gerekli
model <- lm(y ~ x1 + Alet_turu,data = data)
summary(model)
## 
## Call:
## lm(formula = y ~ x1 + Alet_turu, data = data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.4295 -1.9523 -0.2735  1.9542  6.7344 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 35.611851   4.016197   8.867 8.75e-08 ***
## x1          -0.024768   0.005173  -4.788 0.000171 ***
## Alet_turuB  15.486742   1.552521   9.975 1.61e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.457 on 17 degrees of freedom
## Multiple R-squared:  0.871,  Adjusted R-squared:  0.8558 
## F-statistic: 57.38 on 2 and 17 DF,  p-value: 2.758e-08

WeatherAUS Veri Seti - Kapsamlı Analiz

Veri Hazırlama ve Missing Değer İşleme

#Yeni dataset yükledik.

View(weatherAUS)
table(weatherAUS$Location)

# weatherAUS <- read.csv("") içine dosyanın linkini yazıp slashları tersine çevirince de olur.

head(weatherAUS)
# Bu dataset regresyon için uygun değil. Missing değerler için bunu kullanıyoruz.

df_albury <- weatherAUS[weatherAUS$Location == "Albury",]
nrow(df_albury)
names(df_albury)
df_albury <- df_albury[c("Humidity9am","MinTemp","MaxTemp","WindSpeed9am",
                         "Pressure9am","Temp9am","Rainfall")]

View(head(df_albury))

#ilk olarak değişkenler arası ilişkilere bakalım.
cor(df_albury)
# missingler var o yüzden korelasyon yapamadık.
sum(is.na(df_albury))
#83 missing değer var.
# bu missing değerleri atarak bakarsak korelasyona
cor(na.omit(df_albury)) # böyle olabilir fakat bu datada bunu yapmayacağız.

library(mice)
md.pattern(df_albury) # missing değerlere baktık
# Rainfall değişkeninde 16 gözlem eksik buna bakınca mesela.
# MaxTemp de 8 eksik gözlem var.
# mesela en alttan bir üst bölmede 6 farklı değişkende 2şer tane boş var.

# Şimdi bu eksik değerleri dolduracağız.
?mice()# Burada doldurma için kullanabileceğimiz methodlar vardır.
# pmm, logreg, polyreg, polr gibi.

set.seed(123)
imputed <- mice(df_albury, n = 3) # Eğer method uygulamak istersek method = "pmm" gibi yapabiliriz.
# Ama zaten kendi hali ile üretim yaparsak da olur.
names(imputed) # paket bu değerleri atamış
imputed$imp# missing değerler ürettik.
imputed$method # hangi  değişkenlere hangi methodla değer ürettik onu gösterir.

df_albury_Imp <- complete(imputed,3)
View(df_albury_Imp)
md.pattern(df_albury_Imp)
# missing değerleri normal datamıza ekledik ve hiç missing kalmadı.

Model Kurma ve Değerlendirme

# Model oluşturalım.
# eğitim ve test olarak ayıralım.

set.seed(145)
sampleIndex <- sample(1:nrow(df_albury_Imp),
                      size = 0.8*nrow(df_albury_Imp))

trainset <- df_albury_Imp[sampleIndex,]
testset <- df_albury_Imp[-sampleIndex,]
nrow(trainset);nrow(testset)

# modeli yazalım.
cor(df_albury_Imp)
names(df_albury_Imp)
model <- lm(Humidity9am ~ MinTemp + MaxTemp + WindSpeed9am + Pressure9am  + Temp9am + Rainfall,data = trainset)
summary(model)
#  Pressure9am anlamsız çıkaralım.
# çıkarmadan önce de viflerine de bakabiliriz.
vif(model) # 5 den büyük olan değer çıkarılabilir.

model2 <- lm(Humidity9am ~ MinTemp + MaxTemp + WindSpeed9am + Temp9am + Rainfall,data = trainset)
summary(model2)
# model anlamlı.

plot(model2)
# değişen varyans sorunumuz yok. 1. grafik
# normaliteden de sapma olmadığını görüyoruz fazla. 2. grafik
# modelin tespit ettiği aykırı gözlemleri gösterir 4. grafik. Aykırı değerler değilse
# çok küçük değerler olması lazım fakat burada birkaç tane aykırı değer olduğunu görüyoruz.
# Bu gözlemin model üzerinde çok büyük etkiye sahip olduğunu anlıyoruz. etkin bir gözlem.
# çizgiyi kendine doğru çeken bir gözlem varsa o etkin bir gözlemdir. Bu gözlem
# modeldeyken ve modelde değilken incelenmeli.
# kendine çizgi çekmeyip yüksek olan değerler ise aykırı değerler oluyor.

# aykırıları çıkarmadan predictlere bakalım
predictions <- predict(model2, testset)
head(predictions)

library(caret)
R2(predictions,testset$Humidity9am) # tahminlerle gerçek değerler arası ilişkiye bakar.
# 0.71 çok kötü bir ilişkş değil diyebiliriz.
RMSE(predictions, testset$Humidity9am) # y değeri ile çalışıyoruz tabi ki.
# bu değer başka bir değerle kıyaslarsak anlam ifade eder tabi.
MAE(predictions, testset$Humidity9am)

Aykırı Değer Analizi

# Aykırı değer kontrolü

strandardized_residuals <- rstandard(model2)
summary(strandardized_residuals)
olcut1index <- which(abs(strandardized_residuals) > 2)
length(olcut1index)

# Alternatif aykırı belirleme(cook distance)
dist <- cooks.distance(model2)
olcut1 <- mean(dist) * 3
olcut2 <- 4/length(dist)

olcut1index <- which(dist > olcut1)
olcut2index <- which(dist > olcut2)
length(olcut1index);length(olcut2index)

plot(1:length(dist),dist,type = "p", ylim = range(dist)*c(1,1))
plot(1:length(dist),dist,type = "p", ylim = range(dist)*c(1,0.07))
abline(h = olcut1, col = "red")
# üstte kalanlar aykırılar.
abline(h = olcut2, col = "green")

outliers <- which(dist > olcut2 & abs(strandardized_residuals) > 2)
#  Bunları ana setten çıkarmalıyız.
trainsetrem <- trainset[-outliers,]
nrow(trainset);nrow(trainsetrem)

model3 <- lm(Humidity9am ~ MinTemp + MaxTemp + WindSpeed9am + Temp9am + Rainfall,data = trainsetrem)
summary(model3)

prediction3 <- predict(model3,testset)
R2(prediction3,testset$Humidity9am)
RMSE(prediction3,testset$Humidity9am)
MAE(prediction3,testset$Humidity9am)
# Aykırılar olmadan ve aykırılar varkenki tahmminleri kıyasladık aykırıları attığımız model daha iyi.
R2(predictions,testset$Humidity9am)
RMSE(predictions,testset$Humidity9am)
MAE(predictions,testset$Humidity9am)

summary(model3)
summary(model2)
library(car)
vif(model3)
# En yüksek vif değerli olanı direkt çıkarabiliriz.
modelvif_1 <- lm(Humidity9am ~ MinTemp + MaxTemp + WindSpeed9am + Rainfall,data = trainsetrem)
vif(modelvif_1)
# diğerlerinin de sorunlları çözüldü
summary(modelvif_1)
# Nihai modelimiz bu oldu fakat bu model zaten regresyon için iyi bir model olmadığından 
# R^2 hala düşük.
bptest(modelvif_1)

Swiss Veri Seti - Detaylı Analiz

# Farklı bir model ile regresyon yapalım.
data(swiss)
head(swiss)
##              Fertility Agriculture Examination Education Catholic
## Courtelary        80.2        17.0          15        12     9.96
## Delemont          83.1        45.1           6         9    84.84
## Franches-Mnt      92.5        39.7           5         5    93.40
## Moutier           85.8        36.5          12         7    33.77
## Neuveville        76.9        43.5          17        15     5.16
## Porrentruy        76.1        35.3           9         7    90.57
##              Infant.Mortality
## Courtelary               22.2
## Delemont                 22.2
## Franches-Mnt             20.2
## Moutier                  20.3
## Neuveville               20.6
## Porrentruy               26.6
sum(is.na(swiss))
## [1] 0
set.seed(741)
indeks <- sample(1:nrow(swiss), 0.8 * nrow(swiss))
trainset <- swiss[indeks,]
testset <- swiss[-indeks,]
dogurganlik <- lm(Fertility~.,data = trainset)
summary(dogurganlik)
## 
## Call:
## lm(formula = Fertility ~ ., data = trainset)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -14.4835  -3.9833   0.5368   3.9198  11.6717 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      72.43372   11.00335   6.583 2.36e-07 ***
## Agriculture      -0.17999    0.06494  -2.772  0.00934 ** 
## Examination      -0.36110    0.25107  -1.438  0.16039    
## Education        -0.82707    0.18152  -4.556 7.63e-05 ***
## Catholic          0.07113    0.03368   2.112  0.04284 *  
## Infant.Mortality  0.91272    0.38587   2.365  0.02445 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.234 on 31 degrees of freedom
## Multiple R-squared:  0.7485, Adjusted R-squared:  0.708 
## F-statistic: 18.46 on 5 and 31 DF,  p-value: 1.755e-08
dogurganlik2 <- lm(Fertility~. - Examination,data = trainset)
summary(dogurganlik2)
## 
## Call:
## lm(formula = Fertility ~ . - Examination, data = trainset)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -13.6617  -5.6186   0.9912   3.9705  13.2392 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      64.37597    9.62706   6.687 1.51e-07 ***
## Agriculture      -0.15587    0.06377  -2.444  0.02021 *  
## Education        -0.97771    0.15070  -6.488 2.67e-07 ***
## Catholic          0.09972    0.02764   3.607  0.00104 ** 
## Infant.Mortality  0.98137    0.38924   2.521  0.01687 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.337 on 32 degrees of freedom
## Multiple R-squared:  0.7318, Adjusted R-squared:  0.6982 
## F-statistic: 21.82 on 4 and 32 DF,  p-value: 9.134e-09
library(car)
vif(dogurganlik2)
##      Agriculture        Education         Catholic Infant.Mortality 
##         1.959682         1.846227         1.197994         1.129013
plot(vif(dogurganlik2))

library(lmtest)
bptest(dogurganlik2)
## 
##  studentized Breusch-Pagan test
## 
## data:  dogurganlik2
## BP = 5.3653, df = 4, p-value = 0.2518
dwtest(dogurganlik2)
## 
##  Durbin-Watson test
## 
## data:  dogurganlik2
## DW = 2.0881, p-value = 0.6001
## alternative hypothesis: true autocorrelation is greater than 0
sum(dogurganlik$residuals)
## [1] 7.105427e-15
shapiro.test(dogurganlik2$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  dogurganlik2$residuals
## W = 0.95695, p-value = 0.1615
ks.test(dogurganlik2$residuals,"pnorm",mean(dogurganlik2$residuals),sd(dogurganlik2$residuals))
## 
##  Exact one-sample Kolmogorov-Smirnov test
## 
## data:  dogurganlik2$residuals
## D = 0.18125, p-value = 0.1552
## alternative hypothesis: two-sided
# install.packages("nortest")
library(nortest)
ad.test(dogurganlik2$residuals)
## 
##  Anderson-Darling normality test
## 
## data:  dogurganlik2$residuals
## A = 0.8663, p-value = 0.02362
lillie.test(dogurganlik2$residuals)
## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  dogurganlik2$residuals
## D = 0.18125, p-value = 0.003497
cooks.distance(dogurganlik2)
##        Broye      Gruyere      Conthey     Grandson        Nyone Paysd'enhaut 
## 1.640150e-02 4.540309e-03 6.356754e-03 9.773015e-05 1.454079e-02 3.344881e-03 
##     Le Locle       Herens         Oron   Courtelary    Echallens   Neuveville 
## 1.126453e-02 9.656383e-04 2.625916e-04 4.183585e-02 1.559117e-02 5.642959e-02 
##      Payerne  Rive Gauche       Boudry      Monthey ValdeTravers        Aigle 
## 6.456880e-03 1.920035e-01 2.740435e-03 1.491984e-03 5.333781e-02 1.007644e-02 
##      Veveyse     Lausanne      Yverdon       Moudon       Lavaux      Aubonne 
## 1.949550e-02 4.249016e-04 1.600075e-02 5.412886e-02 7.450996e-04 4.166076e-04 
##    La Vallee    Entremont     Cossonay Franches-Mnt     Avenches       Morges 
## 5.601814e-03 2.463409e-02 1.644507e-02 1.097586e-01 7.500332e-03 1.653266e-03 
## V. De Geneve       Sarine     Delemont La Chauxdfnd   Val de Ruz   St Maurice 
## 1.158128e-02 2.320647e-02 8.973038e-03 8.779484e-02 1.278400e-02 3.379410e-02 
##   Porrentruy 
## 2.283504e-01
plot(cooks.distance(dogurganlik2))

distin <- cooks.distance(dogurganlik2)
olcum_birinci <- mean(distin) * 3
olcum_ikinci <- 4/length(distin)

olcut_1_indeks <- which(distin > olcum_birinci)
olcut_2_indeks <- which(distin > olcum_ikinci)
print(olcut_1_indeks);print(olcut_2_indeks)
##  Rive Gauche Franches-Mnt   Porrentruy 
##           14           28           37
##  Rive Gauche Franches-Mnt   Porrentruy 
##           14           28           37
plot(1:length(distin),distin,type = "p",ylim = range(distin) * c(1,1))
abline(h = olcum_birinci, col = "red")
abline(h = olcum_ikinci, col = "green")

standartlastir <- rstandard(dogurganlik2)
outliers <- which(distin > olcum_birinci & abs(standartlastir) > 2)
print(outliers)
## Rive Gauche 
##          14
trainset_rem <- trainset[-outliers,]

aykirisiz <- lm(Fertility ~. - Examination, data = trainset_rem)
summary(aykirisiz)
## 
## Call:
## lm(formula = Fertility ~ . - Examination, data = trainset_rem)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -11.4939  -5.3926   0.9368   3.6175  12.5955 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      63.20339    8.91902   7.086 5.83e-08 ***
## Agriculture      -0.15157    0.05903  -2.568 0.015275 *  
## Education        -0.87502    0.14523  -6.025 1.14e-06 ***
## Catholic          0.10705    0.02574   4.159 0.000234 ***
## Infant.Mortality  0.98386    0.36013   2.732 0.010298 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.863 on 31 degrees of freedom
## Multiple R-squared:  0.7355, Adjusted R-squared:  0.7013 
## F-statistic: 21.55 on 4 and 31 DF,  p-value: 1.386e-08
tahmin <- predict(dogurganlik2,testset)
tahmin2 <- predict(aykirisiz,testset)
library(PerformanceAnalytics)
library(caret)
RMSE(tahmin,testset$Fertility)
## [1] 9.784668
MAE(tahmin,testset$Fertility)
## [1] 8.768165
R2(tahmin,testset$Fertility)
## [1] 0.7227114
RMSE(tahmin2,testset$Fertility)
## [1] 9.604335
MAE(tahmin2,testset$Fertility)
## [1] 8.715065
R2(tahmin2,testset$Fertility)
## [1] 0.7469627

Lojistik Regresyon

Teori ve Kavramlar

Lojistik regresyon, sınıflama problemleri için kullanılan bir regresyon yöntemidir.

Temel Kavramlar:

  • Binary Lojistik Regresyon: İkili sonuç değişkeni için
  • Multinomial Lojistik Regresyon: İkiden fazla kategori için
  • Ordinal Lojistik Regresyon: Sıralı kategoriler için

Odds ve Odds Ratio:

  • Odds = p / (1-p)
  • Odds Ratio (OR):
    • OR = 1: Etki yok
    • OR > 1: Risk artırıcı
    • OR < 1: Risk azaltıcı

Default Veri Seti Analizi

library(glmnet)
library(cli)
library(tidyverse)
library(modelr)
library(broom)
library(ISLR)
# install.packages("pscl")
library(pscl)

data <- Default
View(data)

set.seed(123)
sample1 <- sample(1:nrow(data), size = 0.6*nrow(data))
train <- data[sample1,]
test <- data[-sample1,]
nrow(train);nrow(test)
names(data)
model1 <- glm(default ~ balance, family = "binomial", data = train)
summary(model1)
exp(coef(model1))
# şimdi odds oranlarını elde ettik.
#balance değerinin neredeyse hiç etkili olmadığını gördük.
# balance daki 1 birimlik artış kişinin gecikmeye düşmesi olasılığını 1.0056 kat artırır.
predict(model1, data.frame(balance = c(1000,2000)), type = "response")
# hesabında 1000 dolar olanın gecikmeye düşme olasılığı 0.055, hesabında 2000 doları olanın
#gecikmeye düşme olasılığı 0.596 çıktı. Bunlar y şapkaların olasılıkları.
#Biz ise gecikmeyedüşüyor mu düşmüyor mu diye yorumlamalıyız evet hayır olarak.

# Dummy değişken atayalım
class(data$student) # faktör olmasaydı faktör ataması yapardık.
model2 <- glm(default ~ student, family = "binomial", data = train)
# Şimdi öprenci olan ve olmayanların gecikmeye düşüp düşmemesine bakıcaz.
summary(model2)
exp(coef(model2))
# Öğrenci oluşu öğrenci olmayışına göre gecikmeye düşme ihtimalini 1.48 kat artırıyor.
# Yani öğrenci olmak bir risk.

model3 = glm(default~student+balance+income, family = "binomial", data = train)
summary(model3)
exp(coef(model3))

caret::varImp(model3)

Model Değerlendirme

# Degerlendirmeler:

# 1-) olabilirlik oran testi: herhangi bir modele bagimsiz degiskenler eklemek hemen hemen her zaman model uyumunu 
#iyilestirebilir. yani bos modele (sadece beta0 dan olusan) kiyasla logaritmik olabilirlik artacak ve model sapmasi
#azalacak. ancak model uyumunda gozlenen farkin istatistiksel olarak anlamli olup olmadigi test edilmelidir. bunun icin
#anova kullanilabilir buradaki hipotezler ;
# ho = eklenen bagimsiz degisken modele katki saglamamistir.
# h1 = eklenen bagimsiz degiskenler modele katki saglamistir.

Anova = anova(model1, model3, test = "Chisq")

#ln(ML1) uydurulan model i??in log olabilirlik de??eri.
#ln(ML0) b0 l?? bo?? model i??in log olabilirlik de??eri.
# R2 McFadden = 1 - ln(ML1) / ln(ML0) max de??er 0.6 0.4 ??n ??zeri yeterli.

list(model1=pscl::pR2(model1)["McFadden"],
     model2=pscl::pR2(model2)["McFadden"],
     model3=pscl::pR2(model3)["McFadden"])
# sozde R2'ler arasinda en iyisi model3, sonra model1 en kotusu ise model2. McFadden degeri 0.4 den buyukse guclu model.

# 3-) siniflama oranlari: modeller gelistirilirken en kritik olcu modelin ornek disi ornekler (out of sample) icin hedef
#degiskeni tahmin etmede ne kadar iyi oldugunu belirlemektir.

test.predicted.m1 = predict(model1, newdata = test, type = "response")
test.predicted.m2 = predict(model2, newdata = test, type = "response")
test.predicted.m3 = predict(model3, newdata = test, type = "response")
cbind(test.predicted.m1,test.predicted.m2,test.predicted.m3) # olasılık değerlerini aldık.
# Bu olasılıkları etiketlemeliyiz. (gecikmeye düştü mü düşmedi diye)


table(test$default)
table(train$default) # Dengesiz bir parçalanma var. No grubu çok fazla

predictClass = ifelse(test.predicted.m1>0.5,"Yes","No")
predictClass = as.factor(predictClass)
caret::confusionMatrix(predictClass, reference = test$default, positive = "Yes")
# Baştaki confusion matrix
# Gerçekte no iken no tahmin ettiklerimiz, yes ve no olarak olan değerler yanlış tahminlerdir.
# Accuracy: doğru karar verilen durumların oranı.
# No'ları güzel oranda tespit etmişiz fakat yes'leri pek de iyi tahmin edememişizdir.
# O nedenle accuracy pek işe yaramamış
# No information rate: hiç model uygulamadan sadece dataya bakarak yapılabilecek tahmninin 
# doğruluk oranı. No'lar fazla oldupundan bu değer yüksek çıkmış.
# P değeri 0.05 ten biraz küçük yani H0 reddedilir fakat 0.05'e çok yakın yani sıkıntılı bir model.

# Kappa istatistiği: Gerçek sınıflarla predictionlar incelenerek elde edilen bu
# 2 vektörün ne kadar benzer olduklarını ölçmemizi sağlar.
# 0.1 - 0.2 bandı : Kötü
# 0.2 - 0.4 bandı: Kötü ama idare eder.
# 0.4 - 0.6 bandı: orta
# 0.6 - 0.8 bandı: iyi
# 0.8 - 1 bandı: mükemmel

#Mcnemar testi:

# H0: False pozitif ve False negatif birbirine benzer.
# H1: Birbirine benzer değildir.
# Sonuç olarak burada p 0.05 den küçük. H0 red. Yani birbirlerine benzer değiller.



caret::confusionMatrix(predictClass, reference = test$default, positive = "Yes",
                       mode = "prec_recall")
# Precision değerini aldık. 0.66 yani pek iyi model değil.
# Recall ise doğru pozitiflerin gerçek pozitiflere oranı. 0.28 yani pozitifleri iyi yakalayamıyoruz.


cbind(predicted = as.character(predictClass))