R Dersleri - Ali Onur Gitmez
R Dersi 12
Ali Onur Gitmez
2024-07-15
Regression
İstatistik analizinin temel yapıtaşı olaraksa regresyon var. Yapacağımız neredeyse tüm analizler ve çalışmalar bu yönteme bağlı olacak. Linear regreson temel olarakbağımsız değişken(independent variable) ve bağımlı değişken(dependent variable) arasındaki ilişkiyi ölçer. Amacı X Y düzlemindeki en iyi çizgiyi bulmaktır. Bu çizgi de X kullanarak Y’yi en iyi şekilde tahmin etmemize yardımcı olur. Linear Regresyon metoduna OLS yani Ordinary Least Squares de deriz. Bu observed ve predicted değerlerimizin arasındaki farkın karelerinin alınıp toplanması ile oluşur. Bu karelerin toplamının en düşük olduğu çizgi bizim en iyi fit eden çizgimizdir ve regresyon sonuç olarak bize bunu verir.
Genellikle dependent variableımızın continious veya en azından continiousa yakın discrete olması gerekir. Independent variablelarımız ise discrete veya continious olabilir. Mesela bir kişinin maaşını tahmin etmek istiyorsak bu metodu kullanırız. Kişinin eğitim seviyesine, kaç yıl çalıştığına ve belli sertifikalara sahip olup olmamasına göre biz kişinin maaşını tahmin edebiliriz.
Temel olarak 2 çeşit linear regresyon vardır.
- Simple Linear Regression: 1 bağımlı ve 1 bağımsız değişken arasındaki ilişkiyi ölçer
- Multiple Linear Regression: 1 bağımlı ve 1 üstü bağımsız değişken arasındaki ilişkiyi ölçer.
Formülü şu şekildedir:
\(\hat{Y} = \hat{\beta}_0 + \hat{\beta}_1 x_1 + \hat{\beta}_2 x_2 + \cdots + \hat{\beta}_k x_k + \hat{\epsilon}\)
Bazı alanlarda -özellikle ekonometri- ise şu şekilde yazılır
\(\hat{Y} = \hat{\alpha} + \hat{\beta}_1 x_1 + \hat{\beta}_2 x_2 + \cdots + \hat{\beta}_k x_k + \hat{\epsilon}\)
Bunlar arasında teknik olarak fark olmasa bile ikisini de görmeniz mümkündür.
Bu formülde \(\hat{\beta}_0\) veya \(\hat{\alpha}\) olarak belirttiğimiz şey constant veya intercept olarak adlandırılır. Tüm bağımsız değişkenler 0 iken yani bağımsız değişkenlerimizin bir etkisi olmadığını varsaydığımız durumlarda bağımlı değişkenin değerini gösterir. Bazı durumlarda bağımsız değişkenlerimiz 0 olamayabilir. Bu gibi durumlarda interceptimizin bizim için bir anlamı olmaz ama istatistik hesaplamasında işe yarar.
\(\hat{\beta}_1 x_1\) diye belirttiğimiz şeyler ise bağımsız değişkenlerimiz. X kısmı doğrudan değişkenimizin ne olduğunu bekirtirken \(\hat{\beta}\) kısmı da eğrinin eğimini gösteriyor ve birlikte değişkenimizin etkisini görebiliyoruz.
\(\hat{y}\) ise bağımlı değişkenimiz. \(\hat{\epsilon}\) da modelimiz gerçek dünyada bütün veriyi açıklayamayacağı için error termümüz. Yani predicted ve observed değerlerimiz arasındaki fark. Bunu değişkenlerimiz tarafından açıklanmayan kısmı beirtmek için kullanıyoruz.
##Assumptionlar:
Linearity: Bağımlı ve bağımsız değişken arasındaki ilişki linear olmalı. Linear olmayan ilişkileri linear regreson ile tespit edemeyiz. Farklı yöntemlere bakmamız gerekir. Bunu scatterplot ile kontrol edebiliriz. Eğer sorun varsa başka bir regresyon yöntemi yapmamız gerekir.
No autocoreelation: Errorler birbirinden bağımsız şekilde dağılamlı. Bu özellikle time-series verilerinde bir sorun oluşturur. Mesela borsa fiyatlarını tahmin etmeye çalışıyorsak bugünün borsa fiyatları dünden, yarının ise bugünden etkilenecek. Buna da dikkat etmemiz lazım. Durbin-Watson test uygularız. Bunu yöntem değiştirerek çözebiliriz. Cross sectional verilerde çok problem olmayabilir.
Homoscedasticity: Residualların yani beklenen ve gözlemlenen veri arasındaki farkın varyanslarının bağımsız değişkenin tüm seviyeleri için eşit dağılması gerekir ama bazen bu gerçekleşmez. Bunu da residual ve fitted valueların plotuna bakarak yapabiliriz. Mesela tecrübe ve maaş arasındaki ilişkiyi inceliyoruz. İşe yeni girenler arasında farklar büyük olmayablir ancak 30 yıl çalışmış kişiler için sektöre ve pozisyona bağlı büyük farklılıklar olabilir. Eğer biz bunu önemsemeden analiz yaparsak varyans farklı seviyeler için çok farklı olacaktır. Bunu çözmek için farklı regresyon metodları veya log transformasyonu uygulayabiliriz. Başka örnek olarak gelir ve ev için harcamalar verilebilir. Sizce burada mümkün olan ne?
Residualarımızın dağılımı da normal olmalıdır. Bunu q-q plot ile görebiliriz. Çözüm olarak log transformasyonu yapabiliriz. Çok uç değerleri çıkarmamız da bir çözüm olabilir.
En tartışmaları konulardan birisi şimdi bahsedeceğim. Multicollinearity bağımsız değişkenlerimizin birbiri ile yüksek derecede ilişkili olması. Variance Inflation Factor testi ile ölçebiliriz. Eğer veriler çok ilişkili ise hangisinin etkisi olduğunu bilmemiz zorlaşır. Mesela Amerikan futbolunu düşünelim…
Ancak bazıları da eğer 2 bağımsız değişken arasında perfect ilişki yoksa bunun sorun olmaduğını söyler. Çünkü her açıklama daha iyi bir tahmine yardımcı olur ve bunları çıkarmak çok da gerekli değildir. Zaten homoscedasticity gibi analizi bias etmez sadece bağımsız değişkenlerin etkisini tam olarak anlamayı zorlaştırır. Ondan dolayı umursamadan devam edebilirsiniz diyenler de var. Kişinin hangi eğitimi aldığına göre değişiyor. Zaten iki değişken birbiri ile 1 korelasyona sahipse istatistik yazılımları değişkenleri drop eder.
Bu assumptionları karşılarsak modelimiz BLUE yani Best Linear Unbiased Estimator olur. Burada asıl önemli olan unbiased ve linear assumptionlarının bozulmaması. Diğerleri de iyi bir model için önemli olsa bile bunlar bozulursa modelimiz çöp olabilir.
##Korelasyon
Korelasyon ve Regresyon.
Correlation does not imply causation çok klasik bir laf. Beğenmesem ve hatalı anlaşılsa bile doğru bir noktaya da değiniyor. Korelasyondan nedensellik çıkartmak mümkün değil ve sadece iki değişken beraber hareket ediyor mu ona bakıyor. Regresyon ise bu ilişkinin yönünün yanısıra gücünü de söylüyor ve bu modelin ne kadar başarılı olduğunu da regresyon ile öğrenebiliriz. Regresyon yine tam olarak nedensellik belirtmese bile bizi o yönde ilerletir. Şimdi bir örnek ile regresyon R’da nasıl yapılır inceleyelim.
Regresyon içinse kodumuz:
Burada lm linear model anlamına geliyor. Analiz metodumuzun linear regresyon olduğunu bu şekilde belirtiyoruz. İlk başta dependent variable yani Y değerini alacak değişkenimizi yazıyoruz. Ondan sonra “tilda” adı verilen ~ işaretimizi kullanıyoruz. Bu işaret R’da bağımlı ve bağımsız değişkenlerimizi ayırmaya yarıyor. Ondan sonraysa bağımsız değişkenlerimiz gelecek. En son olaraksa data kısmındaysa verisetinin adını giriyoruz.
Linear Regresyonu Yorumlama
##
## Call:
## lm(formula = mpg ~ wt, data = mtcars)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.5432 -2.3647 -0.1252 1.4096 6.8727
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 37.2851 1.8776 19.858 < 2e-16 ***
## wt -5.3445 0.5591 -9.559 1.29e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.046 on 30 degrees of freedom
## Multiple R-squared: 0.7528, Adjusted R-squared: 0.7446
## F-statistic: 91.38 on 1 and 30 DF, p-value: 1.294e-10Modelimizin sonuçları bu şekilde. Şimdi yorumlama kısmına başlayalım. Bu analiz aracın galon başına kaç bil gidebileceğini aracın ağırlığına bağlı olarak tahmin etmeye çalışıyor. Mesela bu örnekte intercept bizim için mantılı değil çünkü aracın ağırlığının 0 olma ihtimali yok. Ondan dolayı intercept kısmına hiç bakmadan geçebiliriz. Weight ise istatistiksel olarak anlamlı çıkmış. Bunu da 0.001 seviyesinde reddetmeyi başarmışız. Yani aracın ağırlığının aracın yakıt tüketimine etkisi yoktur null hipotezini reddetmeyi başarıyoruz.
Buna ek olarak estimate dediğimiz kısım X değişkeninin Y üzerine olan etkisi. Bu da negatif görünüyor. Yani aracın ağırlığı arttıkça aracın 1 galon benzinle gidebileceği mesafe azalıyor diyebiliriz. Burdaki sayıyı yorumlamak içinse bağımsız değişkenimizin unitini bilmemiz gerekir. Yani 100 kilo mu 1 kilo mu hangi birimdeyse ona göre yorum yapabilirz. Bu örnekte ise ağırlık birimi 1000 pound. Yani bu bize aracın ağırlığında olan her 1000 poundlık artışın aracın galon başına gidebileceği mil mesafesini 5.34 azaltmış.
Aşağıda görebileceğiniz başka bir değer var. R-squared dediğimiz. Bu modelimizin bağımlı değişkendeki varyasyonu ne derece açıkladığına bakıyor. Modelimiz yaklaşık %75 bir açıklama gücüne sahip ki bu çok iyi bir oran. Ancak hala açıklayamadığımız %25 civarı bir varyasyon var.
Modelin p değeri de 0.001 altında. Bu da bize yine boş bir modele göre aracın ağırlığını eklemenin tahmin gücümüzü arttırdığını gösteriyor. Şimdiyse model ile alakalı görsellere bakıp çıkarımlar yapalım ve assumptionları inceleyelim.
Scatter plot ile:
plot(mtcars$wt, mtcars$mpg,
main = "Scatter plot with Regression Line",
xlab = "Weight (1000 lbs)",
ylab = "Miles per Gallon")
abline(model_1, col = "blue", lwd = 2)
Bütün plotları açıklamak gerekirse:
Residuals vs Fitted ise linearity assumptionı karşılanyıor mu diye bakar. Beklenen çizginin düz çizgiyi takip etmesidir.
Q-Q Plot Residualların dağılımına bakar. Eğer diagonel çizgi etrafında dağılıyorlarsa normal dağılıma sahipler deriz.
Scale-Location homoscedasticity gösterir. İyi bir modelde çizginin düz olmasını bekleriz.
Residuals vs Leverage ise modelin coefficentlarını etkileyebilecek verileri gösterir. Çizgilerin dışındaki yerlere bakmamız gerekir.
hist(residuals(model_1),
main = "Histogram of Residuals",
xlab = "Residuals",
col = "lightblue",
border = "white")
Normality Test
##
## Shapiro-Wilk normality test
##
## data: residuals(model_1)
## W = 0.94508, p-value = 0.1044Autocorrelation
## lag Autocorrelation D-W Statistic p-value
## 1 0.3628798 1.251727 0.014
## Alternative hypothesis: rho != 0Homoscedasticity Test
##
## studentized Breusch-Pagan test
##
## data: model_1
## BP = 0.040438, df = 1, p-value = 0.8406Şu ana kadar tek değişkenli analiz yaptık. Ancak gerçek bir analiz yaparken genellikle birden fazla değişken kullanamız gerekecek. Bunun da bir örneğini inceleyelim.
##
## Call:
## lm(formula = mpg ~ wt + am + hp, data = mtcars)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.4221 -1.7924 -0.3788 1.2249 5.5317
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 34.002875 2.642659 12.867 2.82e-13 ***
## wt -2.878575 0.904971 -3.181 0.003574 **
## am 2.083710 1.376420 1.514 0.141268
## hp -0.037479 0.009605 -3.902 0.000546 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.538 on 28 degrees of freedom
## Multiple R-squared: 0.8399, Adjusted R-squared: 0.8227
## F-statistic: 48.96 on 3 and 28 DF, p-value: 2.908e-11Burada bir öncekiyle aynı olarak mpg ve wt değerleri var. am ve hp ise automatic-manuel ile beygir değişkenleri olarak eklendi. Bu sonuçlara baktığımız zaman ağırlık değişkeni 0.01 seviyesinde beygir is 0.001 seviyesinde anlamlı çıktı. Yani ağırlığın etkisi yok ve beygirin etkisi yok null hipotezlerini reddetmeyi başardık. Vites tipi ise anlamlı bir etkiye sahip değil. Ondan dolayı bunu yorumlamıyoruz. Hatta estimate değişkenlerinin yönünü ve gücünü yorumlarken de onu geçebiliriz çünkü bizim için anlamlı değil.
Bu sefer ağırlığın etkisi -2.88. Yani diğer değişkenler sabit tutulduğunda aracın ağırlığında yaşanan her 1 birimlik değişim mgp değişkenini 2.88 azaltmış. Beygirde ise diğer değişkenler sabit tutulduğunda her 1 birimlik artış mpg’yi 0.04 azaltmış. Bu da beygirin ve ağırlığın birimlerini düşündüğümüz zaman normal bir sayı.
RSE değeri yine model fitini ölçmek için kullanabileceğimiz bir değer. Bu tahminlerimizin ortalama olarak 2.54 mpg sapmaya sahip olduğunu gösteriyor. mpg değerinin ortalamasını ve range değerlerini de hesaba katınca %10 civarı bir sapma oluyor. Uluslararası bir ölümü yok ancak veriye ve alana göre değişir.
Ek olarak \(R^2\) değeri de 0.83 civarı. Bu hem çok yüksek bir oran hem de modelimizin iyi bir fite sahip olduğunu gösteriyor.
Burada wt’nin değeri azaldı. Bu bize sadece wt’nin olduğu model ile açıkladığımız değişimin aslında hp ve am ile de açıklanabileceğini gösteriyor. Yani ilk modelimiz wt’nin yanısıra hp ve am değişkenlerinin de confounding dediğimiz etkisi ölçüyordu. Bunları ayrı ölçtüğümüz zaman wt’nin etkisi azalmış oldu çünkü artık önceden omit ettiğimiz değişkenlerin de etkisini ölçmeye başlamış olduk.
## Likelihood ratio test
##
## Model 1: mpg ~ wt
## Model 2: mpg ~ wt + am + hp
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 3 -80.015
## 2 5 -73.067 2 13.895 0.0009612 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Likelihood Ratio Test ise iki modeli karışılaştırmamıza yarıyor. Burada dikkat edeceğimiz nokta p değeri ve Log Likelihood Ratio. Daha büyük bir LogLik değeri modelin daha iyi bir fit’e sahip olduğunu gösterir. P değerinin de anlamlı olması bize ikinci modelin ilk modele göre daha anlamlı ve yardımcı olduğunu gösteriyor.
Şimdi bu modelin de diagnostic değerlerine bakalım.
hist(residuals(model_2),
main = "Histogram of Residuals",
xlab = "Residuals",
col = "lightblue",
border = "white")
Autocorrealation
## lag Autocorrelation D-W Statistic p-value
## 1 0.232331 1.43292 0.038
## Alternative hypothesis: rho != 0Homoscedasticity
##
## studentized Breusch-Pagan test
##
## data: model_2
## BP = 5.534, df = 3, p-value = 0.1366Son olaraksa multiple regression modeline özel olarak vif testi var. Genelde 5 veya 10 üstü değerler multicollinearity göstergesi denir. Ancak dediğim gibi modern analizde çok sorun teşkil etmediği söyleniyor.
## wt am hp
## 3.774838 2.271082 2.088124Eğer verimizde outlier değerler varsa da log transformation kullanarak bunların yaratabileceği sorunlardan kurtulabiliriz.
Mesela hp değerlerine bakalım:
hist(mtcars$hp,
main = "Original Distribution of Horsepower (hp)",
xlab = "Horsepower (hp)",
col = "lightblue",
breaks = 10)
boxplot(mtcars$hp,
main = "Boxplot of Original Horsepower (hp)",
ylab = "Horsepower (hp)",
col = "lightblue")
Özellikle yukarıda bir uç değer görüyoruz. Verimiz küçük olduğu için tek nokta ama büyük verilerde daha fazla olacaktır.
par(mfrow = c(1, 2))
hist(mtcars$hp,
main = "Original Distribution of Horsepower (hp)",
xlab = "Horsepower (hp)",
col = "lightblue",
breaks = 10)
hist(mtcars$log_hp,
main = "Log-Transformed Distribution
of Horsepower (log(hp))",
xlab = "Log(Horsepower)",
col = "lightgreen",
breaks = 10)
par(mfrow = c(1, 2))
boxplot(mtcars$hp,
main = "Boxplot of Original Horsepower (hp)",
ylab = "Horsepower (hp)",
col = "lightblue")
boxplot(mtcars$log_hp,
main = "Boxplot of Log-Transformed Horsepower (log(hp))",
ylab = "Log(Horsepower)",
col = "lightgreen")
Bu grafiklerde de gördüğümüz üzere log transformasyonu sonrasında veri dağılımı daha normal oldu ve outlier değer boxplottan kayboldu. Bu durumların büyük sorun yaşatabileceği durumlar olduysa bunlardan kurtulup rahat ve daha iyi bir analiz sonucu elde etmemize yarar. Böylece linear regresyon haftasını bitirmiş olduk. Haftaya bağımlı değişkenimizin continious değil de 0 1 veya 0 1 2 gibi değerler alabildiği logistic regresyonu öğreneceğiz. Böylece de R ile yapılan temel analizleri öğrenmiş olacağız.
Dummy variable ise 2 farklı değer alabilen değişkenlere verilen bir isimdir. Mesela kadın-erkek veya beyaz-siyah gibi değişkenleri dummy variable olarak alabiliriz. Bunları analize kattığımız durumda birinci kategorinin ikinciye göre farkını gösterir. Örnekle bakacak olursak:
##
## Call:
## lm(formula = mpg ~ am, data = mtcars)
##
## Residuals:
## Min 1Q Median 3Q Max
## -9.3923 -3.0923 -0.2974 3.2439 9.5077
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 17.147 1.125 15.247 1.13e-15 ***
## am 7.245 1.764 4.106 0.000285 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.902 on 30 degrees of freedom
## Multiple R-squared: 0.3598, Adjusted R-squared: 0.3385
## F-statistic: 16.86 on 1 and 30 DF, p-value: 0.000285Bu bize iki grup arasındaki farkın 7.245 olduğunu gösteriyor. Bunu aslında geçen hafta da öğrenmiştik.
##
## Welch Two Sample t-test
##
## data: mpg by am
## t = -3.7671, df = 18.332, p-value = 0.001374
## alternative hypothesis: true difference in means between group 0 and group 1 is not equal to 0
## 95 percent confidence interval:
## -11.280194 -3.209684
## sample estimates:
## mean in group 0 mean in group 1
## 17.14737 24.39231difference_in_means <- t_test_result$estimate[2] - t_test_result$estimate[1]
print(difference_in_means)## mean in group 1
## 7.244939Ve birebir aynı sonuçlar karşımıza çıkıyor. Bu da t-testin aslında bir regresyon varyasyanı veya daha doğrusu bir linear model olduğunun göstergesi. Eğer iki ortak grubu karşılaştırıyorsak bunu t-test yerine regresyon ile de yapabiliriz.
Simpson’s Paradox
Simpson paradoksu, bir veri kümesinin tamamında görülen bir ilişkinin, bu veri kümesi alt gruplara bölündüğünde her bir alt grupta farklı, hatta tam tersi bir ilişki göstermesi durumudur. Yani, bütün olarak bakıldığında bir trend görülürken, daha detaylı incelendiğinde bu trendin tam tersinin ortaya çıkması olayıdır.
Paradoksun ortaya çıkmasının nedeni genellikle veri setindeki gizli bir faktördür. Bu gizli faktör, genel trendi alt gruplardaki trendlerden farklı kılar. Bu paradoks, veri analizinde dikkatli olmamız gerektiğini gösterir. Sadece genel eğilimlere bakmak yanıltıcı olabilir; veriyi daha detaylı incelemek ve olası alt grupları dikkate almak önemlidir.
set.seed(123)
n_total <- 500
neighborhoods <-
rep(c("A", "B", "C"),
c(200, 150, 150))
base_sizes <-
c(1800, 1600, 1400)[match(neighborhoods,
c("A", "B", "C"))]
base_prices <-
c(280000, 320000, 360000)[match(neighborhoods,
c("A", "B", "C"))]
price_slopes <-
c(100, 120, 140)[match(neighborhoods,
c("A", "B", "C"))]
sizes <- base_sizes + rnorm(n_total, 0, 300)
prices <- base_prices +
price_slopes * (sizes - base_sizes) +
rnorm(n_total, 0, 15000) -
80 * (sizes - mean(sizes))
house_data <- data.frame(
Neighborhood = neighborhoods,
Size = sizes,
Price = prices
)## [1] -0.2649864library(tidyverse)
ggplot(house_data, aes(x = Size, y = Price)) +
geom_point(alpha = 0.7) +
geom_smooth(method = "lm", se = FALSE) +
geom_smooth(aes(color = NULL), method = "lm", se = FALSE,
color = "black", linetype = "dashed") +
theme_minimal() +
labs(title = "House Prices vs Size: Simpson's Paradox",
x = "House Size (sq ft)", y = "Price ($)")
neighborhood_cors <-
sapply(unique(house_data$Neighborhood),
function(nbhd) {
subset <- house_data$Neighborhood == nbhd
cor(house_data$Size[subset],
house_data$Price[subset])
})
neighborhood_cors## A B C
## 0.3754936 0.5501883 0.7347689ggplot(house_data, aes(x = Size,
y = Price,
color = Neighborhood)) +
geom_point(alpha = 0.7) +
geom_smooth(method = "lm",
se = FALSE) +
geom_smooth(aes(color = NULL),
method = "lm", se = FALSE,
color = "black",
linetype = "dashed") +
theme_minimal() +
labs(title =
"House Prices vs Size: Simpson's Paradox",
x = "House Size (sq ft)",
y = "Price ($)")## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'