Hafta 8 - OLC733 Öğrenme Günlüğü
Tugay Kaçak
2025-04-19
library(haven)
library(naniar)
library(dplyr)
library(knitr)
library(ggplot2)
library(GGally) ## değişkenler arasındaki ilişkiyi inceleme
library(lavaan) ## Yem analizleri
library(broom) ## sonuçların düzgün elde edilmesi
library(DT) ## interaktif tablolar
library(dplyr) ## veri düzenleme
library(knitr) ## rapor ve tablo oluşturma
library(semPlot) ## görselleştirme
library(semptools) ## anlamlı katsayaları gösterme
library(semoutput) ## çıktıları tablo formatında verir.
## diğer paketlerden farklı olarak devtools::install_github("dr-JT/semoutput") komutu ile yüklenmektedir
library(stargazer) ## regresyon çıktılarını düzenler
library(multilevel) ## sobel testi
library(bda) ## sobel testi-anlamlılık değeri
library(mediation)
library(gvlma) ## varsayım testi
library(rockchalk)
library(tidyr)#Bu kod baya işlevsel oldu
setwd(dirname(rstudioapi::getActiveDocumentContext()$path))1 Giriş
Bu derste yol analizine odaklandık.
- Model-veri uyumu indeksleri
- \(\chi^2\), TLI, CFI, RMSEA, SRMR, p-değeri, df
- Lavaan uygulaması ve grafikler
- Model tanımlama
- Değişken tipleri ?, confounding, suppressor gibi gibi
- Aracılık Modelleri, tekli, paralel ve serial
- Varsayımları
- Sobel testi ve Mediation paketi ile karşılaştırıldı.
2 Yol Analizi
Yol analizi için ICILS2023 Fransa öğrenci veri setinde yer alan “bilişim teknolojilerinin genel kullanımına yönelik özyeterlik algıları indeksi” (S_GENCLASS), “gelecekte bilişim teknolojilerini kullanma isteği indeksi” (S_ICTFUT), “bilişim teknolojilerini güvenli kullanma indeksi” (S_LRNSAFE) ve “bilişim teknolojilerine yönelik genel özyeterlik indeksi” (S_GENEFF) değişkenlerini kullandım. Bu değişkenler için halihazırda kodlanmış 8, 9, 998….. veya 999…. gibi “Not Reached” ya da “Omit” kodlu etiketleri “NA” olacak şekilde manipule ettim.
load("BSGFRAI3.Rdata")
data <- BSGFRAI3 %>%
dplyr::select(S_GENCLASS,S_ICTFUT,S_GENEFF,S_LRNSAFE) %>%
rename(
SINIFGENEL = S_GENCLASS,
GENELOZ = S_GENEFF,
GELECEK = S_ICTFUT,
GUVENLI = S_LRNSAFE
) %>%
expss::drop_var_labs() %>%
mutate(across(everything(), ~ ifelse(. %in% c(8, 9, 998.000000000,999.00000), NA, .)))
knitr::kable(miss_var_summary(data),caption = "Kayıp Verilerin Özeti")| variable | n_miss | pct_miss |
|---|---|---|
| GELECEK | 481 | 13.0 |
| GENELOZ | 314 | 8.50 |
| SINIFGENEL | 276 | 7.47 |
| GUVENLI | 173 | 4.68 |
En çok kayıp verinin GELECEK değişkeninde olduğunu görüyoruz. Bu değişkende kayıp veri, tüm verinin %13’ünü oluşturmaktadır. Kayıp verinin mekanizmasının MCAR’dan farklılaşıp farklılaşmadığını Little’ın MCAR testi ile inceleyelim.
naniar::mcar_test(data)## # A tibble: 1 × 4
## statistic df p.value missing.patterns
## <dbl> <dbl> <dbl> <int>
## 1 76.3 28 0.00000238 16MCAR testi sonucunda kayıp veri mekanizmasının MCAR’dan istatistiksel olarak anlamlı düzeyde farklılaştığı yorumu yapılabilir. Sürekli veri setleri için uygun olabilecek atama yöntemlerinden biri olarak EM ataması yaparak veri setindeki kayıp verileri doldurdum.
data <- missMethods::impute_EM(data, stochastic = F,maxits = 100,verbose = T)## The missing values of following rows were imputed with (parts of) mu: 173, 186, 187, 202, 206, 383, 386, 406, 413, 414, 456, 458, 463, 513, 562, 734, 803, 804, 806, 810, 814, 829, 838, 867, 1003, 1021, 1048, 1055, 1056, 1336, 1434, 1512, 1520, 1528, 1559, 1641, 1644, 1647, 1654, 1656, 1715, 1719, 1732, 1812, 1829, 1848, 1851, 1855, 1896, 2002, 2004, 2040, 2046, 2047, 2191, 2391, 2445, 2519, 2531, 2533, 2543, 2544, 2568, 2571, 2584, 2596, 2708, 2836, 2848, 2906, 2912, 2926, 2942, 2946, 2947, 2968, 2969, 3008, 3043, 3044, 3050, 3051, 3052, 3054, 3055, 3059, 3196, 3210, 3277, 3406, 3416, 3425, 3436, 3448, 3459, 3467, 3521, 3551, 3653, 3663, 3673, 3692, 36943 Betimsel İstatistikler
Veri setine ilişkin betimsel istatistikleri inceleyelim. Tek değişkenli normallik açısından değişkenler incelendiğinde, değişkenlerin çarpıklık katsayılarının [-1,1] arasında yer aldığı ve normal dağılımdan sapmadığı görülmektedir. Bununla birlikte basıklık katsayıları incelendiğinde ise SINIFGENEL ve GELECEK değişkenlerinin normal dağılımdan saptığı yorumu yapılabilir.
kable(psych::describe(data,fast = T),digits = 2)| vars | n | mean | sd | median | min | max | range | skew | kurtosis | se | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| SINIFGENEL | 1 | 3694 | 48.48 | 7.86 | 49.20 | 26.62 | 78.12 | 51.50 | -0.08 | 2.35 | 0.13 |
| GELECEK | 2 | 3694 | 49.08 | 9.16 | 48.68 | 23.17 | 72.64 | 49.47 | 0.25 | 1.21 | 0.15 |
| GENELOZ | 3 | 3694 | 51.69 | 8.95 | 50.56 | 14.31 | 71.57 | 57.26 | 0.27 | 0.67 | 0.15 |
| GUVENLI | 4 | 3694 | 46.28 | 9.00 | 46.09 | 22.86 | 67.29 | 44.42 | 0.23 | 0.56 | 0.15 |
hist(data$SINIFGENEL)
hist(data$GELECEK)
hist(data$GENELOZ)
hist(data$GUVENLI)
4 Çok Değişkenli Normallik Varsayımının İncelenmesi
Çok değişkenli normallik Mardia’nın Testi ile incelenmiştir.
library(mvnormalTest)## Warning: package 'mvnormalTest' was built under R version 4.4.2mardia(data)## $mv.test
## Test Statistic p-value Result
## 1 Skewness 299.548 0 NO
## 2 Kurtosis 43.7068 0 NO
## 3 MV Normality <NA> <NA> NO
##
## $uv.shapiro
## W p-value UV.Normality
## SINIFGENEL 0.9234 0 No
## GELECEK 0.9545 0 No
## GENELOZ 0.9645 0 No
## GUVENLI 0.9655 0 NoMardia’nın testi sonucunda veri setinin çok değişkenli normallik varsayımını sağlamadığı görülmektedir.
5 Uç Değerleri İnceleyelim
5.1 Tek Değişkenli Uç Değerler
Tek değişkenli uç değerler z-puanına göre değerlendirilmiş olup [-3,3] arasında yer almayan bireyler analizden çıkarılmıştır.
data_clean <- data %>%
mutate(across(
c(SINIFGENEL, GENELOZ, GELECEK, GUVENLI),
~ as.numeric(scale(.x)),
.names = "z{.col}"
)) %>%
filter(if_all(starts_with("z"), ~ between(.x, -3.001, 3.001)))5.2 Çok Değişkenli Uç Değerler
Çok değişkenli uç değerler Mahalanobis Uzaklığı kullanılarak veri setinden çıkarılmıştır.
#mahalanobis
data_clean <- data_clean[,1:4]
md <- mahalanobis(data_clean, center = colMeans(data_clean), cov = cov(data_clean))
alpha <- .001
cutoff <- (qchisq(p = 1 - alpha, df = ncol(data_clean)))
ucdegerler <- which(md > cutoff)
data_clean <- data_clean[-ucdegerler, ]6 Başlangıç Modelinin Kurulması
Bu çalışmada yol analizi için başlangıç modelinde öğrencilerin SINIFGENEL indeksleri yordanan değişken, GENELOZ indeksi yordayıcı değişken, GELECEK ve GUVENLI indeksleri ise aracı değişkenler olarak tanımlanmıştır.
7 Değişkenler Arasındaki Korelasyonların İncelenmesi
data <- data_clean
library(ggcorrplot)
corp <- cor_pmat(data)
ggcorrplot(cor(data),lab = T, method = "circle", title = "Değişkenler Arasındaki Korelasyon", p.mat = corp)
Veri setinde alan değişkenler arasındak ilişkilerin tamamı istatistiksel olarak anlamlı olup 0.06 - 0.27 arasında değişmektedir.
8 Çoklu Bağlantı Varsayımının İncelenmesi
library(olsrr)
model = lm(SINIFGENEL ~ GENELOZ + GELECEK + GUVENLI, data = data)
ols_vif_tol(model = lm(SINIFGENEL ~ GENELOZ + GELECEK + GUVENLI, data = data)) %>% kable(.,digits = 2)| Variables | Tolerance | VIF |
|---|---|---|
| GENELOZ | 0.91 | 1.10 |
| GELECEK | 0.93 | 1.08 |
| GUVENLI | 0.97 | 1.03 |
VIF ve TOL değerleri incelendiğinde veri setinde çoklu bağlantı problemi olmadığı yorumu yapılabilir.
9 Modelin Tanımlanması
library(lavaan)
library(semPlot)
yol_model <- 'GUVENLI ~ GENELOZ
GELECEK ~ GENELOZ
SINIFGENEL ~ GENELOZ + GELECEK + GUVENLI'
yol_fit <- sem(yol_model, data)
semPaths(yol_fit,rotation=1, curvePivot = TRUE,
sizeMan = 6, sizeInt = 1,
sizeLat = 5,
edge.label.cex = 1.8,
pastel=TRUE,
nCharNodes = 0, nCharEdges = 0)
# summary(yol_fit, fit.measures = TRUE)
fitmeasures(yol_fit,fit.measures = c("chisq" ,"df" ,"pvalue",
"cfi","tli","rmsea",
"rmsea.ci.lower","rmsea.ci.upper"
,"srmr"))## chisq df pvalue cfi tli
## 9.421 1.000 0.002 0.982 0.894
## rmsea rmsea.ci.lower rmsea.ci.upper srmr
## 0.048 0.024 0.078 0.015Model veri uyumları incelendiğinde \(\chi^2\) değeri 9.421, p-değeri 0.002, RMSEA 0.048[0.024-0.078], CFI 0.982, TLI 0.894 ve SRMR 0.015 olarak bulunmuştur. Bu sonuçlar modelin veri setine iyi uyum sağladığı yorumu yapılabilir.
10 Artıkları Değerlendir
resid(yol_fit, type='normalized')## $type
## [1] "normalized"
##
## $cov
## GUVENL GELECE SINIFG GENELO
## GUVENLI 0.000
## GELECEK 2.905 0.000
## SINIFGENEL 0.105 0.384 0.020
## GENELOZ 0.000 0.000 0.000 0.00011 Modifikasyon İndekslerinin Değerlendirilmesi
Modifikasyon indeksleri değerlendirildiğinde herhangi bir düzeltmenin \(\chi^2\) değerini 10’un üzerinde düzeltmeyeceği görülmüştür. Ayrıca eklenecek herhangi bir modifikasyon ise modelin serbestlik derecesinin 0 olmasına, yani kestirimlerin yapılamamasın neden olacaktır. Bu nedenle herhangi bir model yeniden tanımlama işlemi yapılmamıştır.
modindices(yol_fit, sort = TRUE)## lhs op rhs mi epc sepc.lv sepc.all sepc.nox
## 14 GUVENLI ~ SINIFGENEL 9.409 1.762 1.762 1.425 1.425
## 16 GELECEK ~ SINIFGENEL 9.409 0.469 0.469 0.378 0.378
## 10 GUVENLI ~~ GELECEK 9.409 3.801 3.801 0.051 0.051
## 13 GUVENLI ~ GELECEK 9.409 0.052 0.052 0.052 0.052
## 15 GELECEK ~ GUVENLI 9.409 0.050 0.050 0.050 0.050
## 19 GENELOZ ~ SINIFGENEL 0.000 0.000 0.000 0.000 0.000
## 17 GENELOZ ~ GUVENLI 0.000 0.000 0.000 0.000 0.000
## 18 GENELOZ ~ GELECEK 0.000 0.000 0.000 0.000 0.000yol_model_v1 <-
' GUVENLI ~ GENELOZ
GELECEK ~ GENELOZ
SINIFGENEL ~ GENELOZ + GELECEK + GUVENLI'
yol_fit_v1 <- sem(yol_model_v1, data)
semPaths(yol_fit_v1,rotation=2, curvePivot = TRUE,
sizeMan = 12, sizeInt = 1,
sizeLat = 4,
edge.label.cex = 1.8,
pastel=TRUE,
nCharNodes = 0, nCharEdges = 0)
p_pa <- semPaths(yol_fit_v1, whatLabels = "est",
sizeMan = 10,
edge.label.cex = 1.15,
style = "lisrel",layout = "tree" ,
nCharNodes = 0, nCharEdges = 0)
p_pa_2 <- semptools::mark_sig(p_pa, yol_fit_v1)
plot(p_pa_2)
12 Manuel Tablo
m <- matrix(c(NA, NA, "GUVENLI", NA, NA,
NA, NA, NA, NA, NA,
"GENELOZ", NA, NA, NA, "SINIFGENEL",
NA, NA, NA, NA, NA,
NA, NA, "GELECEK", NA, NA
), byrow = TRUE, 5, 5)
p_pa <- semPaths(yol_fit_v1, whatLabels = "est",
sizeMan = 10,
edge.label.cex = 1.15,
style = "ram",
nCharNodes = 0, nCharEdges = 0,
layout=m)
p_pa_2 <- semptools::mark_sig(p_pa, yol_fit_v1)
plot(p_pa_2)
Görsel 2 incelendiğinde tüm yolların istatistiksel olarak anlamlı olduğu
görülmektedir.
13 Ozet Bilgi
summary(yol_fit_v1, fit.measures = TRUE, standardized = T)## lavaan 0.6-19 ended normally after 1 iteration
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 8
##
## Number of observations 3616
##
## Model Test User Model:
##
## Test statistic 9.421
## Degrees of freedom 1
## P-value (Chi-square) 0.002
##
## Model Test Baseline Model:
##
## Test statistic 482.279
## Degrees of freedom 6
## P-value 0.000
##
## User Model versus Baseline Model:
##
## Comparative Fit Index (CFI) 0.982
## Tucker-Lewis Index (TLI) 0.894
##
## Loglikelihood and Information Criteria:
##
## Loglikelihood user model (H0) -38051.931
## Loglikelihood unrestricted model (H1) -38047.221
##
## Akaike (AIC) 76119.863
## Bayesian (BIC) 76169.408
## Sample-size adjusted Bayesian (SABIC) 76143.988
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.048
## 90 Percent confidence interval - lower 0.024
## 90 Percent confidence interval - upper 0.078
## P-value H_0: RMSEA <= 0.050 0.475
## P-value H_0: RMSEA >= 0.080 0.041
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.015
##
## Parameter Estimates:
##
## Standard errors Standard
## Information Expected
## Information saturated (h1) model Structured
##
## Regressions:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## GUVENLI ~
## GENELOZ 0.180 0.017 10.599 0.000 0.180 0.174
## GELECEK ~
## GENELOZ 0.277 0.017 16.629 0.000 0.277 0.267
## SINIFGENEL ~
## GENELOZ 0.045 0.014 3.134 0.002 0.045 0.054
## GELECEK 0.030 0.014 2.147 0.032 0.030 0.037
## GUVENLI 0.107 0.013 7.922 0.000 0.107 0.132
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .GUVENLI 75.922 1.786 42.521 0.000 75.922 0.970
## .GELECEK 73.148 1.720 42.521 0.000 73.148 0.929
## .SINIFGENEL 49.929 1.174 42.521 0.000 49.929 0.974standardizedSolution(yol_fit_v1)## lhs op rhs est.std se z pvalue ci.lower ci.upper
## 1 GUVENLI ~ GENELOZ 0.174 0.016 10.845 0.000 0.142 0.205
## 2 GELECEK ~ GENELOZ 0.267 0.015 17.567 0.000 0.237 0.296
## 3 SINIFGENEL ~ GENELOZ 0.054 0.017 3.140 0.002 0.020 0.088
## 4 SINIFGENEL ~ GELECEK 0.037 0.017 2.148 0.032 0.003 0.070
## 5 SINIFGENEL ~ GUVENLI 0.132 0.017 7.989 0.000 0.100 0.164
## 6 GUVENLI ~~ GUVENLI 0.970 0.006 174.529 0.000 0.959 0.981
## 7 GELECEK ~~ GELECEK 0.929 0.008 114.868 0.000 0.913 0.945
## 8 SINIFGENEL ~~ SINIFGENEL 0.974 0.005 187.787 0.000 0.964 0.984
## 9 GENELOZ ~~ GENELOZ 1.000 0.000 NA NA 1.000 1.00014 Std. Artık
# Bunu güncelle
resid(yol_fit_v1)## $type
## [1] "raw"
##
## $cov
## GUVENL GELECE SINIFG GENELO
## GUVENLI 0.000
## GELECEK 3.801 0.000
## SINIFGENEL 0.112 0.406 0.024
## GENELOZ 0.000 0.000 0.000 0.00015 Kestirim
parameterEstimates(yol_fit_v1,standardized = T)## lhs op rhs est se z pvalue ci.lower ci.upper std.lv
## 1 GUVENLI ~ GENELOZ 0.180 0.017 10.599 0.000 0.147 0.213 0.180
## 2 GELECEK ~ GENELOZ 0.277 0.017 16.629 0.000 0.244 0.310 0.277
## 3 SINIFGENEL ~ GENELOZ 0.045 0.014 3.134 0.002 0.017 0.074 0.045
## 4 SINIFGENEL ~ GELECEK 0.030 0.014 2.147 0.032 0.003 0.056 0.030
## 5 SINIFGENEL ~ GUVENLI 0.107 0.013 7.922 0.000 0.080 0.133 0.107
## 6 GUVENLI ~~ GUVENLI 75.922 1.786 42.521 0.000 72.422 79.421 75.922
## 7 GELECEK ~~ GELECEK 73.148 1.720 42.521 0.000 69.776 76.519 73.148
## 8 SINIFGENEL ~~ SINIFGENEL 49.929 1.174 42.521 0.000 47.628 52.231 49.929
## 9 GENELOZ ~~ GENELOZ 72.945 0.000 NA NA 72.945 72.945 72.945
## std.all std.nox
## 1 0.174 0.020
## 2 0.267 0.031
## 3 0.054 0.006
## 4 0.037 0.037
## 5 0.132 0.132
## 6 0.970 0.970
## 7 0.929 0.929
## 8 0.974 0.974
## 9 1.000 72.94516 Manuel Indirect Hesaplama
yol_model <- '
# Paths from predictors to mediators
GELECEK ~ a1*GENELOZ
GUVENLI ~ a2*GENELOZ
# Paths from mediators to outcome
SINIFGENEL ~ b1*GELECEK + b2*GUVENLI
# Direct effects to outcome
SINIFGENEL ~ c1*GENELOZ
# Indirect effects: GENELOZ
ind_GELECEK := a1 * b1
ind_GUVENLI := a2 * b2
total_indirect_GENELOZ := ind_GELECEK + ind_GUVENLI
total_GENELOZ := c1 + total_indirect_GENELOZ
gelecek_ratio := ind_GELECEK / total_GENELOZ
guvenli_ratio := ind_GUVENLI / total_GENELOZ
ratio := total_indirect_GENELOZ / total_GENELOZ
'
fsem1 <- sem(yol_model,data)
standardizedsolution(fsem1)## lhs op rhs
## 1 GELECEK ~ GENELOZ
## 2 GUVENLI ~ GENELOZ
## 3 SINIFGENEL ~ GELECEK
## 4 SINIFGENEL ~ GUVENLI
## 5 SINIFGENEL ~ GENELOZ
## 6 GELECEK ~~ GELECEK
## 7 GUVENLI ~~ GUVENLI
## 8 SINIFGENEL ~~ SINIFGENEL
## 9 GENELOZ ~~ GENELOZ
## 10 ind_GELECEK := a1*b1
## 11 ind_GUVENLI := a2*b2
## 12 total_indirect_GENELOZ := ind_GELECEK+ind_GUVENLI
## 13 total_GENELOZ := c1+total_indirect_GENELOZ
## 14 gelecek_ratio := ind_GELECEK/total_GENELOZ
## 15 guvenli_ratio := ind_GUVENLI/total_GENELOZ
## 16 ratio := total_indirect_GENELOZ/total_GENELOZ
## label est.std se z pvalue ci.lower ci.upper
## 1 a1 0.267 0.015 17.567 0.000 0.237 0.296
## 2 a2 0.174 0.016 10.845 0.000 0.142 0.205
## 3 b1 0.037 0.017 2.148 0.032 0.003 0.070
## 4 b2 0.132 0.017 7.989 0.000 0.100 0.164
## 5 c1 0.054 0.017 3.140 0.002 0.020 0.088
## 6 0.929 0.008 114.868 0.000 0.913 0.945
## 7 0.970 0.006 174.529 0.000 0.959 0.981
## 8 0.974 0.005 187.787 0.000 0.964 0.984
## 9 1.000 0.000 NA NA 1.000 1.000
## 10 ind_GELECEK 0.010 0.005 2.131 0.033 0.001 0.019
## 11 ind_GUVENLI 0.023 0.004 6.402 0.000 0.016 0.030
## 12 total_indirect_GENELOZ 0.033 0.006 5.641 0.000 0.021 0.044
## 13 total_GENELOZ 0.087 0.016 5.269 0.000 0.055 0.119
## 14 gelecek_ratio 0.112 0.057 1.976 0.048 0.001 0.224
## 15 guvenli_ratio 0.264 0.063 4.205 0.000 0.141 0.387
## 16 ratio 0.376 0.096 3.930 0.000 0.189 0.564Modele ilişkin doğrudan ve dolaylı etkilerin tamamının anlamlı olduğu görülmektedir. Toplam dolaylı etkiler, doğrudan etkilerin %37.6’sını açıklamaktadır. Ayrıca GENELOZ –> GELECEK –> SINIFGENEL yolu toplam etkilerin %11.2’sini, GENELOZ –> GUVENLI –> SINIFGENEL yolu ise toplam etkilerin %26.4’ünü açıklamaktadır. Standartlaştırılmış beta katsayıları incelendiğinde GUVENLI kullanım aracı değişkeninin GELECEK aracı değişkenine göre daha güçlü bir aracı değişken olduğu görülmektedir.
16.0.1 Durbin-Watson
Hataların ilişkililiğini test eder.
library(car)
model<-lm(SINIFGENEL~ GENELOZ + GELECEK, data=data)
durbinWatsonTest(model)## lag Autocorrelation D-W Statistic p-value
## 1 0.09946575 1.800187 0
## Alternative hypothesis: rho != 0Durbin-Watson testine göre değişkenler arasında otokorelasyon yoktur/vardır. (?) Bu istatistiği daha önce mi incelemeliydik ?
16.0.2 Yöntem 1 - Baron & Kenny
GENELOZ –> GELECEK –> SINIFGENEL yolunun aracılık etkisini incelemek için Baron & Kenny yöntemi kullanılmıştır. Bu yöntemde 4 temel adım bulunmaktadır.
library(stargazer) # Kullanışlı regresyon tabloları
data <- data_clean
#1. toplam etki
fit <- lm(SINIFGENEL ~ GENELOZ, data=data)
#2. YOL A (X on M)
fita <- lm(GELECEK ~ GENELOZ, data=data)
#3. YOL B' (M on Y, X kontrol edildiğinde )
fitb <- lm(SINIFGENEL ~ GELECEK + GENELOZ, data=data)
#4. YOL C - ters yol (Y on X, M kontrol edildiğinde)
fitc <- lm(GENELOZ ~ SINIFGENEL + GELECEK, data=data)
# Ö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:
## -------------------------------------------------------------------------------------------------
## SINIFGENEL GELECEK SINIFGENEL GENELOZ
## (1) (2) (3) (4)
## ---------------------------------------------------------------------------------------------------------------------
## SINIFGENEL 0.08***
## (0.02)
##
## GELECEK 0.04** 0.25***
## (0.01) (0.02)
##
## GENELOZ 0.07*** 0.28*** 0.06***
## (0.01) (0.02) (0.01)
##
## Constant 44.57*** 34.72*** 43.35*** 35.27***
## (0.73) (0.87) (0.87) (1.16)
##
## ---------------------------------------------------------------------------------------------------------------------
## Observations 3,616 3,616 3,616 3,616
## R2 0.01 0.07 0.01 0.08
## Adjusted R2 0.01 0.07 0.01 0.08
## Residual Std. Error 7.14 (df = 3614) 8.56 (df = 3614) 7.13 (df = 3613) 8.21 (df = 3613)
## F Statistic 27.42*** (df = 1; 3614) 276.35*** (df = 1; 3614) 16.93*** (df = 2; 3613) 148.45*** (df = 2; 3613)
## =====================================================================================================================
## Note: *p<0.1; **p<0.05; ***p<0.01#Sobel Test
library(multilevel)
sobel(data$GENELOZ, data$GELECEK, data$SINIFGENEL)## $`Mod1: Y~X`
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 44.56950967 0.72777719 61.240597 0.000000e+00
## pred 0.07275435 0.01389302 5.236755 1.727491e-07
##
## $`Mod2: Y~X+M`
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 43.35231959 0.87213565 49.708230 0.0000000000
## pred 0.06304789 0.01440368 4.377207 0.0000123612
## med 0.03505246 0.01386341 2.528415 0.0115001236
##
## $`Mod3: M~X`
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 34.7248146 0.87259053 39.79509 1.592835e-287
## pred 0.2769124 0.01665746 16.62393 7.407413e-60
##
## $Indirect.Effect
## [1] 0.009706461
##
## $SE
## [1] 0.0038831
##
## $z.value
## [1] 2.499668
##
## $N
## [1] 3616library(bda)
attach(data)
mediation.test(iv = c(GENELOZ),mv = c(GELECEK),dv = c(SINIFGENEL))## Sobel Aroian Goodman
## z.value 2.49966822 2.4952596 2.50410030
## p.value 0.01243097 0.0125865 0.0122763216.0.3 Illness Veri Setinin Baron & Kenny Yaklaşımı ile İncelenmesi
illness <- read.table("illness.dat")
colnames(illness) <- c("form", "stres", "hastalik", "egzersiz", "dayaniklilik")
library(stargazer) # Kullanışlı regresyon tabloları
#1. toplam etki
fit <- lm(hastalik ~ egzersiz, data=illness)
#2. YOL A (X on M)
fita <- lm(form ~ egzersiz, data=illness)
#3. YOL B' (M on Y, X kontrol edildiğinde )
fitb <- lm(hastalik ~ form + egzersiz, data=illness)
#4. YOL C - ters yol (Y on X, M kontrol edildiğinde)
fitc <- lm(egzersiz ~ form + hastalik, data=illness)
# Ö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:
## ---------------------------------------------------------------------------------------
## hastalik form hastalik egzersiz
## (1) (2) (3) (4)
## -----------------------------------------------------------------------------------------------------------
## form -0.54*** 0.70***
## (0.08) (0.09)
##
## egzersiz -0.06 0.21*** 0.05
## (0.04) (0.03) (0.05)
##
## hastalik 0.06
## (0.06)
##
## Constant 182.68*** 237.11*** 310.96*** 32.43
## (11.18) (6.55) (21.98) (29.90)
##
## -----------------------------------------------------------------------------------------------------------
## Observations 400 400 400 400
## R2 0.01 0.14 0.11 0.14
## Adjusted R2 0.003 0.14 0.10 0.14
## Residual Std. Error 61.65 (df = 398) 36.15 (df = 398) 58.54 (df = 397) 65.02 (df = 397)
## F Statistic 2.03 (df = 1; 398) 63.80*** (df = 1; 398) 23.34*** (df = 2; 397) 32.50*** (df = 2; 397)
## ===========================================================================================================
## Note: *p<0.1; **p<0.05; ***p<0.01EGZERSİZ –> FORM –> HASTALIK yolu Baron & Kenny yöntemi ile test edildiğinde aracı değişken olan FORMUN anlamlı bir aracı değişken olabileceğine ilişkin beklenti gelişir. Bu çıkarım, bağımlı değişken HASTALIK ve yordayıcı değişken form, egzersizin ise kontrol edildiği durumda elde edilen yolun anlamlı olması ve (\(\Beta = -0.54, p<0.001\)), ters yol için aracı değişken olan FORMun anlamlı olması (\(\Beta = 0.70, p<0.001\)) bu değişkenin anlamlı bir aracı etkiye sahip olabileceğine işaret eder.
16.0.4 Yöntem 2 - Mediation Bootstrap
library(mediation)
library(gvlma)
fitM <- lm(GELECEK ~ GENELOZ, data=data)
fitY <- lm(SINIFGENEL ~ GENELOZ + GELECEK, data=data)
gvlma(fitM) # veri birçok varsayımı sağlamıyor##
## Call:
## lm(formula = GELECEK ~ GENELOZ, data = data)
##
## Coefficients:
## (Intercept) GENELOZ
## 34.7248 0.2769
##
##
## 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 229.3034 0.0000000 Assumptions NOT satisfied!
## Skewness 7.4805 0.0062372 Assumptions NOT satisfied!
## Kurtosis 208.7870 0.0000000 Assumptions NOT satisfied!
## Link Function 12.1096 0.0005016 Assumptions NOT satisfied!
## Heteroscedasticity 0.9264 0.3358125 Assumptions acceptable.gvlma(fitY)##
## Call:
## lm(formula = SINIFGENEL ~ GENELOZ + GELECEK, data = data)
##
## Coefficients:
## (Intercept) GENELOZ GELECEK
## 43.35232 0.06305 0.03505
##
##
## 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 560.9783 0.00000 Assumptions NOT satisfied!
## Skewness 178.1040 0.00000 Assumptions NOT satisfied!
## Kurtosis 376.8062 0.00000 Assumptions NOT satisfied!
## Link Function 0.8936 0.34449 Assumptions acceptable.
## Heteroscedasticity 5.1746 0.02292 Assumptions NOT satisfied!fitMed <- mediate(fitM, fitY, treat=c("GENELOZ"), mediator=c("GELECEK"))
summary(fitMed)##
## Causal Mediation Analysis
##
## Quasi-Bayesian Confidence Intervals
##
## Estimate 95% CI Lower 95% CI Upper p-value
## ACME 0.00981 0.00223 0.02 0.012 *
## ADE 0.06272 0.03330 0.09 <2e-16 ***
## Total Effect 0.07253 0.04321 0.10 <2e-16 ***
## Prop. Mediated 0.13333 0.02870 0.29 0.012 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Sample Size Used: 3616
##
##
## Simulations: 1000#Bootstrap
fitMedBoot <- mediate(fitM, fitY, boot=TRUE, sims=999, treat="GENELOZ", mediator="GELECEK")## Running nonparametric bootstrapsummary(fitMedBoot)##
## Causal Mediation Analysis
##
## Nonparametric Bootstrap Confidence Intervals with the Percentile Method
##
## Estimate 95% CI Lower 95% CI Upper p-value
## ACME 0.00971 0.00209 0.02 0.016 *
## ADE 0.06305 0.03134 0.09 <2e-16 ***
## Total Effect 0.07275 0.04204 0.10 <2e-16 ***
## Prop. Mediated 0.13341 0.02470 0.30 0.016 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Sample Size Used: 3616
##
##
## Simulations: 999Bootstrap, Baron & Kenny yöntemi ile GENELOZ –> GELECEK –> SINIFGENEL yolu incelendiğinde Doğrudan etkilerin, toplam etkinni ve aracılık etkisinin anlamlı bir yordayıcı olduğu görülmektedir. Bootstrap yöntemi için dolaylı etki toplam etkinin %13.3’ünü açıklamakta olup, Baron & Kenny yöntemi için (??) hangi kısımı almalıyız. axb/c+(axb) olarak hesaplamak mı gerekirdi emin olamadım. Sonuç olarak GELECEK değişkeninin, GENELOZ ve SINIFGENEL arasındaki ilişkiye aracılık ettiği gözlenmiştir.
Not: 2 günde, toplam 3 versiyon halinde tamamlayabildim. 5-6 saatimi aldı.