Veri setinin yüklenmesi

data(iris)
head(iris)
##   Sepal.Length Sepal.Width Petal.Length Petal.Width Species
## 1          5.1         3.5          1.4         0.2  setosa
## 2          4.9         3.0          1.4         0.2  setosa
## 3          4.7         3.2          1.3         0.2  setosa
## 4          4.6         3.1          1.5         0.2  setosa
## 5          5.0         3.6          1.4         0.2  setosa
## 6          5.4         3.9          1.7         0.4  setosa

Dummy Değişken Oluşturma

library(fastDummies)   #seçilen sütunda dummy değişken oluşturdum.

dat <- dummy_cols(iris, select_columns = 'Species')
table(dat$Species_setosa)
## 
##   0   1 
## 100  50
table(dat$Species_versicolor)  #frekansları inceledim
## 
##   0   1 
## 100  50
table(dat$Species_virginica)
## 
##   0   1 
## 100  50

Model Kurulması

#Çanak yaprak uzunluğu bağımlı değişken ve iki farklı tür bağımsız değişken o.ü modeli kurdum

model_species <- lm(Sepal.Length ~ Species_setosa + Species_versicolor , 
                  data=dat)
model_species       
## 
## Call:
## lm(formula = Sepal.Length ~ Species_setosa + Species_versicolor, 
##     data = dat)
## 
## Coefficients:
##        (Intercept)      Species_setosa  Species_versicolor  
##              6.588              -1.582              -0.652

#Virginica türündeki bir çiçeğin ortalama çanak yaprak uzunluğu ≈ 6.588 cm’dir.

#Setosa türündeki bir çiçeğin ortalama sepal length’i ≈ 6.588 − 1.582 = 5.006 cm

#Versicolor türü bir çiçeğin ortalama sepal length’i ≈ 6.588 − 0.652 = 5.936 cm

library(broom)
tidy(model_species)   #benzer değerleri elde ettim
## # A tibble: 3 × 5
##   term               estimate std.error statistic   p.value
##   <chr>                 <dbl>     <dbl>     <dbl>     <dbl>
## 1 (Intercept)           6.59     0.0728     90.5  8.51e-131
## 2 Species_setosa       -1.58     0.103     -15.4  2.21e- 32
## 3 Species_versicolor   -0.652    0.103      -6.33 2.77e-  9
model_1 <- lm(Sepal.Length ~ Species , 
                  data=dat)
library(broom)
tidy(model_1)    #farklı sıralarda verdi ama olsun farketmez ikiside aynı 
## # A tibble: 3 × 5
##   term              estimate std.error statistic   p.value
##   <chr>                <dbl>     <dbl>     <dbl>     <dbl>
## 1 (Intercept)           5.01    0.0728     68.8  1.13e-113
## 2 Speciesversicolor     0.93    0.103       9.03 8.77e- 16
## 3 Speciesvirginica      1.58    0.103      15.4  2.21e- 32

İki yönlü etkileşim etkisi

library(MASS)
data(Boston)
dat1 <- Boston[, c( "medv", "rm", "lstat")]
head (dat1)
##   medv    rm lstat
## 1 24.0 6.575  4.98
## 2 21.6 6.421  9.14
## 3 34.7 7.185  4.03
## 4 33.4 6.998  2.94
## 5 36.2 7.147  5.33
## 6 28.7 6.430  5.21
library(GGally)
## Zorunlu paket yükleniyor: ggplot2
## Registered S3 method overwritten by 'GGally':
##   method from   
##   +.gg   ggplot2
ggpairs(dat1)

# oda sayısı (rm) ve düşük sosyoekonomik durum (lstat) ile medyan konut fiyatı (medv) arasındaki ilişkiler şekilde gösterilmiştir. Aslında her iki bağımsız değişkende medyan konut fiyatını etkileyen birer bağımsız değişkendir. Ancak ben korelasyon değeri düşük ve anlamsız olmayacak şekilde veri seti bulamadığım için bu iki değişken üzerinden analizlere devam edeceğim.

#oda sayısı ve düşük sosyoekonomik durum etkileşimini test etmek istiyorum. Bu iki değişkenin çarpımıyla yeni bir değişken oluşturacağım çoklu bağlantı problemi için de bu değişkenleri merkezleyeceğim (ortalamayı sütundaki her bir değerden çıkaracağım). Bu yüzden 0 değerleri ortalama düzeylerini temsil edecekler.

attach(dat1)
merkez_rm   <- c(scale(rm, center=TRUE, scale=FALSE)) 
merkez_lstat    <- c(scale(lstat, center=TRUE, scale=FALSE))
dat2 <- data.frame(cbind(dat1, cross = rm*lstat, 
                                merkez_rm, merkez_lstat, cross_merkez = merkez_rm*merkez_lstat))
cormat <- cor(dat2[, c(1: 4)])
# Correlation matrix plot:
ggcorrplot::ggcorrplot(corr = cormat, # correlation matrix
                       type = "lower", # print only the lower part of the correlation matrix
                       hc.order = TRUE, # hierarchical clustering
                       show.diag = TRUE, # show the diagonal values of 1
                       lab = TRUE, # add correlation values as labels
                       lab_size = 3) # Size of the labels

# merkezlenmemiş verilerde lstat ile cross değişkeni arasıdaki korelasyon çok yüksek çıktı rm de o kadar yüksek değil açıkçası.

cormat <- cor(dat2[, c(1,5: 7)])
# Correlation matrix plot:
ggcorrplot::ggcorrplot(corr = cormat, # correlation matrix
                       type = "lower", # print only the lower part of the correlation matrix
                       hc.order = TRUE, # hierarchical clustering
                       show.diag = TRUE, # show the diagonal values of 1
                       lab = TRUE, # add correlation values as labels
                       lab_size = 3) # Size of the labels

merkezlenmiş değerler ile yapılan analizlerde korelasyon değerleri oldukça düşük çıktı. Son çalışmalar VİF olsa da analizler yapılabilir diyor.

Etkileşim etkisi modelleri

normal_model <- lm(medv  ~ merkez_rm  + merkez_lstat, data=dat2  )
cross_model <- lm(medv  ~ merkez_rm  + merkez_lstat + merkez_rm*merkez_lstat, data=dat2)  
glance(normal_model)$r.squared
## [1] 0.6385616
glance(cross_model)$r.squared
## [1] 0.7402461

#rm ve lstat etkileşim etkisi açıklanan varyans oranı %10 oranında arttırdı. etkileşimli modeli tercih edebiliriz.

glance(normal_model)$adj.r.squared
## [1] 0.6371245
glance(cross_model)$adj.r.squared
## [1] 0.7386937

#benzer durum düzenlenmiş R2 değerleri için de geçerli.

Etkileşimli modelin Yorumlanması

tidy(cross_model)
## # A tibble: 4 × 5
##   term                   estimate std.error statistic   p.value
##   <chr>                     <dbl>     <dbl>     <dbl>     <dbl>
## 1 (Intercept)              21.0      0.234      89.7  2.90e-311
## 2 merkez_rm                 3.57     0.393       9.08 2.48e- 18
## 3 merkez_lstat             -0.854    0.0401    -21.3  2.92e- 72
## 4 merkez_rm:merkez_lstat   -0.485    0.0346    -14.0  6.50e- 38

# rm, lstat ve cross değişkenlerinin katsayıları sırasıyla 3,56, -0,85 ve -0,48 olarak bulunmuştur. Üç değişken de istatistiksel olarak anlamlıdır.

Etkileşim etkisinin grafikle incelenmesi

min(dat2$merkez_rm, na.rm = TRUE)
## [1] -2.723634
# Maksimum değer
max(dat2$merkez_rm, na.rm = TRUE)
## [1] 2.495366

#-1,5, 0 ve 1, değerleri düşük orta ve yüksek olarak belirlenebilir.

library(dplyr)
## 
## Attaching package: 'dplyr'
## The following object is masked from 'package:MASS':
## 
##     select
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
dat2 <- dat2 %>% 
  mutate(
    merkez_rm_kategor = case_when(
      merkez_rm <= -1.5 ~ "dusuk",          # -1.5 ve altı
      merkez_rm > -1.5 & merkez_rm < 1.5 ~ "orta",  # -1.5 ile +1.5 arası
      merkez_rm >= 1.5 ~ "yuksek"           # 1.5 ve üzeri
    )
  ) %>%
  arrange(merkez_lstat)
library(ggplot2)

ggplot(dat2, aes(x = merkez_lstat, y = medv, color = merkez_rm_kategor)) +
  geom_smooth(method = "lm", se = TRUE) +  
  labs(
    x = "LSTAT",
    y = "MEDV",
    title = "LSTAT ve MEDV İlişkisi (RM Kategorilerine Göre)",
    color = "RM Kategori"
  ) +
  theme_minimal()
## `geom_smooth()` using formula = 'y ~ x'

  • rm (oda sayısı) değişkeninin farklı kategorilerinde ev fiyatları değişkeninde bariz bir değişiklik gözlenmemiştir. Orta düzeyde ki rm değişkenlerinde rekabetçi bir etki gözlenmiştir.

Etkileşim etkisinin farklı paketlerle incelenmesi

library(stargazer)
## 
## Please cite as:
##  Hlavac, Marek (2022). stargazer: Well-Formatted Regression and Summary Statistics Tables.
##  R package version 5.2.3. https://CRAN.R-project.org/package=stargazer
# Model özetlerini stargazer ile yazdırma
stargazer(normal_model, cross_model, type = "text", title = "Duzencisi Etkisi")
## 
## Duzencisi Etkisi
## ========================================================================
##                                       Dependent variable:               
##                        -------------------------------------------------
##                                              medv                       
##                                  (1)                      (2)           
## ------------------------------------------------------------------------
## merkez_rm                      5.095***                 3.565***        
##                                (0.444)                  (0.393)         
##                                                                         
## merkez_lstat                  -0.642***                -0.854***        
##                                (0.044)                  (0.040)         
##                                                                         
## merkez_rm:merkez_lstat                                 -0.485***        
##                                                         (0.035)         
##                                                                         
## Constant                      22.533***                21.042***        
##                                (0.246)                  (0.234)         
##                                                                         
## ------------------------------------------------------------------------
## Observations                     506                      506           
## R2                              0.639                    0.740          
## Adjusted R2                     0.637                    0.739          
## Residual Std. Error        5.540 (df = 503)         4.701 (df = 502)    
## F Statistic            444.331*** (df = 2; 503) 476.866*** (df = 3; 502)
## ========================================================================
## Note:                                        *p<0.1; **p<0.05; ***p<0.01

#Etki büyüklüğünün incelenmesi de oldukça önemlidir. Düzenleyicilik etki analizini hesaplamada genellikle Cohen f^2 kullanılır ve küçük, orta, büyük olarak değerlendirmede .005 , .01, .025 kesme değerleri kullanılır (Aguinis vd., 2015; Kenny, 2018)

library(effectsize)
cohens_f_squared(cross_model)
## # Effect Size for ANOVA (Type I)
## 
## Parameter              | Cohen's f2 (partial) |      95% CI
## -----------------------------------------------------------
## merkez_rm              |                 1.86 | [1.59, Inf]
## merkez_lstat           |                 0.60 | [0.47, Inf]
## merkez_rm:merkez_lstat |                 0.39 | [0.30, Inf]
## 
## - One-sided CIs: upper bound fixed at [Inf].

Cohen f^2 değeri eklenen değişken için oldukça yeterli bir etkileşim etkisine sahip ve elde edilen katsayı değerimizde orta düzeyde.

library(rockchalk)
## 
## Attaching package: 'rockchalk'
## The following object is masked from 'package:effectsize':
## 
##     standardize
## The following object is masked from 'package:dplyr':
## 
##     summarize
## The following object is masked from 'package:MASS':
## 
##     mvrnorm
ps  <- plotSlopes(cross_model, plotx="merkez_lstat", modx="merkez_rm", xlab = "ekonomik duzey", ylab = "Ev fıyatları", modxVals = c("std.dev"))

# oda sayısı artıkça ev fiyatlarında artış görülmektedir. Sosyo ekonomik düzey azaldıkça ev fiyatlarında bir azalma mevcut. Ancak bu istenmeyen bir durum böyle durumlarda o bölgenin ileride yatırım potansiyeli yüksek mahallelerolabileceği söylenebilir.

library(car)
## Zorunlu paket yükleniyor: carData
## 
## Attaching package: 'car'
## The following object is masked from 'package:dplyr':
## 
##     recode
avPlots(cross_model)

###oda sayısı ve sosyo ekonomik düzeyin etkisi oldukça belirgin.

#Anscombe’un Dörtlüsü göz önünde bulundurularak mutlaka veri seti incelenmeli ve grafikler üzerinden de kontroller yapılmamlıdır.