Öğrenme Günlüğü 4/Hafta 5
dilek vural
2026-03-22
library(moderndive)
library(ggplot2)
library(dplyr)
library(lavaan)
library(semPlot)Ders Notları ve Kodlar
Standartlaştırılmış artık varyansı: bir içsel değişkendeki yordayıcılar tarafından açıklanamayan varyans miktarıdır.
Form değişkeninde disturbance değerini 1261.54 bulmuştuk, kendi varyansı da 1508.781. 1261.54 / 1508.781 = 0.836 = form değişkende açıklanamayan varyans yaklaşık %84’tür.
library(readr)
veri <- read_table("import/illness.dat", col_names = FALSE)
colnames(veri) <- c("form", "stres", "hastalik", "egzersiz", "dayaniklilik")yol_model_v2 <-
'stres ~ dayaniklilik
hastalik ~ form + stres
form ~ egzersiz + dayaniklilik
stres ~ form
egzersiz ~~ dayaniklilik'
yol_fit_v2 <- sem(yol_model_v2, veri)fitted(yol_fit_v2)$cov## stres hastlk form dynkll egzrsz
## stres 4914.301
## hastalik 1788.526 3800.642
## form -626.495 -766.862 1508.781
## dayaniklilik -833.343 -354.767 242.743 1496.498
## egzersiz -341.945 -488.496 1008.959 9.118 4883.673# aynı sonuç için diğer kod= inspectSampleCov(yol_fit,data=veri)Elde edilen değerler farklı ölçek düzeylerinde olduğu için, yorumlamayı standartlaştırmak adına değerlerin z puanı cinsinden bir puana dönüştürürüz.
parameterEstimates(yol_fit_v2,standardized = TRUE)inspect(yol_fit_v2, "rsquare") ## stres hastalik form
## 0.128 0.225 0.164# std.all sütünündaki en büyük katsayı 0.371 ile form ~ egzersiz arasında.
# Modeldeki en büyük etki egzersizin forma olan etkisi diyebiliriz.
# Her bir içsel değişken için 𝑅2 değeri ve standartlaştırılmış artık varyansının toplamı “1”e eşit olması gerektiğinden bu çıktıda forma ait olan standartlaştırılmış artık varyans 0.836 olduğundan 𝑅2 değeri 1- 0.836 =0.164.Form değişkenindeki varyansın yaklaşık %16’sının yordayıcılar tarafından açıklandığını söyleyebiliriz.Model Sonuçlarının Rapor Edilmesi
library(knitr)
standardizedsolution(yol_fit_v2) %>%
filter(op == "~") %>%
select('Bağımlı Değişkenler'=lhs, Gosterge=rhs,
B=est.std, SE=se, Z=z, 'p-value'=pvalue) %>%
knitr::kable(digits = 3, booktabs=TRUE, format="markdown",
caption="Factor Loadings")| Bağımlı Değişkenler | Gosterge | B | SE | Z | p-value |
|---|---|---|---|---|---|
| stres | dayaniklilik | -0.277 | 0.046 | -6.085 | 0 |
| hastalik | form | -0.238 | 0.044 | -5.375 | 0 |
| hastalik | stres | 0.359 | 0.043 | 8.434 | 0 |
| form | egzersiz | 0.371 | 0.043 | 8.721 | 0 |
| form | dayaniklilik | 0.160 | 0.045 | 3.544 | 0 |
| stres | form | -0.185 | 0.047 | -3.977 | 0 |
# YORUMLAR
# Tablodaki en yüksek pozitif katsayı 0.371 ile egzersiz ve form arasındadır.
# Stres ve Hastalık (0.359): Pozitif bir etki. Yani stres arttıkça hastalık bulguları da artıyor. Bu, modeldeki ikinci en güçlü etkidir.
# Form ve Hastalık (-0.238): Negatif bir etki. Form düzeyi arttıkça hastalık düzeyi azalıyor.
# Dayanıklılık ve Stres (-0.277): Dayanıklılık, stresi doğrudan ve anlamlı bir şekilde azaltıyor (Z = -6.085).
# Form ve Stres (-0.185): Fiziksel formda olmak da stresi azaltıyor.Model karşılaştırma
devtools::install_github("dr-JT/semoutput") # Önceden kurduğumuz model 1
yol_model_v1 <-
'stres ~ egzersiz + dayaniklilik
hastalik ~ egzersiz + dayaniklilik + form + stres
form ~ egzersiz + dayaniklilik
stres ~ form
egzersiz ~~ dayaniklilik'
yol_fit_v1 <- sem(yol_model_v1, veri)library(semoutput)
sem_anova(yol_fit_v2,yol_fit_v1)| term | df | AIC | BIC | statistic | Chisq.diff | RMSEA | Df.diff | p.value |
|---|---|---|---|---|---|---|---|---|
| 1 | 3 | 21418.20 | 21466.10 | 1.354 | 1.354 | 0 | 3 | 0.716 |
| 2 | 0 | 21422.85 | 21482.72 | 0.000 | NA | NA | NA | NA |
# ikinci alternatif kod= sem_modelcomp(yol_fit_v2,yol_fit_v1)
# yol_fit_v2 seçmeliyim çünkü; serbestlik derecesi (df) 3, diğer modelin sıfır, AIC ve BIC değerleri daha düşük ve p = 0.716, ve anlamlı bir fark yaratmıyor. Doğrudan, Dolaylı ve Toplam Etkiler
Toplam etki: bir değişken bir birim değiştiğinde diğer bir değişkenin ne kadar değişeceğini belirtir.
Toplam etkinin iki kısmı var; doğrudan ve dolaylı etki
Bir değişkenin diğer bir değişken üzerindeki doğrudan etkisi yol modelindeki ağırlığıyla belirtilir.
Dolaylı etkiler doğrudan etkilerin çarpımları olarak istatistiksel olarak kestirilir.
yol_model <- '
# Regresyon Yolları
stres ~ s_d*dayaniklilik + s_f*form
hastalik ~ h_f*form + h_s*stres
form ~ f_e*egzersiz + f_d*dayaniklilik
# Kovaryans
egzersiz ~~ dayaniklilik
# Dolaylı Etkiler (Mediation)
ind_egz_form_hast := f_e * h_f # Egzersiz -> Form -> Hastalık
ind_day_stres_hast := s_d * h_s # Dayanıklılık -> Stres -> Hastalık
ind_day_form_hast := f_d * h_f # Dayanıklılık -> Form -> Hastalık
# Toplam Dolaylı Etki
total_indirect := ind_egz_form_hast + ind_day_stres_hast + ind_day_form_hast
# Total Effects
tot_egz_hast := f_e*h_f + f_e*s_f*h_s
tot_day_hast := s_d*h_s + f_d*h_f
'
fsem1 <- sem(yol_model,veri)standardizedsolution(fsem1)library(knitr)
standardizedsolution(fsem1) %>%
filter(op %in% c("~", ":=", "r2")) %>%
select(Hedef=lhs,
Kaynak=rhs,
Islem = op,
B=est.std,
SE=se,
Z=z,
'p-value'=pvalue) %>%
knitr::kable(digits = 3, booktabs=TRUE, format="markdown",
caption="Yol Analizi (SEM) Final Sonuç Tablosu")| Hedef | Kaynak | Islem | B | SE | Z | p-value |
|---|---|---|---|---|---|---|
| stres | dayaniklilik | ~ | -0.277 | 0.046 | -6.085 | 0.000 |
| stres | form | ~ | -0.185 | 0.047 | -3.977 | 0.000 |
| hastalik | form | ~ | -0.238 | 0.044 | -5.375 | 0.000 |
| hastalik | stres | ~ | 0.359 | 0.043 | 8.434 | 0.000 |
| form | egzersiz | ~ | 0.371 | 0.043 | 8.721 | 0.000 |
| form | dayaniklilik | ~ | 0.160 | 0.045 | 3.544 | 0.000 |
| ind_egz_form_hast | f_e*h_f | := | -0.088 | 0.020 | -4.489 | 0.000 |
| ind_day_stres_hast | s_d*h_s | := | -0.100 | 0.021 | -4.824 | 0.000 |
| ind_day_form_hast | f_d*h_f | := | -0.038 | 0.013 | -2.949 | 0.003 |
| total_indirect | ind_egz_form_hast+ind_day_stres_hast+ind_day_form_hast | := | -0.226 | 0.031 | -7.189 | 0.000 |
| tot_egz_hast | f_eh_f+f_es_f*h_s | := | -0.113 | 0.022 | -5.203 | 0.000 |
| tot_day_hast | s_dh_s+f_dh_f | := | -0.138 | 0.023 | -5.960 | 0.000 |
Sonuç: 1.Doğrudan Etkiler Form → Hastalık (-0.238): Fiziksel formun arttıkça hastalık riskin anlamlı şekilde düşüyor. Stres → Hastalık (0.359): Stres arttıkça hastalık belirtileri de artıyor. Egzersiz → Form (0.371): Egzersizin fiziksel form üzerindeki etkisi pozitif.
2.Aracılık (Dolaylı) Etkileri: (ind_egz_form_hast = -0.088): Egzersiz yapmak, formu artırarak hastalığı dolaylı yoldan azaltıyor. (ind_day_stres_hast = -0.100): Psikolojik dayanıklılık, stresi azaltarak hastalık riskini düşürüyor.
3.Toplam Dolaylı Etki (-0.226): total_indirect etki anlamlı (p=0.000): Hem dayanıklılığın stres üzerindeki etkisi hem de egzersizin form üzerindeki etkisi birleşerek, hastalık riskinde toplamda -0.226 birim azalmaya sebep oluyor.
4.Egzersizin Toplam Gücü (tot_egz_hast = -0.113, p=0.000): Egzersizin hastalık üzerinde doğrudan bir etkisi olmasa da, hem form hem de stres üzerinden yaptığı tüm etkilerin toplamı, hastalıkta -0.113 birimlik bir azalma sağlıyor.
5.Dayanıklılığın Toplam Gücü (tot_day_hast = -0.138, p=0.000): Psikolojik dayanıklılık, hastalıkta egzersizden biraz daha yüksek (-0.138) bir düşüşe sebep oluyor. Bu fark, dayanıklılığın hem stresi doğrudan düşürmesinden hem de dolaylı yoldan formu desteklemesinden kaynaklanıyor.
Aracılık Etkisi
library(GGally)
library(lavaan)
library(broom)
library(DT)
library(dplyr)
library(knitr)
library(semPlot)
library(semptools)
library(semoutput)
library(stargazer)
library(multilevel)
library(bda)
library(mediation)
library(gvlma)
library(rockchalk)Bozucu Değişken: Confounding: Bağımsız ve bağımlı değişkenler arasında yanlış bir ilişki oluşturabilen üçüncü bir değişkendir, gerçek ilişkinin yanlış değerlendirilmesine yol açabilir.
Bastırıcı Değişken (suppressor): regresyon analizinde bağımsız değişkenlerden biriyle yüksek korelasyon gösteren ancak bağımlı değişkenle doğrudan ilişki göstermeyen bir değişkendir.
Aracı Değişken (mediator) ve Nedensellik: Bağımsız değişken ile bağımlı değişken arasındaki ilişkiyi açıklayan bir değişkendir.
Varsayımların İncelenmesi
araci <- readRDS("import/araci.Rds")
summary(araci)## X M Y
## Min. :158.8 Min. :108.0 Min. :38.55
## 1st Qu.:171.5 1st Qu.:117.7 1st Qu.:46.05
## Median :175.4 Median :122.7 Median :49.57
## Mean :175.6 Mean :122.4 Mean :49.56
## 3rd Qu.:179.8 3rd Qu.:127.3 3rd Qu.:52.39
## Max. :190.3 Max. :136.8 Max. :63.82# veri setinde yanlış girilmiş bir değer olup olmadığının kontrolülibrary(naniar)
araci <- readRDS("import/araci.Rds")
any_na(araci) # kayıp veri kontrolü## [1] FALSEn_miss(araci) # toplam eksik veri sayısı## [1] 0prop_miss(araci) # eksik veri oranı## [1] 0miss_var_summary(araci) # eksik veri tablosulibrary(car)
model<-lm(Y~X+M, data=araci)
durbinWatsonTest(model)## lag Autocorrelation D-W Statistic p-value
## 1 0.05649613 1.829628 0.38
## Alternative hypothesis: rho != 0# DW testi bağımsızlık varsayımı için yapılır, hata terimlerinin ilişkisiz olduğunu belirlemenin bir yoludur. 1.5-2.5 arası ise hatalr ilişkisiz diyebiliriz. dogrusallik <- plot(model,1)
# Doğrusallık varsayımı: Değişkenler arasında doğrusal bir ilişki olmalıdır. değişken sayısı fazla olduğunda artıkların grafiği incelenir.
# Kırmızı çizginin yatay eksene paralelliği, veriye uygun olduğunu kanıtlar.# Uç değerler kontrolü
ucdegerbsiz<-boxplot(araci$X, ylab= "saat")
ucdegerbli<-boxplot(araci$Y, ylab="uyanıklık derecesi")
ucdegeraraci<-boxplot(araci$M, ylab="kahve tüketimi")
# Tek değişkenli uç değer kontrolü için z puan hesaplama
library(outliers)
z.scores<-scores(araci, type="z")
round(head(z.scores),2)# Çok değişkenli uç değerleri belirlemek için Mahalanobis uzaklığını hesaplanması:
md <- mahalanobis(araci, center = colMeans(araci), cov = cov(araci))
alpha <- .001
cutoff <- (qchisq(p = 1 - alpha, df = ncol(araci))) # kritik değer belirleme
ucdegerler<-which(md>cutoff)
ucdegerler## integer(0)# Normallik : Hata terimlerinin normal dağılıp dağılmadığı incelenir.
plot(model,2)
# ek olarak hist(model$residuals, col='steelblue') bakılabilir.
# Noktalar referans çizgiden sapmadığı için veri setinde hata terimlerinin dağılımı normaldir.# Eş varyanslılık :Varyansların homojenliği (homoscedasticity), bağımlı değişken(ler)in bağımsız değişken(ler)in aralığı boyunca aynı düzeyde varyansa sahip olduğu sayıltısıdır.
# Kırmızı çizgi biraz saptığı için yeterli kanıt yok.
plot(model,3)
library(MVN) # çok değişkenli normalliğin incelenmesi
sonuc <- mvn(data = araci, mvn_test = "mardia")
sonuc$multivariate_normalitycoklubag<-cor(araci)
round(coklubag,2)## X M Y
## X 1.00 0.66 0.21
## M 0.66 1.00 0.43
## Y 0.21 0.43 1.00# Çoklu bağlantılılık : Bağımsız değişkenlerin birbiriyle çok yüksek korelasyona sahip olması durumudur..80 üzerinde değer istemeyiz. library(olsrr)
ols_vif_tol(model)# tolerance depğeri, 0.10'dan küçük olmamalı
# vif değeri: 10'dan büyük olmayacakFarklı Modeller Aracı değişken sayısana göre tekli veya çoklu olabilir.
Eğer bağımlı ve bağımsız değişken arasındaki ilişkiyi tek bir aracı değişken açıklamaya çalışıyorsa tekli aracılık modeli.
Bağımlı ve bağımsız değişken arasındaki ilişkiyi etkileyen birden fazla aracı değişken olması ve bu aracı değişkenlerin nedensel olarak birbirlerini etkilemediği durumda paralel çoklu aracılık modeli.bu iki aracı değişkenler ilişkili olabilir ama nedensel bir ilişki yoktur.
En az iki aracı değişken arasında nedensel bir ilişki olduğu varsayılırsa: Çoklu Aracılık Modeli-Serial.
Tekli Aracılık Modeli Analizi Aracılık analizi, bağımsız değişkenin bağımlı değişken üzerindeki etkisinin üçüncü bir değişken olan aracı üzerinden işleyip işlemediğini test eder.
Yöntem 1: Baron & Kenny:
x:şafaktan itibaren geçen saat sayısı m:kahve tüketimi y:öznel uyanıklık derecelerini
4 adet adımı var:
- adım: X ile Y arasındaki ilişkinin kestirilmesi(ilişki olmalıdır)
- adım: X ile M arasındaki ilişkinin kestirilmesi (ilişki olmalıdır)
- adım: X’i kontrol ederek Y ile M arasındaki ilişkinin kestirilmesi (y ve m arasında ilişki olmalıdır)
- adım: M’yi kontrol ederek Y ile X arasındaki ilişkinin kestirilmesi(Anlamlı olmamalı ve yaklaşık 0 olmalıdır.)
library(stargazer)
#1. toplam etki
fit <- lm(Y ~ X, data=araci)
#2. YOL A (X on M)
fita <- lm(M ~ X, data=araci)
#3. YOL B (M on Y, X kontrol edildiğinde )
fitb <- lm(Y ~ M + X, data=araci)
#4. YOL C - ters yol (Y on X, M kontrol edildiğinde)
fitc <- lm(X ~ Y + M, data=araci)
# Özet tablo
stargazer(fit, fita, fitb, fitc, type = "text", title = "Baron and Kenny Yöntemi ile Aracılık Etkisi",digits = 2,
font.size ="tiny")##
## Baron and Kenny Yöntemi ile Aracılık Etkisi
## =========================================================================================================
## Dependent variable:
## -------------------------------------------------------------------------------------
## Y M Y X
## (1) (2) (3) (4)
## ---------------------------------------------------------------------------------------------------------
## Y -0.11
## (0.10)
##
## M 0.42*** 0.70***
## (0.10) (0.08)
##
## X 0.17** 0.66*** -0.11
## (0.08) (0.08) (0.10)
##
## Constant 19.88 6.04 17.32 96.11***
## (14.26) (13.42) (13.16) (9.28)
##
## ---------------------------------------------------------------------------------------------------------
## Observations 100 100 100 100
## R2 0.04 0.43 0.19 0.44
## Adjusted R2 0.03 0.43 0.18 0.43
## Residual Std. Error 5.16 (df = 98) 4.85 (df = 98) 4.76 (df = 97) 4.82 (df = 97)
## F Statistic 4.34** (df = 1; 98) 75.31*** (df = 1; 98) 11.72*** (df = 2; 97) 38.39*** (df = 2; 97)
## =========================================================================================================
## Note: *p<0.1; **p<0.05; ***p<0.01# Yorumlar:
# 1. (X) ile uyanıklık (Y) arasında anlamlı bir pozitif ilişki var.
# 2. şafaktan beri geçen saatlerin (X) kahve tüketimi (M) ile de pozitif ilişkisi olduğunu göstermektedir.
# 3. kahve tüketiminin (M) şafaktan bu yana geçen saat (X) kontrol edildiğinde uyanıklığı (Y) pozitif olarak yordadığını göstermektedir.
# 4. kahve tüketimi (M) kontrol edildiğinde uyanıklık (Y) şafaktan bu yana geçen süreyi (X) yordamamaktadır.
# Kahve tüketimi kontrol edildiğinde şafaktan bu yana geçen saatler ile uyanıklık arasındaki ilişki artık anlamlı olmadığından, bu durum kahve tüketiminin aslında bu ilişkiye aracılık ettiğini göstermektedir. Bulunan dolaylı etkinin anlamlılığını test etmek için iki temel yöntem var: Sobel testi ve bootstrapping.#Sobel Test
library(multilevel)
sobel(araci$X, araci$M, araci$Y)## $`Mod1: Y~X`
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 19.8836805 14.2637142 1.394004 0.16646905
## pred 0.1689931 0.0811601 2.082220 0.03992761
##
## $`Mod2: Y~X+M`
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 17.3217682 13.16215851 1.316028 1.912663e-01
## pred -0.1117904 0.09949262 -1.123605 2.639537e-01
## med 0.4238113 0.09899469 4.281152 4.371472e-05
##
## $`Mod3: M~X`
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.0449365 13.41692114 0.4505457 6.533122e-01
## pred 0.6625203 0.07634187 8.6783345 8.871741e-14
##
## $Indirect.Effect
## [1] 0.2807836
##
## $SE
## [1] 0.07313234
##
## $z.value
## [1] 3.83939
##
## $N
## [1] 100library(bda)
attach(araci)
mediation.test(M,X,Y)# yorum: şafaktan bu yana geçen saatler ile uyanıklık hissi arasındaki ilişkiye kahve tüketimi önemli ölçüde aracılık etmektedir (z’ = 3.84, p < .001).YÖNTEM 2:
# Analiz öncesi varsayım kontrolü
library(mediation)
library(gvlma)
fitM <- lm(M ~ X, data=araci) # şafaktan beri geçen saatler (X) kahve tüketimini (M) etkiler
fitY <- lm(Y ~ X + M, data=araci) # kahve tüketimi kontrol edilirken, şafaktan bu yana geçen saat uyanıklığı etkiler
gvlma(fitM) # veri pozitif çarpıktır##
## Call:
## lm(formula = M ~ X, data = araci)
##
## Coefficients:
## (Intercept) X
## 6.0449 0.6625
##
##
## ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
## USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
## Level of Significance = 0.05
##
## Call:
## gvlma(x = fitM)
##
## Value p-value Decision
## Global Stat 8.833 0.06542 Assumptions acceptable.
## Skewness 6.314 0.01198 Assumptions NOT satisfied!
## Kurtosis 1.219 0.26949 Assumptions acceptable.
## Link Function 1.076 0.29959 Assumptions acceptable.
## Heteroscedasticity 0.223 0.63674 Assumptions acceptable.gvlma(fitY)##
## Call:
## lm(formula = Y ~ X + M, data = araci)
##
## Coefficients:
## (Intercept) X M
## 17.3218 -0.1118 0.4238
##
##
## ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
## USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
## Level of Significance = 0.05
##
## Call:
## gvlma(x = fitY)
##
## Value p-value Decision
## Global Stat 3.41844 0.4904 Assumptions acceptable.
## Skewness 1.85648 0.1730 Assumptions acceptable.
## Kurtosis 0.77788 0.3778 Assumptions acceptable.
## Link Function 0.71512 0.3977 Assumptions acceptable.
## Heteroscedasticity 0.06896 0.7929 Assumptions acceptable.# Çıktıdaki değerler:
# Global Stat: bağımlı değişken ile bağımsız değişken arasındaki ilişkinin doğrusal olup olmadığını kontrol eder.
# Çarpıklık (Skewness) ve basıklık (Kurtosis) varsayımları, artıkların (hataların) normal dağılıp dağılmadığını gösterir.
# Link Function: bağımlı değişkenin sürekli mi yoksa kategorik mi olduğunu kontrol eder.
# Heteroskedasticity (değişen varyans) varsayımı, artık varyansın değer aralığı boyunca nispeten sabit olduğu boşluğu test eder.plot(gvlma(fitM),1)
# modelin istatistiksel varsayımları (normallik, doğrusallık, uç değerler vb.) karşılayıp karşılamadığını gösteren diagnostik (tanısal) grafikleri elde ederiz.fitMed <- mediate(fitM, fitY, treat="X", mediator="M")
summary(fitMed)##
## Causal Mediation Analysis
##
## Quasi-Bayesian Confidence Intervals
##
## Estimate 95% CI Lower 95% CI Upper p-value
## ACME 0.2785581 0.1542942 0.4322466 <2e-16 ***
## ADE -0.1113693 -0.3181451 0.0765751 0.274
## Total Effect 0.1671888 0.0042728 0.3271031 0.046 *
## Prop. Mediated 1.5625302 0.5124882 10.1156280 0.046 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Sample Size Used: 100
##
##
## Simulations: 1000# ACME: (toplam etki - doğrudan etki)
# Ortalama Doğrudan Etki (ADE)
# Birleştirilmiş Dolaylı ve Doğrudan Etki (Toplam Etki)
# Prop. Mediated: ACME / Total Effect
# Kurulan model, kahve tüketiminin şafaktan bu yana geçen saatler ile uyanıklık hissi arasındaki ilişki üzerinde önemli bir etkisi olduğunu (ACME = .28, p < .001), şafaktan bu yana geçen saatlerin doğrudan etkisi olmadığını (ADE = -0.11, p = .27) ve toplam etkinin önemli olduğunu (.170, p < .05) göstermektedir.#Bootstrap yöntemi
fitMedBoot <- mediate(fitM, fitY, boot=TRUE, sims=999, treat="X", mediator="M")
summary(fitMedBoot)##
## Causal Mediation Analysis
##
## Nonparametric Bootstrap Confidence Intervals with the Percentile Method
##
## Estimate 95% CI Lower 95% CI Upper p-value
## ACME 0.280784 0.146694 0.421051 < 2e-16 ***
## ADE -0.111790 -0.281894 0.104693 0.28629
## Total Effect 0.168993 -0.002348 0.339581 0.05205 .
## Prop. Mediated 1.661509 -1.356902 10.179730 0.05205 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Sample Size Used: 100
##
##
## Simulations: 999# Kurulan model, kahve tüketiminin şafaktan bu yana geçen saatler ile uyanıklık hissi arasındaki ilişki üzerinde önemli bir etkisi olduğunu (ACME = .28, p < .001), şafaktan bu yana geçen saatlerin doğrudan etkisi olmadığını (ADE = -0.11, p = .27) ve toplam etkinin istatistiksel olarak anlamlı olmadığını (1.661, p > .05) göstermektedir.CHAPTER 15-Multiple Regression by Howell
##Çoklu Regresyon
library(haven)
library(dplyr)
library(knitr)
datahw <- read_sav("howell/Tab15-1.sav")
View(datahw)
psych::describe(datahw) %>% kable(digit=3)| vars | n | mean | sd | median | trimmed | mad | min | max | range | skew | kurtosis | se | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| id | 1 | 50 | 25.500 | 14.577 | 25.500 | 25.500 | 18.532 | 1.000 | 50.000 | 49.000 | 0.000 | -1.272 | 2.062 |
| State* | 2 | 50 | 25.500 | 14.577 | 25.500 | 25.500 | 18.532 | 1.000 | 50.000 | 49.000 | 0.000 | -1.272 | 2.062 |
| Expend | 3 | 50 | 5.905 | 1.363 | 5.768 | 5.747 | 1.259 | 3.656 | 9.774 | 6.118 | 1.041 | 0.877 | 0.193 |
| PTratio | 4 | 50 | 16.858 | 2.266 | 16.600 | 16.605 | 1.927 | 13.800 | 24.300 | 10.500 | 1.255 | 2.008 | 0.321 |
| Salary | 5 | 50 | 34.829 | 5.941 | 33.288 | 34.287 | 5.203 | 25.994 | 50.045 | 24.051 | 0.712 | -0.207 | 0.840 |
| PctSAT | 6 | 50 | 35.240 | 26.762 | 28.000 | 34.100 | 29.652 | 4.000 | 81.000 | 77.000 | 0.247 | -1.642 | 3.785 |
| Verbal | 7 | 50 | 457.140 | 35.176 | 448.000 | 456.550 | 42.995 | 401.000 | 516.000 | 115.000 | 0.120 | -1.483 | 4.975 |
| Math | 8 | 50 | 508.780 | 40.205 | 497.500 | 507.025 | 43.737 | 443.000 | 592.000 | 149.000 | 0.310 | -1.172 | 5.686 |
| SATcombined | 9 | 50 | 965.920 | 74.821 | 945.500 | 963.275 | 84.508 | 844.000 | 1107.000 | 263.000 | 0.222 | -1.367 | 10.581 |
| PctACT | 10 | 50 | 40.520 | 28.351 | 47.000 | 40.475 | 39.289 | 2.000 | 83.000 | 81.000 | -0.081 | -1.725 | 4.009 |
| ACTcombined | 11 | 50 | 21.052 | 0.881 | 21.300 | 21.130 | 0.815 | 18.700 | 22.400 | 3.700 | -0.765 | 0.056 | 0.125 |
| LogPctSAT | 12 | 50 | 3.157 | 0.995 | 3.332 | 3.217 | 1.338 | 1.386 | 4.394 | 3.008 | -0.254 | -1.514 | 0.141 |
# İlgili değişkenleri seçip korelasyon matrisi oluşturma
# Expend, PTratio, Salary, PctSAT, SAT, LogPctSAT
cor_data <- datahw %>%
dplyr::select(Expend, PTratio, Salary, PctSAT, SATcombined, LogPctSAT)
cor(cor_data, use = "complete.obs") %>%
kable(digits = 3, caption = "Korelasyon Matrisi")| Expend | PTratio | Salary | PctSAT | SATcombined | LogPctSAT | |
|---|---|---|---|---|---|---|
| Expend | 1.000 | -0.371 | 0.870 | 0.593 | -0.381 | 0.561 |
| PTratio | -0.371 | 1.000 | -0.001 | -0.213 | 0.081 | -0.132 |
| Salary | 0.870 | -0.001 | 1.000 | 0.617 | -0.440 | 0.613 |
| PctSAT | 0.593 | -0.213 | 0.617 | 1.000 | -0.887 | 0.961 |
| SATcombined | -0.381 | 0.081 | -0.440 | -0.887 | 1.000 | -0.926 |
| LogPctSAT | 0.561 | -0.132 | 0.613 | 0.961 | -0.926 | 1.000 |
library(GGally)
ggpairs(datahw, columns = c("Expend", "PTratio", "Salary", "PctSAT", "SATcombined", "LogPctSAT"))
# Expend ve SATcombined arasında negatif korelasyon var -0.381 (p < .01).
# Salary ve Expend: Korelasyon 0.870. Çok yüksek ve pozitif. Yani öğretmen maaşları arttıkça eğitim harcamaları da artıyor.
# PTratio (Öğrenci/Öğretmen Oranı): SAT puanı ile hiç ilişki yok.# Model 1: (Simple Regression)
model1 <- lm(SATcombined ~ Expend, data = datahw)
# Model 2: (Multiple Regression)
model2 <- lm(SATcombined ~ Expend + LogPctSAT, data = datahw)
library(broom)
sqrt(glance(model1)[,1]) #r.squared değerinin karekoku alınır sonuç=0.380sqrt(glance(model2)[,1]) #r.squared değerinin karekoku alınır sonuç= 0.941# Model 1 (R = .381): harcama tek başınayken ilişki zayıf ve (korelasyon olarak bakıldığında) negatiftir.
# Model 2 (R = .941): LogPctSAT eklendiğinde Korelasyon Katsayısı .941 e yükseldi.
library(stargazer)
stargazer(model1, model2, type = "text")##
## =================================================================
## Dependent variable:
## ---------------------------------------------
## SATcombined
## (1) (2)
## -----------------------------------------------------------------
## Expend -20.892*** 11.129***
## (7.328) (3.264)
##
## LogPctSAT -78.203***
## (4.471)
##
## Constant 1,089.294*** 1,147.100***
## (44.390) (16.701)
##
## -----------------------------------------------------------------
## Observations 50 50
## R2 0.145 0.886
## Adjusted R2 0.127 0.881
## Residual Std. Error 69.909 (df = 48) 25.782 (df = 47)
## F Statistic 8.128*** (df = 1; 48) 182.834*** (df = 2; 47)
## =================================================================
## Note: *p<0.1; **p<0.05; ***p<0.01Sonuçlar: Model (1): Sadece harcamaların (Expend) SAT puanları üzerindeki etkisine baktığımızda sonuç: (Estimate): -20.892. Yani eyalet bazında eğitime yapılan harcama arttıkça, SAT puanları yaklaşık 21 puan düşüyor gibi görünüyor.Anlamlılık: *** (p < 0.01) düzeyinde istatistiksel olarak anlamlı.Açıklayıcılık (R2): 0.145. Harcamalar, SAT puanlarındaki değişimin sadece %14.5’ini açıklıyor.
Model (2): Çoklu Regresyon: Sınava katılım oranını (LogPctSAT) dahil ettiğimizde: Katsayı Değişimi (Expend): Harcamanın katsayısı -20.892’den +11.129’a fırladı ve işareti pozitife döndü. Sınava katılım oranları eşit olan (kontrol edilen) eyaletleri kıyasladığımızda, eğitime harcanan her bir birim para, puanları aslında 11 puan artırıyor. Katılımın Etkisi (LogPctSAT): Katsayısı -78.203. Yani sınava katılım yüzdesi arttıkça eyaletin ortalama puanı düşüyor.(R2): 0.886. Bu iki değişken birlikte, eyaletler arası başarı farkının %88.6’sını açıklıyor.
ANOVA tablosundaki F değeri ( 182.834), p < .001 olduğu için regresyon modeli veriye iyi uyum sağlamaktadır.
# Kestirimin standart hatası ikinci modelde düştü.
glance(model1)[,3]glance(model2)[,3]tidy(model2) %>% kable(digit=3)| term | estimate | std.error | statistic | p.value |
|---|---|---|---|---|
| (Intercept) | 1147.100 | 16.701 | 68.684 | 0.000 |
| Expend | 11.129 | 3.264 | 3.409 | 0.001 |
| LogPctSAT | -78.203 | 4.471 | -17.490 | 0.000 |
artıksat1 <- lm(SATcombined ~ LogPctSAT, data = datahw)$residuals # Sat'ın LogPctSAT'tan yordanmayan kısmı
artıkexp1<- lm(Expend ~ LogPctSAT, data = datahw)$residuals # Expend'in LogPctSAT'tan yordanmayan kısmı
lm(artıksat1 ~artıkexp1,data=data.frame(artıksat1,artıkexp1))$coefficients %>% kable(digit=3)| x | |
|---|---|
| (Intercept) | 0.000 |
| artıkexp1 | 11.129 |
# Elde ettiğim katsayı 11.130, daha önce kurduğumuz çoklu regresyon modelindeki Expend değişkenine ait olan katsayısın aynısı. Bu değer, sınava katılım oranı kontrol altına alındıktan sonra, harcama düzeylerindeki bir birimlik artışın SAT puanını 11.129 artırmaya eğilimli olduğunu önermektedir. # Model 3: Üç yordayıcı ile regresyon
model3 <- lm(SATcombined ~ Expend + LogPctSAT + PTratio, data = datahw)
# Sonuçlar ve Çoklu Doğrusallık (VIF) değerleri için
summary(model3)##
## Call:
## lm(formula = SATcombined ~ Expend + LogPctSAT + PTratio, data = datahw)
##
## Residuals:
## Min 1Q Median 3Q Max
## -60.153 -14.430 -2.416 18.020 56.068
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1131.9823 39.7887 28.450 < 2e-16 ***
## Expend 11.6651 3.5329 3.302 0.00186 **
## LogPctSAT -78.3910 4.5332 -17.293 < 2e-16 ***
## PTratio 0.7442 1.7743 0.419 0.67687
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 26.01 on 46 degrees of freedom
## Multiple R-squared: 0.8865, Adjusted R-squared: 0.8791
## F-statistic: 119.8 on 3 and 46 DF, p-value: < 2.2e-16library(car)
# VIF değerlerini al
vif_values <- vif(model3)
# Tolerans değerlerini hesapla (1 / VIF)
tolerance_values <- 1 / vif_values
# Sonuçları tablo gibi gör
data.frame(VIF = vif_values, Tolerance = tolerance_values)# Sonuçlar
# PTratio'nun SAT puanlarını tahmin etmeye hiçbir anlamlı katkısı yok.
# VIF ve Tolerans değerleri, çoklu doğrusallık gösterir; analizindeki düşük VIF (10'un altı) ve yüksek Tolerans (0.10'un üstü) değerleri, modelinin güvenilir ve değişkenlerinin birbirinden yeterince bağımsız olduğunu gösterir.# Standardized Regression Coefficients (Beta)
# B (Unstandardized): "1 dolarlık" artışın etkisini gösterir.
# Beta (Standardized): "1 standart sapmalık" artışın etkisini gösterir.
library(QuantPsyc)
lm.beta(model2) %>% kable(digit=3)| x | |
|---|---|
| Expend | 0.203 |
| LogPctSAT | -1.040 |
# Expend (Beta = 0.203): Harcamadaki 1 standart sapmalık artış, SAT puanını 0.203 standart sapma artırır.
# LogPctSAT (Beta = -1.040): Katılım oranındaki 1 standart sapmalık artış, SAT puanını 1.040 standart sapma düşürür.
# Eyaletler arası SAT puanı farklarını açıklarken, sınava katılım oranı (LogPctSAT), yapılan harcamadan (Expend) daha önemli bir rol oynamaktadır.Çoklu Regresyonda Yeniden Örnekleme (Bootstrapping) Yaklaşımı Geleneksel t-testleri verinin belirli bir teorik dağılıma (normal dağılım gibi) uyduğunu varsayarken; Yeniden Örnekleme (Resampling/Bootstrapping) yöntemi, mevcut veriyi bir “evren” kabul ederek varsayımlardan bağımsız, daha dayanıklı sonuçlar sunar. Bu yöntemde, orijinal veri setinden binlerce kez “yerine koyarak” (with replacement) yeni örneklemler çekilir ve her seferinde regresyon katsayıları (b_j) yeniden hesaplanır. Sonuçlar %95’lik Güven Aralığı ile (Percentile Confidence Interval) incelenir: Eğer hesaplanan aralık 0 (sıfır) değerini içermiyorsa, o değişkenin bağımlı değişken üzerindeki etkisi istatistiksel olarak anlamlı kabul edilir.
anova(model2)# Residual Variance : anova(model2) komutunu kullandığında "Residuals" satırındaki "Mean Sq" değeri 665 değeridir.Modelin açıklayamadığı "hata payının" karesidir. Bu değer ne kadar küçükse, yordayıcılar bağımlı değişkeni o kadar hassas açıklıyor demektir.
summary(model2)##
## Call:
## lm(formula = SATcombined ~ Expend + LogPctSAT, data = datahw)
##
## Residuals:
## Min 1Q Median 3Q Max
## -61.515 -13.616 -2.572 16.541 54.901
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1147.100 16.701 68.684 < 2e-16 ***
## Expend 11.129 3.264 3.409 0.00135 **
## LogPctSAT -78.203 4.471 -17.490 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 25.78 on 47 degrees of freedom
## Multiple R-squared: 0.8861, Adjusted R-squared: 0.8813
## F-statistic: 182.8 on 2 and 47 DF, p-value: < 2.2e-16# Standard Error of Estimate: summary(model2) çıktısındaki Residual standard error değeridir.Tahminlerimizin gerçek değerlerden ortalamada ne kadar saptığını (puan bazında) gösterir. Yani bu modelle bir eyaletin SAT puanını tahmin ederken yaklaşık 25.78 puanlık bir yanılma payımız var demektir.𝑅2 değeri çoklu korelasyon katsayısı (multiple correlation coefficient) olup bağımlı değişkenin gözlenen değerleri ile bağımsız değişkenlerin en iyi doğrusal kombinasyonu arasındaki korelasyondur. adj𝑅2 değeri, 𝑅2 değerinin modeldeki bağımsız değişken sayısı için modifiye edilmiş versiyonudur. adj𝑅2 değeri, yeni eklenen bağımsız değişken modeli şans eseri beklenenden daha fazla geliştirirse artar, daha az geliştirirse azalır.
Örneklem Büyüklüğü Kuralları: Her yordayıcı için en az 10 gözlem Harris’in Kuralı: N, p’den en az 50 fazla olmalı (N > p + 50). Güç Analizi (Cohen): Belirli bir istatistiksel güç (.80) için çok daha fazla (yüzlerce) katılımcı gerekebilir. Darlington’un Kuralı: Ne kadar çok, o kadar iyi (More is better).
Suppressor Variables: Normal şartlarda, bir değişken bağımlı değişkenle pozitif korelasyona sahipse, regresyon katsayısının da pozitif olmasını bekleriz. Ancak bazen bir değişken, modele girdiğinde beklenmedik bir işaret alabilir. Darlington’ın “Tarih Sınavı” örneği: Bir tarih sınavı yapıyorsunuz ama süre çok kısıtlı. Sınav puanı hem tarih bilgisini hem de okuma hızını ölçüyor. Okuma hızı aslında tarih bilgisinden bağımsızdır ama eğer modele “Okuma Hızı” değişkenini eklerseniz, model hızlı okuduğu için yüksek puan alanları “cezalandırır”, yavaş okuduğu için düşük alanları “ödüllendirir”. Böylece okuma hızının yarattığı hatayı bastırır ve saf tarih bilgisini ortaya çıkarır.
Regression Diagnostics Cook’s D: Bağımlı değişkenlerdeki potansiyel uç değerlerin belirlenmesinde kullanışlı bir istatistiktir. Uzaklık için en yaygın ölçüm artıktır.
DFBETA: Cook’s D etkinin genel bir ölçümü olarak düşünülebilir. Gözlemin eklenmesiyle her katsayının nasıl değiştiğini ölçen daha spesifik bir ölçüm ele alınabilir.
Leverage (hi): Bağımsız değişkenlerdeki potansiyel uç değerlerin belirlenmesinde kullanışlı bir istatistiktir. Levarage bir gözlemin bir bağımsız değişkene, 𝑋𝑗, göre olağan dışı olma derecesini ölçer.
Influence (Etki): Etkili bir gözlem uzaklık ve/veya leverage için yüksek değere sahip olan ve modelin kesişim ve eğim katsayılarını anlamlı olarak etkileyen bir gözlemdir.
Varsayımların Testi: Homoscedasticity (Eş Varyanslılık)
library(olsrr)
ols_plot_cooksd_bar(model)
# Cook’s D için kesme noktası 4/50# DFBETA için kesme noktası ise 2/√50
# hat değerleri ise levarge a karşılık geliyor.
# DFBETA
ols_plot_dfbetas(model2)
influence.measures(model2, infl = influence(model2))## Influence measures of
## lm(formula = SATcombined ~ Expend + LogPctSAT, data = datahw) :
##
## dfb.1_ dfb.Expn dfb.LPSA dffit cov.r cook.d hat inf
## 1 -0.03619 0.01562 0.01477 -0.0416 1.122 5.88e-04 0.0512
## 2 0.13365 -0.16679 0.05058 -0.1896 1.210 1.22e-02 0.1321 *
## 3 0.00734 -0.00931 0.00620 0.0126 1.116 5.41e-05 0.0448
## 4 -0.44713 0.11127 0.29632 -0.5457 0.854 9.22e-02 0.0610
## 5 -0.01047 0.02146 -0.02131 -0.0296 1.135 2.98e-04 0.0607
## 6 0.09115 -0.11334 0.09938 0.2365 0.965 1.83e-02 0.0271
## 7 -0.07270 0.06883 0.00214 0.0933 1.198 2.96e-03 0.1132 *
## 8 -0.00691 0.00263 0.00705 0.0143 1.116 6.95e-05 0.0448
## 9 -0.01719 0.07115 -0.10472 -0.1512 1.071 7.69e-03 0.0392
## 10 -0.05894 0.19442 -0.23333 -0.2855 1.085 2.72e-02 0.0771
## 11 0.00552 0.02463 -0.05416 -0.0769 1.101 2.01e-03 0.0403
## 12 -0.02478 0.02204 -0.00548 -0.0298 1.126 3.02e-04 0.0534
## 13 0.04213 0.11674 -0.16013 0.2507 0.988 2.06e-02 0.0348
## 14 -0.00337 0.04908 -0.08136 -0.1080 1.101 3.95e-03 0.0464
## 15 0.08419 0.06793 -0.16682 0.1964 1.125 1.30e-02 0.0788
## 16 0.05776 0.06648 -0.12959 0.1730 1.074 1.01e-02 0.0458
## 17 -0.08631 0.00988 0.06149 -0.1330 1.064 5.95e-03 0.0321
## 18 -0.04732 0.01501 0.02485 -0.0595 1.107 1.20e-03 0.0417
## 19 -0.01502 -0.01103 0.04367 0.0627 1.110 1.34e-03 0.0447
## 20 -0.03118 0.01893 0.02062 0.0562 1.113 1.07e-03 0.0459
## 21 -0.11435 0.04904 0.10042 0.2042 1.076 1.40e-02 0.0540
## 22 0.00674 -0.03876 0.03827 -0.0507 1.153 8.75e-04 0.0767
## 23 0.09970 0.18875 -0.31029 0.4009 0.926 5.13e-02 0.0501
## 24 -0.52567 0.12428 0.37505 -0.6271 0.907 1.23e-01 0.0881
## 25 0.04270 0.01070 -0.05079 0.0782 1.100 2.08e-03 0.0397
## 26 0.07803 -0.02365 -0.00633 0.2105 0.952 1.44e-02 0.0205
## 27 0.02143 0.03398 -0.05897 0.0770 1.112 2.01e-03 0.0484
## 28 -0.07490 0.10202 -0.08226 -0.1620 1.055 8.80e-03 0.0352
## 29 -0.00608 -0.26400 0.45383 0.5633 0.818 9.70e-02 0.0571
## 30 0.41952 -0.47386 0.10563 -0.5408 1.219 9.70e-02 0.1918 *
## 31 0.02998 -0.01638 -0.00667 0.0354 1.110 4.27e-04 0.0406
## 32 0.39707 -0.43394 0.07866 -0.5070 1.198 8.53e-02 0.1761 *
## 33 -0.05399 0.15105 -0.17055 -0.2132 1.115 1.53e-02 0.0765
## 34 0.24827 0.01023 -0.25769 0.3608 1.025 4.27e-02 0.0695
## 35 0.00168 0.00610 -0.00389 0.0256 1.088 2.24e-04 0.0212
## 36 -0.01370 0.00361 0.00806 -0.0178 1.111 1.08e-04 0.0407
## 37 -0.06004 -0.01183 0.14015 0.2611 0.968 2.22e-02 0.0324
## 38 0.06236 -0.02620 -0.05867 -0.1222 1.097 5.05e-03 0.0467
## 39 0.05741 -0.03736 -0.03176 -0.0936 1.115 2.97e-03 0.0530
## 40 -0.16086 0.38493 -0.39698 -0.5029 0.993 8.16e-02 0.0888
## 41 -0.04798 -0.00198 0.04980 -0.0697 1.141 1.65e-03 0.0695
## 42 0.29227 -0.20063 -0.01398 0.3372 0.960 3.68e-02 0.0454
## 43 -0.02983 0.06962 -0.07630 -0.1047 1.113 3.72e-03 0.0536
## 44 -0.04123 0.01626 0.02130 -0.0451 1.180 6.92e-04 0.0973
## 45 -0.02768 0.00126 0.04392 0.0735 1.106 1.84e-03 0.0433
## 46 0.03045 -0.11450 0.14462 0.1772 1.118 1.06e-02 0.0709
## 47 0.00684 -0.07508 0.13400 0.2035 1.028 1.38e-02 0.0353
## 48 -0.06374 -0.14771 0.18258 -0.4127 0.735 5.08e-02 0.0254 *
## 49 -0.02134 0.19492 -0.20902 0.2603 1.119 2.27e-02 0.0888
## 50 -0.05439 -0.16401 0.23637 -0.3127 0.991 3.20e-02 0.0483# leverage
library(olsrr)
ols_plot_resid_lev(model2)
# influence
ols_plot_dffits(model2)
# Modelin tanı grafikleri(4 ana grafik)
par(mfrow = c(2, 2))
plot(model2)
MODEL KARŞILAŞTIRMA
model_full <- lm(SATcombined ~ Expend + LogPctSAT + PTratio + Salary, data = datahw)
model_red <- lm(SATcombined ~ Expend + LogPctSAT, data = datahw)
# İki modeli kıyasla (F-testi) # MODEL 1 daha iyi
anova(model_red, model_full)# Yuvalanmamış Modeller İçin (AIC):
model_a <- lm(SATcombined ~ Expend + PctSAT, data = datahw)
model_b <- lm(SATcombined ~ Expend + LogPctSAT, data = datahw)
# AIC değerlerini gör (Küçük olan DAHA İYİ, model b)
AIC(model_a, model_b)Regresyon Denklemi Oluşturma Değişken Seçim Yöntemleri: Tüm Alt Küme Regresyonu (All Subsets), Geriye Doğru Eleme (Backward Elimination), Aşamalı Regresyon (Stepwise).
Importance of Variables: Beta (b_j) Yanılgısı: Katsayılar örneklemden örnekleme çok değişebilir.
Kullanışlılık (Usefulness): Bir değişkenin önemini ölçmenin en iyi yolu, o değişken çıkarıldığında R2’de meydana gelen düşüştür.
Kayıp Verilerle Başa Çıkma: Listwise Deletion, Pairwise Deletion, Imputation (Atama).
Çapraz Geçerlilik (Cross-Validation): Geliştirdiğimiz bir regresyon denklemi, kendi verimizde çok iyi çalışabilir (R2 yüksek çıkar) ama yeni bir veriye uygulandığında performansı düşebilir.
Mediating and Moderating Relationships
Aracılık analizi
X (Bağımsız): Anne şefkati (Maternal Care)
M (Aracı): Özsaygı (Self-esteem)
Y (Bağımlı): Anne özyeterliliği (Self-efficacy)
library(haven)
MediationData <- read_sav("howell/MediationData.sav")
View(MediationData)# Adım 1: Toplam Etki (X -> Y)
# Bağımsız değişkenin bağımlı değişken üzerinde anlamlı bir etkisi olmalı (c yolu)
modelm1 <- lm(Efficacy ~ MatCare, data = MediationData)
# Adım 2: X'in Aracı Değişken Üzerindeki Etkisi (X -> M)
# Bağımsız değişken aracı değişkeni anlamlı etkilemeli (a yolu)
modelm2 <- lm(Esteem ~ MatCare, data = MediationData)
#Adım 3. YOL B (M on Y, X kontrol edildiğinde )
modelm3 <- lm(Efficacy ~ Esteem + MatCare, data=MediationData)
# Adım 4. YOL C - ters yol (Y on X, M kontrol edildiğinde)
modelm4 <- lm(MatCare ~ Efficacy + Esteem, data = MediationData)
# Özet tablo
stargazer(modelm1, modelm2, modelm3, modelm4, type = "text", title = "Baron and Kenny Yöntemi ile Aracılık Etkisi",digits = 2,
font.size ="tiny")##
## Baron and Kenny Yöntemi ile Aracılık Etkisi
## =======================================================================================================
## Dependent variable:
## -----------------------------------------------------------------------------------
## Efficacy Esteem Efficacy MatCare
## (1) (2) (3) (4)
## -------------------------------------------------------------------------------------------------------
## Efficacy 0.33
## (0.25)
##
## Esteem 0.15*** 0.39***
## (0.05) (0.11)
##
## MatCare 0.11** 0.36*** 0.06
## (0.04) (0.09) (0.04)
##
## Constant 3.27*** 2.26*** 2.94*** 0.78
## (0.14) (0.29) (0.17) (0.85)
##
## -------------------------------------------------------------------------------------------------------
## Observations 92 92 92 92
## R2 0.07 0.16 0.16 0.18
## Adjusted R2 0.06 0.15 0.14 0.16
## Residual Std. Error 0.24 (df = 90) 0.50 (df = 90) 0.23 (df = 89) 0.55 (df = 89)
## F Statistic 6.92** (df = 1; 90) 17.45*** (df = 1; 90) 8.34*** (df = 2; 89) 9.63*** (df = 2; 89)
## =======================================================================================================
## Note: *p<0.1; **p<0.05; ***p<0.01- Anne şefkati, özyeterliliği tek başına anlamlı ve pozitif yönde etkiliyor. İlk şart sağlandı.
- Anne şefkati, özsaygıyı anlamlı bir şekilde yorduyor. İkinci şart sağlandı.
- Özsaygı, anne şefkati kontrol edildiğinde özyeterlilik üzerinde anlamlı bir etkiye sahip. Üçüncü şart sağlandı.
- Aracı değişken modele girdiğinde, anne şefkatinin özyeterlilik üzerindeki doğrudan etkisi istatistiksel olarak anlamlı değil.
#Sobel Test
library(bda)
attach(MediationData)
mediation.test(Esteem, MatCare, Efficacy)# Anne şefkati ile özyeterlilik arasındaki dolaylı etkinin anlamlılığı Sobel testi ile doğrulanmıştır (z = 2.448, p = 0.014). Özsaygı değişkeninin aracı rolünün istatistiksel olarak güvenilir olduğunu kanıtlamaktadır. #Bootstrap öncesi varsayım testi
library(mediation)
library(gvlma)
fit1 <- lm(Esteem ~ MatCare, data=MediationData)
fit2 <- lm(Efficacy ~ MatCare + Esteem, data=MediationData)
gvlma(fit1)##
## Call:
## lm(formula = Esteem ~ MatCare, data = MediationData)
##
## Coefficients:
## (Intercept) MatCare
## 2.2574 0.3637
##
##
## ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
## USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
## Level of Significance = 0.05
##
## Call:
## gvlma(x = fit1)
##
## Value p-value Decision
## Global Stat 16.7804 0.002132 Assumptions NOT satisfied!
## Skewness 10.2318 0.001380 Assumptions NOT satisfied!
## Kurtosis 3.9237 0.047611 Assumptions NOT satisfied!
## Link Function 2.5246 0.112082 Assumptions acceptable.
## Heteroscedasticity 0.1004 0.751371 Assumptions acceptable.gvlma(fit2)##
## Call:
## lm(formula = Efficacy ~ MatCare + Esteem, data = MediationData)
##
## Coefficients:
## (Intercept) MatCare Esteem
## 2.93609 0.05653 0.14597
##
##
## ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
## USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
## Level of Significance = 0.05
##
## Call:
## gvlma(x = fit2)
##
## Value p-value Decision
## Global Stat 13.17229 0.010464 Assumptions NOT satisfied!
## Skewness 7.55963 0.005969 Assumptions NOT satisfied!
## Kurtosis 0.04858 0.825551 Assumptions acceptable.
## Link Function 5.32459 0.021027 Assumptions NOT satisfied!
## Heteroscedasticity 0.23949 0.624575 Assumptions acceptable.plot(gvlma(fit1),1)
#Bootstrap
med_boot<- mediate(fit1, fit2, boot=TRUE, sims=999, treat="MatCare", mediator="Esteem")
summary(fitMedBoot)##
## Causal Mediation Analysis
##
## Nonparametric Bootstrap Confidence Intervals with the Percentile Method
##
## Estimate 95% CI Lower 95% CI Upper p-value
## ACME 0.280784 0.146694 0.421051 < 2e-16 ***
## ADE -0.111790 -0.281894 0.104693 0.28629
## Total Effect 0.168993 -0.002348 0.339581 0.05205 .
## Prop. Mediated 1.661509 -1.356902 10.179730 0.05205 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Sample Size Used: 100
##
##
## Simulations: 999ACME: Güven aralığı sıfırı (0) içermiyor Ayrıca p-değeri (< 2e-16) anlamlı. Bu, aracılık etkisinin istatistiksel olarak var olduğunu kanıtlar.
ADE: Bu, aracı değişken (M) kontrol edildikten sonra X’in Y üzerindeki doğrudan etkisini gösterir. İstatistiksel olarak anlamlı değil.
Total Effect: Toplam Etki: Toplam etki istatistiksel olarak anlamlı değil.
Prop. Mediated: Aracılık Oranı Toplam etkinin ne kadarının aracı değişken üzerinden geçtiğini yüzde olarak belirtir.İstatistiksel olarak anlamlı değil.
Sonuç: Yapılan 999 simülasyonlu Bootstrap analizi sonucunda, dolaylı etkinin (ACME) istatistiksel olarak anlamlı olduğu saptanmıştır (b = 0.28, p < 0.001, %95 CI [0.14, 0.43]).
Moderasyon
Zorluklar (Hassles) arttıkça psikolojik belirtiler de artar. Sosyal desteği yüksek olan kişilerde, zorluklar artsa bile belirtiler çok az artar. Ancak sosyal desteği düşük olanlarda, zorluklar arttıkça belirtiler çok daha hızlı ve şiddetli artar.
Değişken çarpımı öncesinde merekezleme işlemi yapıyoruz çünkü çoklu Doğrusallığı önlemek (Multicollinearity) istiyoruz ve yorumlama kolaylığı sağlar.
library(haven)
Moderating <- read_sav("howell/Moderating.sav")
View(Moderating)# Değişkenler arası ilişkiler(merkezlenmiş ve merkezlenmemiş)
# Kitaptaki değişkenler
secilen_degiskenler <- Moderating[, c("junehass", "junesupp", "junesymp",
"hassupp", "chassles", "csupport", "chassupp")]
# Korelasyon matrisini hesapla
kor_matrisi <- cor(secilen_degiskenler, use = "complete.obs")
# Tabloyu daha okunabilir kılmak için yuvarla
round(kor_matrisi, 3)## junehass junesupp junesymp hassupp chassles csupport chassupp
## junehass 1.000 0.130 0.452 0.784 1.000 0.130 0.401
## junesupp 0.130 1.000 0.077 -0.450 0.130 1.000 0.168
## junesymp 0.452 0.077 1.000 0.403 0.452 0.077 -0.010
## hassupp 0.784 -0.450 0.403 1.000 0.784 -0.450 -0.002
## chassles 1.000 0.130 0.452 0.784 1.000 0.130 0.401
## csupport 0.130 1.000 0.077 -0.450 0.130 1.000 0.168
## chassupp 0.401 0.168 -0.010 -0.002 0.401 0.168 1.000model_moderasyon <- lm(junesymp ~ chassles + csupport + chassupp, data = Moderating)
summary(model_moderasyon)##
## Call:
## lm(formula = junesymp ~ chassles + csupport + chassupp, data = Moderating)
##
## Residuals:
## Min 1Q Median 3Q Max
## -38.734 -13.379 -0.444 8.689 35.521
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 89.584951 2.291529 39.094 < 2e-16 ***
## chassles 0.085952 0.019211 4.474 4.21e-05 ***
## csupport 0.146381 0.305239 0.480 0.6336
## chassupp -0.005066 0.002363 -2.144 0.0367 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 16.89 on 52 degrees of freedom
## (10 observations deleted due to missingness)
## Multiple R-squared: 0.3885, Adjusted R-squared: 0.3532
## F-statistic: 11.01 on 3 and 52 DF, p-value: 1.045e-05Sonuç: R2: 0.3885. Model (Zorluklar, Destek ve Etkileşim) bağımlı değişken olan psikolojik belirtilerin (Symptoms) %38.8’ini açıklıyor. Zorluklar arttıkça psikolojik belirtiler anlamlı şekilde artıyor. Sosyal desteğin belirtiler üzerinde doğrudan, tek başına bir etkisi yok. Etkileşim katsayısı istatistiksel olarak anlamlı olması, sosyal destek arttıkça, zorlukların belirtiler üzerindeki o yıkıcı etkisinin azaldığını gösterir.
#Plotting
# Kitaptakinden faklı çıktı görsel anlamadım, sanırım biz bu data üzerinde geçen sene çalışmıştık diye hatırlıyorum ama hatayı anlamadım açıkcası.
library(rockchalk)
ps <- plotSlopes(model_moderasyon,
plotx = "chassles",
modx = "csupport",
modxVals = "std.dev", # Düşük, orta, yüksek destek için +/- 1 SD
xlab = "Stres (Hassles)",
ylab = "Belirtiler (Symptoms)",
main = "Sosyal Desteğin Moderasyon Etkisi",
xlim = c(-200, 200))
CHAPTER 3- Modern Psychometrics With R
Çok değişkenli regresyon (multivariate regression) yol modeli kurma yöntemleri:
Klasik doğrusal model bağlamında kalmak: Her bir yanıt değişkeni için ayrı bir çoklu regresyon modeli tahmin etmek ve ardından MANOVA’dan (çok değişkenli varyans analizi) bilinen çeşitli ölçümleri hesaplayarak çıkarımda bulunmak.
Yol analizi (path-analytic) çerçevesi: Modeli bu yapı üzerinden kurmak.
Karma etkili modeller (mixed-effects models): Bağımlı değişkenler için rastgele bir etki (random effect) belirlemek.
Uygulama Örneği
library("MPsychoR")
data("Bergh")#multivariate regression
library("lavaan")
mvreg.model <- '
EP ~ b11*A1 + b12*A2 + b13*O1 + b14*O2
DP ~ b21*A1 + b22*A2 + b23*O1 + b24*O2'
fitmvreg2 <- sem(mvreg.model, data = Bergh)
# EP ve DP: Etnik önyargı ve engelli bireylere yönelik önyargıyı temsil eden iki ayrı bağımlı değişkenin var.
# A1 ve A2 (Agreeableness / Uyumluluk) ile ilgili iki madde
# O1 ve O2 (Openness to Experience / Deneyime Açıklık) ile ilgili iki madde
# b11,A1 (Uyumluluk 1),EP (Etnik Önyargı),A1'in Etnik Önyargı üzerindeki etkisi.
# b12,A2 (Uyumluluk 2),EP (Etnik Önyargı),A2'in Etnik Önyargı üzerindeki etkisi.
# b13,O1 (Açıklık 1),EP (Etnik Önyargı),O1'in Etnik Önyargı üzerindeki etkisi.
# b14,O2 (Açıklık 2),EP (Etnik Önyargı),O2'nin Etnik Önyargı üzerindeki etkisi.
# b21,A1 (Uyumluluk 1),DP (Engelli Önyargısı),A1'in Engelli Önyargısı üzerindeki etkisi.
# b22,A2 (Uyumluluk 2),DP (Engelli Önyargısı),A2'nin Engelli Önyargısı üzerindeki etkisi.
# b23,O1 (Açıklık 1),DP (Engelli Önyargısı),O1'in Engelli Önyargısı üzerindeki etkisi.
# b24,O2 (Açıklık 2),DP (Engelli Önyargısı),O2'nin Engelli Önyargısı üzerindeki etkisi.parameterEstimates(fitmvreg2)[c(1:8, 11),]# Sonuç: A1 (ilk uyumluluk maddesi) değişkeninin ne EP ne de DP üzerinde anlamlı bir etkisi yoktur.A1 dışındaki tüm yordayıcılar (A2, O1, O2) her iki önyargı türü üzerinde de anlamlı bir etkiye sahip.library("semPlot")
semPaths(fitmvreg2, what = "est", edge.label.cex = 1, layout = "tree", residuals = FALSE, edge.color = 1,
esize = 1, rotation = 3, sizeMan = 8, asize = 2.5, fade = FALSE, optimizeLatRes = TRUE)
parameterEstimates(fitmvreg2)[c(1:8, 11),]Moderator and Mediator Models Düzenleyici (moderator); bağımsız bir değişken ile bağımlı bir değişken arasındaki ilişkinin yönünü ve/veya gücünü etkileyen nitel (örneğin cinsiyet, ırk, sınıf) veya nicel (örneğin ödül seviyesi) bir değişkendir.
Aracı (mediator) değişkenin ise; bağımsız değişken ile bağımlı değişken arasındaki ilişkiyi açıklayan/karşılayan bir işlev gördüğü söylenebilir.
Düzenleyici, doğrudan X ile etkileşime girer (Örneğin: X = alkol tüketimi, Z = sosyal bağlam, Y = sosyal kabul). Aracı, X ile Y arasında “köprü” görevi görür (Örneğin: X = alkol tüketimi, M = eve araçla dönmek, Y = hapishanede uyanmak).
Düzenleme (moderation) bir etkileşim (interaction) olarak ifade edilir: Bir regresyon modelinde Z değişkeninin düzenleyici etkisini hesaba katmak için, X ve Z yordayıcıları arasında bir etkileşime izin vermemiz gerekir.
data("Paskvan")
wintense.c <- scale(Paskvan$wintense, scale = FALSE) ## merkezleme
fit.YX <- lm(cogapp ~ wintense.c, data = Paskvan) ## Y on X: "Bilişsel Değerlendirme" (Y), İş Yoğunluğu (X) tarafından yordanıyor mu?
round(summary(fit.YX)$coefficients, 4)## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.7265 0.0227 164.4470 0
## wintense.c 0.5458 0.0237 22.9954 0pclimate.c <- scale(Paskvan$pclimate, scale = FALSE) ## merkezleme
fit.YZ <- lm(cogapp ~ pclimate.c, data = Paskvan) ## Y on Z: Bilişsel Değerlendirme" (Y), İş yerindeki katılımcı iklim (X) tarafından yordanıyor mu?
round(summary(fit.YZ)$coefficients, 4)## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.7265 0.0272 136.8341 0
## pclimate.c -0.3324 0.0304 -10.9408 0# iş yoğunluğu arttıkça, bilişsel değerlendirme artıyor.
# Katılımcı iklim arttıkça, bilişsel değerlendirme azalıyor.library("QuantPsyc")
fit.mod <- moderate.lm(x = wintense, z = pclimate, y = cogapp, data = Paskvan)
round(summary(fit.mod)$coefficients, 4)## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.7097 0.0227 163.0954 0.0000
## mcx 0.5003 0.0241 20.7653 0.0000
## mcz -0.1790 0.0257 -6.9765 0.0000
## mcx:mcz -0.0663 0.0236 -2.8094 0.0051fit.ss <- sim.slopes(fit.mod, Paskvan$pclimate) # Basit Eğim Analizi (Simple Slopes)
round(fit.ss, 4)# Düzenleyici değişkeni (Z) üç seviyeye ayırır: Düşük (Ortalamanın 1 standart sapma altı), Orta (Ortalama) ve Yüksek (Ortalamanın 1 standart sapma üstü).Her bir seviye için X'in Y üzerindeki etkisini (eğimini) ayrı ayrı hesaplar.
# moderate.lm: Bu fonksiyon, X (İş yoğunluğu) ve Z (Katılımcı iklim) değişkenlerini otomatik olarak merkezler, bir etkileşim terimi (X*Z) oluşturur ve regresyon modelini kurar.Yorum: Merkezlenmiş iş yoğunluğu (mcx) ile merkezlenmiş katılımcı iklim (mcz) arasındaki etkileşim anlamlıdır.
at zLow (Düşük İklim): Katsayı 0.5598. Katılımcı iklim zayıfken, iş yoğunluğundaki artış bilişsel değerlendirmeyi en şiddetli şekilde artırıyor.
at zMean (Orta İklim): Katsayı 0.5003. Etki biraz daha hafifliyor.
at zHigh (Yüksek İklim): Katsayı 0.4407. Katılımcı iklim en yüksek seviyedeyken, iş yoğunluğunun bilişsel değerlendirme üzerindeki olumsuz etkisi en düşük seviyesine iniyor.
Mediator Model
Aracılık analizi için önce regresyon analizi yapmak gerekir:
fit.MX: Bağımsız değişkenin (X: iş yoğunluğu), aracı değişken (M: bilişsel değerlendirme) üzerindeki etkisini ölçer.
fit.YXM: Hem X’in hem de M’in, bağımlı değişken (Y: duygusal tükenmişlik) üzerindeki etkisini aynı anda ölçer.
mediate() Fonksiyonu ile Analiz: Bu fonksiyon, yukarıdaki iki modeli birleştirerek aracılık etkisini hesaplar. Analizi 999 kez tekrarlar (bootstrapping)
ACME Dolaylı Etki (Aracılık Etkisi): X →M → Y yolu. ADE Doğrudan Etki: Aracı değişken devreden çıkarıldığında X’in Y üzerindeki kalan etkisi. Total Effect: Toplam Etki: ACME + ADE toplamı. Prop. Mediated: Aracılık Oranı: Toplam etkinin ne kadarının aracı değişken üzerinden geçtiğini gösterir.
Sonuçlar: Bu çıktıda ACME (ortalama nedensel aracılık etkisi), 95% bootstrap güven aralığı (CI) ile birlikte ab dolaylı etkisini temsil eder. ACME = 0 (sıfır hipotezi) güven aralığının içinde yer almadığı için, bilişsel değerlendirmenin (M), iş yoğunluğu (X) ile duygusal tükenmişlik (Y) arasındaki ilişkiye aracılık ettiği sonucuna varırız.
ADE, modelimizdeki c’ katsayısı olan ortalama doğrudan etkiyi ifade eder. Toplam etki ise basitçe ab + c’ toplamıdır. Son satır, dolaylı etkinin toplam etkiye bölünmesiyle elde edilen aracılık edilen etkinin oranını rapor eder. Bu örnekte, iş yoğunluğunun (X) duygusal tükenmişlik (Y) üzerindeki toplam etkisinin %57’sine bilişsel değerlendirme (M) aracılık etmektedir.
library("mediation")
fit.MX <- lm (cogapp ~ wintense, data = Paskvan)
fit.YXM <- lm(emotion ~ wintense + cogapp, data = Paskvan)
set.seed(123)
fitmed <- mediation::mediate(fit.MX, fit.YXM, treat = "wintense", mediator = "cogapp",
sims = 999, boot = TRUE, boot.ci.type = "bca")
summary(fitmed)##
## Causal Mediation Analysis
##
## Nonparametric Bootstrap Confidence Intervals with the BCa Method
##
## Estimate 95% CI Lower 95% CI Upper p-value
## ACME 0.33370 0.26244 0.42043 < 2.2e-16 ***
## ADE 0.24926 0.13986 0.37070 < 2.2e-16 ***
## Total Effect 0.58297 0.50044 0.67547 < 2.2e-16 ***
## Prop. Mediated 0.57242 0.43215 0.75245 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Sample Size Used: 803
##
##
## Simulations: 999Aynı modeli lavaan kullanarak kuralım. Önce iki regresyonu belirtiriz, ardından dolaylı etkiyi ve eğer gerekiyorsa toplam etki ile aracılık oranını tanımlarız. Her bir parametre için %95 bootstrap güven aralığı elde ederiz.
library("lavaan")
med.model <- '
emotion ~ c*wintense + b*cogapp
cogapp ~ a*wintense
ind := a*b
tot := ind+c
prop := ind/tot'
set.seed(123)
fitmedsem <- lavaan::sem(med.model, Paskvan, se = "bootstrap", bootstrap = 999)
parameterEstimates(fitmedsem, zstat = FALSE, pvalue = FALSE, boot.ci.type = "bca.simple")[c(7,1,8,9),]Sonuç: dolaylı etki (ab), doğrudan etki (c’), toplam etki ve aracılık oranını, mediate paketinin çıktısıyla aynı sırada olacak şekilde dışarı aktarır. Bunların tamamı için bootstrap yöntemiyle hesaplanmış güven aralıkları (CI) elde ederiz. Bu sonuçlar, bootstrap işleminin rastgele doğasından kaynaklanan çok küçük sapmalar dışında, mediation paketinden elde edilen sonuçlarla eşleşmektedir.
Öğrenme Günlüğü
Bu haftaki dersimize önceki hafta öğrendiğimiz yol analizi konusunu tekrar ederek başladık. Sobel testinin modası geçmiş olduğunu onun yerine bootstrapping yönteminin tavsiye edildiğini öğrendim. Ödev için iki farklı kitapta bulunan farklı veri setleri üzerinde uygulamalar yaptım ve konuyu biraz daha pekiştirmiş oldum. Aslında genel olarak konuyu anladım ama ara ara kafam çok karıştı, sanırım öğrenilen bu kadar fazla bilgiyi sindirmek için biraz zamana ihtiyacım var. Modern Psychometrics With R kitabının anlatımını çok beğendim, gerçekten tane tane ve başarılı.