3. ve 5. Haftaların Öğrenme Günlüğü (Aracı Etki Analizi Denemesi)
Samet EKER
20 MART 2025
R ile Çok Değişkenli İstatistik ve Psikometri dersi 3. ve 5. hafta için yapacağım öğrenme günlüğünde PISA verileri ile bir aracı etki analizi deneyeceğim.
Verinin Yüklenmesi ve Temizlenmesi
Gerekli paketlerin yüklenmesi
library(tidyverse)
library(haven)
library(sjPlot)
library(sjmisc)
library(sjlabelled)
library(sjstats)
library(vtable)
library(kableExtra)
library(dplyr)
library(knitr)
library(psych)
library(Hmisc)
library(naniar)
library(ggplot2)
library(mice)
library(tidyverse)
library(stevemisc)
library(knitr)
library(haven)
library(summarytools)
library(outliers)
library(ggplot2)
library(plotly)
library(ggpmisc)
library(psych)
library(sur)
library(moments)
library(corrplot)
library(olsrr)
library(magrittr)
library(skimr)
library(DT)
library(robustbase)
library(car)
library(MVN)
library(mediation)
library(lm.beta)
library(interactions)Veri dosyasının R ortamına yüklenmesi
Resmi siteden .sav uzantılı indirilen ham veriler.
İndirilen dosyada 1278 değişken, 613744 gözlem bulunmaktadır.
Türkiye verilerinin seçilmesi
Türkiye verilerinde 1278 değişken, 7250 gözlem bulunmaktadır.
Bu değişkenlerden sadece bazılarını kullanacağım. Kullanacağım değişkenler: Zorbalığa Maruz Kalma (BULLIED), Okulda Aidiyet Duygusu (BELONG) Güvende Hissetme (FEELSAFE)
library(dplyr)
pisa_turkiye <- pisa_turkiye %>% dplyr::select(BULLIED, BELONG, FEELSAFE)
head(pisa_turkiye, 5)## # A tibble: 5 × 3
## BULLIED BELONG FEELSAFE
## <dbl+lbl> <dbl+lbl> <dbl+lbl>
## 1 -0.267 -0.907 -0.0323
## 2 1.41 -1.92 -2.19
## 3 -1.23 -0.340 -0.756
## 4 0.927 -0.326 -0.756
## 5 -1.23 -0.313 -0.756Seçimden sonra 3 değişken, 7250 gözlem bulunmaktadır.
Veri Seti Ön İnceleme
## BULLIED BELONG FEELSAFE
## Min. :-1.2280 Min. :-3.2583 Min. :-2.7886
## 1st Qu.:-1.2280 1st Qu.:-0.9024 1st Qu.:-1.1042
## Median :-0.2670 Median :-0.3485 Median :-0.7560
## Mean :-0.1990 Mean :-0.3077 Mean :-0.4216
## 3rd Qu.: 0.6017 3rd Qu.: 0.1773 3rd Qu.: 0.4417
## Max. : 4.6939 Max. : 2.7562 Max. : 1.1687
## NA's :29 NA's :40 NA's :26sumtable fonksiyonu ile değişkenlerin frekans
tablolarını inceleyelim.
| Variable | NotNA | Min | Max |
|---|---|---|---|
| BULLIED | 7221 | -1.2 | 4.7 |
| BELONG | 7210 | -3.3 | 2.8 |
| FEELSAFE | 7224 | -2.8 | 1.2 |
Detaylı bilgi sunan skimr paketini kullanarak veri setini inceleyelim.
| skim_type | skim_variable | n_missing | complete_rate | numeric.mean | numeric.sd | numeric.p0 | numeric.p25 | numeric.p50 | numeric.p75 | numeric.p100 | numeric.hist |
|---|---|---|---|---|---|---|---|---|---|---|---|
| numeric | BULLIED | 29 | 0.9960000 | -0.1989759 | 1.0513621 | -1.2280 | -1.2280 | -0.2670 | 0.6017 | 4.6939 | ▇▅▂▁▁ |
| numeric | BELONG | 40 | 0.9944828 | -0.3076687 | 0.9702173 | -3.2583 | -0.9024 | -0.3485 | 0.1773 | 2.7562 | ▁▃▇▂▁ |
| numeric | FEELSAFE | 26 | 0.9964138 | -0.4215754 | 1.0304068 | -2.7886 | -1.1042 | -0.7560 | 0.4417 | 1.1687 | ▁▃▇▂▆ |
Kayıp Veri Analizi
Kayıp veri analizi için önce eksik verilerin özetine bakalım.
| variable | n_miss | pct_miss |
|---|---|---|
| BELONG | 40 | 0.552 |
| BULLIED | 29 | 0.4 |
| FEELSAFE | 26 | 0.359 |
Kayıp veriyi görselleştirmek için mice paketinin
md.pattern() fonksiyonu da oldukça kullanışlı.
## FEELSAFE BULLIED BELONG
## 7191 1 1 1 0
## 21 1 1 0 1
## 11 1 0 1 1
## 1 1 0 0 2
## 6 0 1 1 1
## 3 0 1 0 2
## 2 0 0 1 2
## 15 0 0 0 3
## 26 29 40 95Ön incelemelerde BULLIED değişkeninde 29 (% 0.40), BELONG değişkeninde 40 (% 0.55), FEELSAFE değişkeninde 26 (% 0.36) eksik gözlem olduğu görülmektedir. Bu eksik gözlemlerin tamamen rast gele dağılıp dağılmadığını incelemek için little MCAR testi yapalım.
MCAR Testi
mcar_testi <- naniar::mcar_test(pisa_turkiye)
kable(mcar_testi, caption = "BULLIED, BELONG ve FEELSAFE Değişkenlerinin Ortalamaları", align = "c") %>%
kable_styling(full_width = TRUE, position = "center")| statistic | df | p.value | missing.patterns |
|---|---|---|---|
| 8.968269 | 9 | 0.4402084 | 8 |
MCAR testi sonuçları incelendiğinde p değerinin 0.05’ten BÜYÜK olduğu görülmektedir. Bu durumda eksik verilerin tamamen rastgele (MCAR) olduğu hipotezini reddedemiyoruz. Yani, eksik verilerin belirli bir desen veya sistematik bir neden olmaksızın rastgele dağıldığını söyleyebiliriz.
MCAR varsayımı sağlandığı için, eksik veri işlemleri için liste bazlı silme (listwise deletion), ortalama ile doldurma (mean imputation) gibi basit atama yöntemleri de güvenilir olabilir. Ancak daha güvenilir sonuçlar verdiği için çoklu atama (multiple imputation) yöntemi tercih edilmelidir.
Önceki öğrenme günlüğünde multiple imputation uygulamıştım. Bu ödevde veri setindeki kayıp verilerin oranları % 5 ’ten az olduğu için listwise deletion yöntemini uygulayacağım.
Listwise Deletion
Kayıp veri var mı?
Listewise silme işlemi sonrası kayıp veri olup olmadığını kontrol edelim.
## BULLIED BELONG FEELSAFE
## 0 0 0Kayıp veri olmadığı görülmektedir.
Kaç gözlem silindi?, Kaç gözlem kaldı?
Silinen ve kalan gözlem sayılarını kontrol edelim.
## [1] "Silinen toplam gözlem sayısı: 59"## [1] "Kalan toplam gözlem sayısı: 7191"Satır bazlı silme işlemi sonrası veri seti incelendiğinde 7250 gözlemden 7191 gözlem kaldığı görülmektedir. Yani, 59 gözlem silinmiştir.
Uç Değerler
Tek Değişkenli Uç Değer Analizi
Tek değişkenli uç değer analizi için z puanlarına bakılır. PISA benzeri çok büyük veri setlerinde normalde +- 3 arası olması beklenen değerlerin, +-4 arasında olması da uygun kabul edilebilir.
Kullanılan değişkenler PISA tarafından ortalaması sıfır ve standart sapması 1 olacak şekilde standartlaştırıldığı için ölçekleme yapmadan devam edeceğim.
Uç değerlerin olup olmadığını min ve max değerlere bakarak kontrol edebiliriz.
| BULLIED | BELONG | FEELSAFE | |
|---|---|---|---|
| Min. :-1.2280 | Min. :-3.2583 | Min. :-2.7886 | |
| 1st Qu.:-1.2280 | 1st Qu.:-0.8994 | 1st Qu.:-1.1042 | |
| Median :-0.2670 | Median :-0.3485 | Median :-0.7560 | |
| Mean :-0.1987 | Mean :-0.3067 | Mean :-0.4226 | |
| 3rd Qu.: 0.6017 | 3rd Qu.: 0.1773 | 3rd Qu.: 0.4417 | |
| Max. : 4.6939 | Max. : 2.7562 | Max. : 1.1687 |
Sonuçlar incelendiğinde, veri setinde bazı uç değerler olduğu düşünülmektedir. -3, +3 aralığında olması beklenen değerlerin büyük veri setlerinde +-4 aralığında olabileceği kabul edilen bir durum olsa da BULLIED değişkeninde +4’ün üzerinde değerler olduğu anlaşılmıştır. Ancak bu değerlerin uç değer olmadığı ve anlamlı gözlemler olduğu değerlendirilmiş ve analizden çıkarılmamıştır (Zorbalığa uğramadığını bildirenlerin sayısı oldukça fazladır).
Çok Değişkenli Uç Değerler
Çok değişkenli uç değerlerin tespiti için Mahalanobis uzaklığı hesaplanır. ki-kare dağılımı kullanılarak değerlendirilebilen bu değerin gözlenme olasılığı 0.001 veya daha küçükse gözlem uç değer olarak tespit edilir.
robustbasepaketinde mahaalanobis fonksiyonu kullanılır.- Sonuçlar kritik ki-kare değeri ile karşılaştırılır.
DTpaketi ile de görselleştirme yapıldığında uç değerlerin kolaylıkla tespit edilebileceği anlaşılmıştır.
library(robustbase)
cov_robust <- covMcd(pisa_turkiye_listwise)
cdud <- mahalanobis(pisa_turkiye_listwise,
center = cov_robust$center,
cov = cov_robust$cov)
alpha <- 0.001
cutoff <- qchisq(p = 1 - alpha, df = ncol(pisa_turkiye_listwise))
outliers <- which(cdud > cutoff)
outlier_table <- data.frame(Observation = outliers, Mahalanobis_Distance = cdud[outliers])
library(DT)
datatable(outlier_table,
options = list(pageLength = 10, autoWidth = TRUE),
caption = "Çok Değişkenli Aykırı Gözlemler")Analizler sonucunda, 774 adet çok boyutlu uç değer tespit edilmiştir. Bu değerlerin silinmesi analizlerin güvenilirliğini artırabilir.
Çok Değişkenli Uç Değerlerin Silinmesi
pisa_turkiye_listwise_clean <- pisa_turkiye_listwise[-outliers, ]
print(paste("Silinen toplam gözlem sayısı:", length(outliers)))## [1] "Silinen toplam gözlem sayısı: 770"Kaç gözlem kaldı?
## [1] "Kalan toplam gözlem sayısı: 6421"Çok değişkenli uç değerlerin silinmesi sonrası 7191 gözlemden 6417 gözlem kaldığı görülmektedir. Yani, 774 gözlem silinmiştir.
Normallik
Tek değişkenli ve çok değişkenli normallğe bakılması gerekmektedir. Tek değişkenli normallik için skewness ve kurtosis değerlerine, çok değişkenli normallik için Henze-Zirkler çok değişkenli normallik testi sonucuna bakacağım.
Tek Değişkenli Normallik
Önce skewness ve kurtosis değerlerine bakalım.
pisa_turkiye_listwise_clean %>%
dplyr::select(BULLIED, BELONG, FEELSAFE) %>%
psych::describe() %>%
kable() %>%
kable_styling(full_width = TRUE, position = "center")| vars | n | mean | sd | median | trimmed | mad | min | max | range | skew | kurtosis | se | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| BULLIED | 1 | 6421 | -0.2090119 | 1.0202158 | -0.2670 | -0.3145767 | 1.4247786 | -1.2280 | 3.4060 | 4.6340 | 0.5025799 | -0.8854824 | 0.0127318 |
| BELONG | 2 | 6421 | -0.3623587 | 0.7178521 | -0.3485 | -0.3773060 | 0.7295875 | -2.8882 | 2.2967 | 5.1849 | 0.1493175 | 0.1539826 | 0.0089585 |
| FEELSAFE | 3 | 6421 | -0.5248476 | 0.9415930 | -0.7560 | -0.5512762 | 0.8564980 | -2.7886 | 1.1246 | 3.9132 | 0.3619832 | -0.5316392 | 0.0117506 |
skewness ve kurtosis değerlerine bakıldığında, tüm değişkenlerin +1 -1 aralığında değer aldığı ve veri setinin tek değişkenli normal dağılıma sahip olduğunun kabul edilebileceği anlaşılmıştır.
Çok Değişkenli Normallik
Henze-Zirkler ve Q-Q plot ile, çok değişkenli normallik testi yaparak bu durumu değerlendirebiliriz.
henze_zirkler <- MVN::mvn(data = pisa_turkiye_listwise_clean, mvnTest = "hz", multivariatePlot = "qq")
## $multivariateNormality
## Test HZ p value MVN
## 1 Henze-Zirkler 75.38457 0 NO
##
## $univariateNormality
## Test Variable Statistic p value Normality
## 1 Anderson-Darling BULLIED 338.6127 <0.001 NO
## 2 Anderson-Darling BELONG 16.4306 <0.001 NO
## 3 Anderson-Darling FEELSAFE 249.4568 <0.001 NO
##
## $Descriptives
## n Mean Std.Dev Median Min Max 25th 75th
## BULLIED 6421 -0.2090119 1.0202158 -0.2670 -1.2280 3.4060 -1.2280 0.5869
## BELONG 6421 -0.3623587 0.7178521 -0.3485 -2.8882 2.2967 -0.8709 0.1033
## FEELSAFE 6421 -0.5248476 0.9415930 -0.7560 -2.7886 1.1246 -1.1042 0.1368
## Skew Kurtosis
## BULLIED 0.5025799 -0.8854824
## BELONG 0.1493175 0.1539826
## FEELSAFE 0.3619832 -0.5316392Q-Q plottaki çizgiden sapmalar ve Henze-Zirkler testi sonucuna göre, veri setinin çok değişkenli normal dağılım göstermediği söylenebilir.
Henze-Zirkler testinde p < 0.05 olduğu için çok değişkenli normallik reddedilmiştir. Bu testte p değerinin .05’ten büyük çıkması beklenir. Bu durumda, ya parametrik olmayan testler ya da bootstrap yöntemi kullanılmalıdır. (Bootstrapping yöntemini kullanacağım)
Değişkenlerin histogramlarını çizerek dağılımlarını görselleştirelim.
Zorbalığa maruz kalma;
pisa_turkiye_listwise %>% ggplot(aes(x = BULLIED)) +
geom_histogram(fill = "lightblue", color = "black", bins = 25) +
geom_vline(aes(xintercept = mean(BULLIED), color = "mean"), linetype = 'solid', size = 1) +
geom_vline(aes(xintercept = median(BULLIED), color = "median"), linetype = "solid", size = 1) +
scale_color_manual(values = c("mean" = "red", "median" = "blue")) +
labs(title = "Zorbalığa Maruz Kalma Değişkeninin Dağılımı", x = "Zorbalığa Maruz Kalma", y = "Frekans") +
theme_minimal()## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
Bu değişkenin normal dağılıma uymadığı gafikten de görülmektedir.
Okulda aidiyet duygusu;
pisa_turkiye_listwise %>% ggplot(aes(x = BELONG)) +
geom_histogram(fill = "lightblue", color = "black", bins = 30) +
geom_vline(aes(xintercept = mean(BELONG), color = "mean"), linetype = "solid", size = 1) +
geom_vline(aes(xintercept = median(BELONG), color = "median"), linetype = "solid", size = 1) +
scale_color_manual(values = c("mean" = "red", "median" = "blue")) +
labs(title = "Okulda Aidiyet Duygusu Değişkeninin Dağılımı", x = "Okulda Aidiyet Duygusu", y = "Frekans") +
theme_minimal()
Güvende Hissetme;
pisa_turkiye_listwise %>% ggplot(aes(x = FEELSAFE)) +
geom_histogram(fill = "lightblue", color = "black", bins = 30) +
geom_vline(aes(xintercept = mean(FEELSAFE), color = "mean"), linetype = "solid", size = 1) +
geom_vline(aes(xintercept = median(FEELSAFE), color = "median"), linetype = "solid", size = 1) +
scale_color_manual(values = c("mean" = "red", "median" = "blue")) +
labs(title = "Güvende Hissetme Değişkeninin Dağılımı", x = "Güvende Hissetme", y = "Frekans") +
theme_minimal()
Doğrusallık
İki değişken arasında doğrusal bir iliskinin olup olmadığına bakmak için iki değişkenli saçılım grafiği incelenerek değerlendirilebilir. Eğer iki değişken de normal dağılıyorsa ve doğrusal olarak ilişkiliyse, saçılım grafiği oval şeklindedir. 3 değişkenin de ikili grafiklerini aşağıdaki çıktıda görebiliriz.
pisa_turkiye_listwise %>%
ggplot(aes(x = BULLIED, y = BELONG)) +
geom_point(color='red', alpha=0.2) +
labs(title = "Zorbalığa Maruz Kalma ve Okulda Aidiyet Duygusu Arasındaki İlişki", x = "Zorbalığa Maruz Kalma", y = "Okulda Aidiyet Duygusu") +
theme_minimal()
pisa_turkiye_listwise %>%
ggplot(aes(x = FEELSAFE, y = BELONG)) +
geom_point(color='darkblue', alpha=0.2) +
labs(title = "Okulda Aidiyet Duygusu ve Güvende Hissetme Arasındaki İlişki", x = "Okulda Aidiyet Duygusu", y = "Güvende Hissetme") +
theme_minimal()
Grafiklerin oval olduğu (ben ovale de benzetmiş olabilirim hocam:) görülmektedir.
Bir de pairs fonksiyonu ile değişkenler arasındaki
ilişkiyi inceleyelim.
Doğrusallık varsayımının karşılandığı söylenebilir.
Varyansların homojenliği
Homoskedastisite, regresyon modelinde artıkların (residuals) sabit bir varyansa sahip olması anlamına gelir. Yani, bağımsız değişken arttıkça artıkların (hataların) saçılımının değişmemesi gerekir. Bu varsayım ihlal edilirse heteroskedastisite ortaya çıkar ve modelin güvenilirliği azalır. Bu varsayımı grafiksel ve istatistiksel testlerle kontrol edebiliriz (Umarım doğru yorumlamışımdır hocam.)
model_vhtest <- lm(BELONG ~ BULLIED + FEELSAFE, data = pisa_turkiye_listwise_clean)
plot(model_vhtest$fitted.values, residuals(model_vhtest),
xlab = "Tahmin Edilmiş Değerler",
ylab = "Artıklar",
main = "Artıklar - Tahmin Edilmiş Değerler")
abline(h = 0, col = "red")
Bu grafik yorumlanırken eğer noktalar rastgele ve simetrik dağılmışsa, homoskedastisite sağlanmıştır yorumu yapılmaktadır. Grafikte görüldüğü üzere, artıkların rastgele ve simetrik bir şekilde dağıldığı görülmektedir. Bu durumda homoskedastisite varsayımı sağlanmıştır.
Bir de Breusch-Pagan testi kullanılmaktadır. Bu test sonucunda p değeri 0.05’ten küçük olursa modelde değişen varyans (heteroskedastisite) problemi olduğu söylenebilir, Bu da regresyon katsayılarının güvenilirliğini etkileyebilir.
Levene testi ile varyans homojenliğini test edelim.
library(car)
pisa_turkiye_listwise_clean$BELONG_GROUP <- ifelse(pisa_turkiye_listwise_clean$BELONG < median(pisa_turkiye_listwise_clean$BELONG),"Düşük", "Yüksek")
leveneTest(BULLIED ~ BELONG_GROUP, data = pisa_turkiye_listwise_clean)## Warning in leveneTest.default(y = y, group = group, ...): group coerced to
## factor.## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 1 104.17 < 2.2e-16 ***
## 6419
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Levene testi sonucu, p-değeri 0.05’ten oldukça küçük. Bu durumda, gruplar arasında varyansların eşit olduğu (homojen olduğu) varsayımı reddedilir.
Multicollinearity
Değişkenler arasındaki ilişkiyi incelemek için korelasyon analizi yapalım. Çoklu bağlantı sorunu olup olmadığını kontrol etmek için VIF değerlerine bakalım. VIF değeri 10’dan büyük olan değişkenlerin çoklu bağlantı sorunu olduğu kabul edilir. Tolerans değeri 0.1’den küçük olan değişkenlerin de çoklu bağlantı sorunu olduğu kabul edilir.
pisa_numeric <- pisa_turkiye_listwise_clean %>%
mutate(across(where(is.factor), as.numeric)) %>%
select_if(is.numeric)
cor(pisa_numeric)## BULLIED BELONG FEELSAFE
## BULLIED 1.0000000 -0.2989147 -0.2469997
## BELONG -0.2989147 1.0000000 0.3895851
## FEELSAFE -0.2469997 0.3895851 1.0000000## Warning: package 'PerformanceAnalytics' was built under R version 4.4.3## Loading required package: xts## Loading required package: zoo## Warning: package 'zoo' was built under R version 4.4.3##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric##
## ######################### Warning from 'xts' package ##########################
## # #
## # The dplyr lag() function breaks how base R's lag() function is supposed to #
## # work, which breaks lag(my_xts). Calls to lag(my_xts) that you type or #
## # source() into this session won't work correctly. #
## # #
## # Use stats::lag() to make sure you're not using dplyr::lag(), or you can add #
## # conflictRules('dplyr', exclude = 'lag') to your .Rprofile to stop #
## # dplyr from breaking base R's lag() function. #
## # #
## # Code in packages is not affected. It's protected by R's namespace mechanism #
## # Set `options(xts.warn_dplyr_breaks_lag = FALSE)` to suppress this warning. #
## # #
## #################################################################################
## Attaching package: 'xts'## The following objects are masked from 'package:dplyr':
##
## first, last##
## Attaching package: 'PerformanceAnalytics'## The following objects are masked from 'package:moments':
##
## kurtosis, skewness## The following object is masked from 'package:graphics':
##
## legenddf <- pisa_turkiye_listwise_clean[, c("BULLIED", "BELONG", "FEELSAFE")]
chart.Correlation(df, histogram=TRUE, pch=19)
model_mediation <- lm(BELONG ~ BULLIED + FEELSAFE, data = pisa_turkiye_listwise_clean)
vif_table <- ols_vif_tol(model_mediation)
vif_table %>%
kable(digits = 2, caption = "VIF ve Tolerans Değerleri") %>%
kable_styling(full_width = TRUE, position = "center")| Variables | Tolerance | VIF |
|---|---|---|
| BULLIED | 0.94 | 1.06 |
| FEELSAFE | 0.94 | 1.06 |
Korelasyon matrisi incelendiğinde, değişkenler arasında .001 düzeyinde anlamlı pozitif ve negatif ilişkiler olduğu görülmektedir.
Korelasyonlari VIF ve Tolerans değerlerini birlikte incelediğimizde, değişkenler arasında çoklu bağlantı sorunu olmadığı görülmektedir.
ARACI ETKİ ANALİZİ DENEMESİ
Aracı etki analizi, bir bağımsız değişkenin bağımlı değişken üzerindeki etkisinin bir kısmının bir aracı değişken aracılığıyla gerçekleşip gerçekleşmediğini inceleyen bir analiz türüdür. Bu analizde, bağımsız değişkenin bağımlı değişken üzerindeki etkisi, aracı değişkenin bağımlı değişken üzerindeki etkisi ve bağımsız değişkenin aracı değişken üzerindeki etkisi hesaplanır. Bu analizde, bağımsız değişkenin bağımlı değişken üzerindeki etkisinin bir kısmının aracı değişken aracılığıyla gerçekleşip gerçekleşmediği test edilir.
Çok değişkenli normallik varsayımı sağlanamadığı için analizlerde bootstrap yöntemini kullanacağım.
Lavaan Paketi ile Aracı Etki Analizi Denemesi
Bootstrap yöntemi ile yapılan analizde 5000 yeniden örnekleme ve %95 güven aralığı kullanılacaktır.
## This is lavaan 0.6-19
## lavaan is FREE software! Please report any bugs.##
## Attaching package: 'lavaan'## The following object is masked from 'package:psych':
##
## cor2covmed_model <- '
FEELSAFE ~ a*BULLIED
BELONG ~ b*FEELSAFE + c*BULLIED
indirect := a*b
total := c + (a*b)
'
fit <- sem(med_model, data=pisa_turkiye_listwise_clean, se="bootstrap", bootstrap=5000)## Warning: lavaan->lav_model_nvcov_bootstrap():
## 7 bootstrap runs failed or did not converge.## lavaan 0.6-19 ended normally after 1 iteration
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 5
##
## Number of observations 6421
##
## Model Test User Model:
##
## Test statistic 0.000
## Degrees of freedom 0
##
## Model Test Baseline Model:
##
## Test statistic 1801.208
## Degrees of freedom 3
## P-value 0.000
##
## User Model versus Baseline Model:
##
## Comparative Fit Index (CFI) 1.000
## Tucker-Lewis Index (TLI) 1.000
##
## Loglikelihood and Information Criteria:
##
## Loglikelihood user model (H0) -14805.466
## Loglikelihood unrestricted model (H1) NA
##
## Akaike (AIC) 29620.933
## Bayesian (BIC) 29654.769
## Sample-size adjusted Bayesian (SABIC) 29638.881
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.000
## 90 Percent confidence interval - lower 0.000
## 90 Percent confidence interval - upper 0.000
## P-value H_0: RMSEA <= 0.050 NA
## P-value H_0: RMSEA >= 0.080 NA
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.000
##
## Parameter Estimates:
##
## Standard errors Bootstrap
## Number of requested bootstrap draws 5000
## Number of successful bootstrap draws 4993
##
## Regressions:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## FEELSAFE ~
## BULLIED (a) -0.228 0.011 -21.502 0.000 -0.228 -0.247
## BELONG ~
## FEELSAFE (b) 0.256 0.009 29.288 0.000 0.256 0.336
## BULLIED (c) -0.152 0.008 -19.037 0.000 -0.152 -0.216
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .FEELSAFE 0.832 0.012 66.836 0.000 0.832 0.939
## .BELONG 0.414 0.008 52.741 0.000 0.414 0.804
##
## Defined Parameters:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## indirect -0.058 0.003 -17.231 0.000 -0.058 -0.083
## total -0.210 0.008 -25.855 0.000 -0.210 -0.299Sonuçları tablo ve grafiklerle görelim
library(lavaan)
library(dplyr)
library(knitr)
library(kableExtra)
results <- parameterEstimates(fit, standardized = TRUE, level = 0.95) %>%
as.data.frame()
results_table <- results %>%
dplyr::select(lhs, op, rhs, label, est, se, z, pvalue, ci.lower, ci.upper, std.all)
kable(results_table, format = "markdown", caption = "Lavaan Model Sonuçları") %>%
kable_styling(full_width = FALSE)| lhs | op | rhs | label | est | se | z | pvalue | ci.lower | ci.upper | std.all |
|---|---|---|---|---|---|---|---|---|---|---|
| FEELSAFE | ~ | BULLIED | a | -0.2279647 | 0.0106021 | -21.50180 | 0 | -0.2484866 | -0.2074794 | -0.2469997 |
| BELONG | ~ | FEELSAFE | b | 0.2563647 | 0.0087532 | 29.28808 | 0 | 0.2389392 | 0.2738236 | 0.3362687 |
| BELONG | ~ | BULLIED | c | -0.1518825 | 0.0079784 | -19.03675 | 0 | -0.1680706 | -0.1364998 | -0.2158564 |
| FEELSAFE | ~~ | FEELSAFE | 0.8323773 | 0.0124541 | 66.83553 | 0 | 0.8082504 | 0.8572185 | 0.9389911 | |
| BELONG | ~~ | BELONG | 0.4144893 | 0.0078590 | 52.74089 | 0 | 0.3988805 | 0.4298745 | 0.8044721 | |
| BULLIED | ~~ | BULLIED | 1.0406782 | 0.0000000 | NA | NA | 1.0406782 | 1.0406782 | 1.0000000 | |
| indirect | := | a*b | indirect | -0.0584421 | 0.0033917 | -17.23093 | 0 | -0.0653339 | -0.0519228 | -0.0830583 |
| total | := | c+(a*b) | total | -0.2103246 | 0.0081349 | -25.85455 | 0 | -0.2265105 | -0.1942606 | -0.2989147 |
Grafik
effects_df <- results %>%
filter(label %in% c("a", "indirect", "total")) %>%
mutate(Effect = case_when(
label == "a" ~ "Doğrudan Etki (Direct)",
label == "indirect" ~ "Dolaylı Etki (Indirect)",
label == "total" ~ "Toplam Etki (Total)"
)) %>%
rename(Estimate = est, Lower_CI = ci.lower, Upper_CI = ci.upper)
ggplot(effects_df, aes(x = Effect, y = Estimate, fill = Effect)) +
geom_bar(stat = "identity", color = "black") +
geom_errorbar(aes(ymin = Lower_CI, ymax = Upper_CI), width = 0.2) +
theme_minimal() +
labs(title = "Aracılık Analizi Sonuçları", x = "Etki Türü", y = "Tahmini Etki") +
theme(axis.text.x = element_text(angle = 15, hjust = 1)) +
guides(fill = "none")
Sonuçlar incelendiğinde, FEELSAFE değişkeninin, BULLIED’in BELONG üzerindeki etkisinin bir kısmına aracılık ettiği görülmektedir.
Doğrudan Etki -0.216 olup, p değeri 0.05’ten küçüktür. Yani Zorbalığa uğramanın okul aidiyetine olan negatif doğrudan etkisi anlamlıdır. Toplam Etki -0.299 olarak hesaplanmıştır ve p değeri 0.05’ten küçüktür, yani toplam etki de anlamlıdır. Anlamlılık güven aralıkları tarafından da desteklenmektedir. Aracı etki ise -0.083 olarak hesaplanmıştır ve p değeri 0.05’ten küçüktür. Yani, FEELSAFE değişkeninin BULLIED’in BELONG üzerindeki etkisinin bir kısmına aracılık ettiği görülmektedir.