Örnek3.1 Kolej Gpa’sının Belirleyicileri

library(wooldridge)
data(gpa1)
library(rmarkdown)
paged_table(gpa1)

Örnek için kullanacağımız değişkenlerin tanımlarını şöyle yapabiliriz.

colGPA: MSU GPA, Michigan State Üniversitesi not ortalaması

hsGPA: high school GPA, lise not ortalaması

ACT: ‘achievement’ score, ACT, bir öğrencinin okulda ne öğrendiğini ölçen bir başarı testidir. SAT daha çok bir yetenek testidir, muhakeme ve sözel yetenekleri test eder. Üniversiteler SAT veya ACT sonuçları ile öğrenci kabul eder.

coklureg <- lm(colGPA~ hsGPA+ACT,data= gpa1)
summary(coklureg)
## 
## Call:
## lm(formula = colGPA ~ hsGPA + ACT, data = gpa1)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.85442 -0.24666 -0.02614  0.28127  0.85357 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 1.286328   0.340822   3.774 0.000238 ***
## hsGPA       0.453456   0.095813   4.733 5.42e-06 ***
## ACT         0.009426   0.010777   0.875 0.383297    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3403 on 138 degrees of freedom
## Multiple R-squared:  0.1764, Adjusted R-squared:  0.1645 
## F-statistic: 14.78 on 2 and 138 DF,  p-value: 1.526e-06

colGPA=1,29+0,453hsGPA+0,0094ACT+u

Peki bu modeli nasıl yorumlayacağız? İlk olarak kesişim katsayısı, hsGPA ve ACT’nin her ikisi de sıfır olursa, öngörülen üniversite GPA’sini göstermektedir ve örneğimizde 1.29 olarak hesaplanmıştır. Üniversiteye giden hiç kimsenin lise not ortalaması sıfır veya başarı testi sıfır olmadığı için, bu denklemdeki kesişim kendi başına anlamlı değildir. Eğim katsayılarına bakınca ikisinin de bekleneceği gibi pozitif olduğunu görüyoruz. Eğer ACT puanını sabit tutarsak, hsGPA’deki 1 birimlik artış colGPA üzerinde 0.453 puanlık bir artış yaratacaktır. Başka bir deyişle, eğer ACT puanları aynı olan iki öğrenci varsa ve bir öğrencinin lise not ortalaması diğer öğrenciden bir puan yüksekse, bu öğrencinin üniversite not ortalamasının yarım puan daha fazla olmasını tahmin edebiliriz.

Üniversite not ortalamalarının 4 üzerinden, ACT sonuçlarının 36 puan üzerinden hesaplandığı düşünülürse, ACT’nin eğimi 0.0094’ün üniversite not ortalaması üzerinde fazla bir etkisi olmadığı düşünülebilir. Hata paylarını göz önüne aldığımızda ACT’nin sadece küçük değil aynı zamanda istatiksel olarak da anlamlı olmadığını göreceğiz.

Basit regresyon örneğimizi hatırlarsak

basitreg <- lm(colGPA~ACT,data = gpa1)
summary(basitreg)
## 
## Call:
## lm(formula = colGPA ~ ACT, data = gpa1)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.85251 -0.25251 -0.04426  0.26400  0.89336 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2.40298    0.26420   9.095  8.8e-16 ***
## ACT          0.02706    0.01086   2.491   0.0139 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3656 on 139 degrees of freedom
## Multiple R-squared:  0.04275,    Adjusted R-squared:  0.03586 
## F-statistic: 6.207 on 1 and 139 DF,  p-value: 0.0139

ACT’nin katsayısı çoklu regresyona göre bir hayli büyük kalıyor (0.02706/0.009426=2.87 kat daha fazla). ACT değerinin çoklu regresyonda lise not ortalaması eklenince azaldığını gözlemliyoruz. Çünkü lise not ortalaması ve ACT arasında yüksek bir korelasyon var. Bu durum bize aslında basit regresyon kullanıldığında temel varsayımı ihlal ettiğimize dair güçlü bir ipucu veriyor.

Bu örnekte son olarak basit ve iki bağımsız değişkenli regresyon grafiklerini gösterelim. Basit regresyon örneğimiz için, dağılım grafiğine bir regresyon çizgisi eklememizin yeterli olduğunu hatırlayalım.

plot(gpa1$ACT,gpa1$colGPA)
abline(basitreg)

Bu grafik sadece iki boyutlu olduğundan, iki boyutlu bilgisayar ekranlarımızda göstermek daha kolaydır. Y eksenine üniversite not ortalaması, X eksenine ACT puanlarını koyuyoruz. Son olarak bulduğumuz eğim ve kesen katsayılarını grafiğimizde temsil eden doğrumuzu çiziyoruz. Ancak çok değişkenli regresyon grafiklerini, iki boyutlu bilgisayar ekranlarda göstermek biraz daha zor olabilir. En azından örneğimizdeki gibi üç boyutlu bir grafiğe ihtiyacımız var. Önce bütün değişkenler için üç boyutlu bir dağılım grafiği çizelim daha sonra doğru yerine bir düzlem eklememiz gerekecek.

library(plotly)
## Loading required package: ggplot2
## 
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
## 
##     last_plot
## The following object is masked from 'package:stats':
## 
##     filter
## The following object is masked from 'package:graphics':
## 
##     layout
library(reshape2)

cf.mod <- coef(coklureg)
hsGPA <- seq(min(gpa1$hsGPA),max(gpa1$hsGPA),length.out=50)
ACT <- seq(min(gpa1$ACT),max(gpa1$ACT),length.out=50)
colGPA <- t(outer(hsGPA, ACT, function(x,y) cf.mod[1]+cf.mod[2]*x+cf.mod[3]*y))
p <- plot_ly(x=~hsGPA, y=~ACT, z=~colGPA,type="surface") %>%
  add_trace(data=gpa1, x=gpa1$hsGPA, y=gpa1$ACT, z=gpa1$colGPA, mode="markers", type="scatter3d", opacity=0.5)
p

Üç boyutlu grafiğimiz olduğu için artık üzerine gidip sağa sola çevirebilirsiniz. Tahmin edebileceğiniz gibi üç boyutlu grafikleri algılamamız daha büyük boyutlu grafiklerden daha kolay. Şimdi bu noktalarda oluşacak hata payını minimize edecek bir düzleme ihtiyacımız var. Bunun için çoklu regresyonumuzdan yardım alacağız.

Sizin şimdilik bu grafikleri çizerken kullandığımız kodları anlamanıza gerek yok. Sadece grafiklerin şekline odaklanın. Bu grafik bize yapmış olduğumuz minimizasyon işleminin neden basit regresyona bir değişken daha eklediğimizde farklılaştığı ve ek değişkenin toplamda neden daha çok varyasyon açıklamaya yardım ettiği hakkında daha fazla bir içgüdü sunabilir.

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
stargazer(list(coklureg,basitreg),type = "text")
## 
## =================================================================
##                                  Dependent variable:             
##                     ---------------------------------------------
##                                        colGPA                    
##                               (1)                    (2)         
## -----------------------------------------------------------------
## hsGPA                      0.453***                              
##                             (0.096)                              
##                                                                  
## ACT                          0.009                 0.027**       
##                             (0.011)                (0.011)       
##                                                                  
## Constant                   1.286***               2.403***       
##                             (0.341)                (0.264)       
##                                                                  
## -----------------------------------------------------------------
## Observations                  141                    141         
## R2                           0.176                  0.043        
## Adjusted R2                  0.164                  0.036        
## Residual Std. Error    0.340 (df = 138)       0.366 (df = 139)   
## F Statistic         14.781*** (df = 2; 138) 6.207** (df = 1; 139)
## =================================================================
## Note:                                 *p<0.1; **p<0.05; ***p<0.01

iki bağımsız değişkenin birbirleri arasındaki basit regresyon ilişkisinden de bulunabilir

ACT=S0+S1hsGPA

sapmareg<-lm(hsGPA~ACT,data = gpa1)
stargazer(list(sapmareg),type = "text")
## 
## ===============================================
##                         Dependent variable:    
##                     ---------------------------
##                                hsGPA           
## -----------------------------------------------
## ACT                          0.039***          
##                               (0.009)          
##                                                
## Constant                     2.463***          
##                               (0.218)          
##                                                
## -----------------------------------------------
## Observations                    141            
## R2                             0.120           
## Adjusted R2                    0.113           
## Residual Std. Error      0.301 (df = 139)      
## F Statistic           18.879*** (df = 1; 139)  
## ===============================================
## Note:               *p<0.1; **p<0.05; ***p<0.01

Regresyon sonucuyla bulduğumuz ve hesapladığımız S1 aynı olduğuna dikkat edin. Basit regresyon yaptığımızda gözlemleyemediğimiz değişken yüzünden olan sapma B2S1=0,453.0.0397=0.017 olur. Bu fark B1-B1 arasındaki farka eşittir.

Örneğimizi kullanarak bu parçaları çıkaralım.

standarthata<-summary(coklureg)$sigma
standarthata
## [1] 0.3403158

standart hatanın karesi bize regresyonumuzun hata payının varyansını verecektir.

sigmakare<-standarthata^2
sigmakare
## [1] 0.1158148

Regresörlerimiz arasındaki R2 hesaplamak için daha önce hesapladığımız sapmareg’i kullanabiliriz. Bu R2 bizim için lise not ortalamasıyla ACT sonuçları arasındaki ilişkiyi gösterecektir. Eğer R2 bir’e yakınsa çoklu eş doğrusallık muhtemel olabilir.

Rjkare<-summary(sapmareg)$r.squared
Rjkare
## [1] 0.1195815

İki bağımsız değişkenimizin basit regresyonundan elde ettiğimiz R2j yaklaşık olarak yüzde 12’dir.1/1-R2j varyans enflasyon faktörü (Variance inflation factor (VIF))’i hesaplayalım.

VIF<-1/(1-Rjkare)
VIF
## [1] 1.135823

son olarak n’i bulmalıyız.

n<-nobs(coklureg)
n
## [1] 141

son olarak var(xj) ’yi bulmamız lazım.

varxj<-var(gpa1$hsGPA)*(n-1)/n
varxj
## [1] 0.1016267

Şimdi formülümüzü kullanarak var(Bj)’yi hesaplayabiliriz.

var(Bj)=1/n.o”2/var(xj).1/1-R”2j

varbetaj<-(1/n)*(sigmakare/varxj)*VIF
varbetaj
## [1] 0.009180115

Bj’nin standart hatası var(Bj)’nin karaköküdür.

sdbetaj<-varbetaj^0.5
sdbetaj
## [1] 0.09581292

Tablodan hsGPA’in standart hatasına bakarsanız bulduğumuz rakamla aynı olduğunu göreceksiniz.

require(car)
## Loading required package: car
## Loading required package: carData
vif(coklureg)
##    hsGPA      ACT 
## 1.135823 1.135823

1.1358’in bizim daha önce hesapladığımız VIF değeriyle aynı olduğuna dikkat edin.