Varsayılım ki elimizde 3’er gözlemli iki değişken var. Değişkenler X ve Y olarak adlandırılmıştır.
Aşağıda kullanılan cbind () fonksiyonu, sütun bağlama işlevi, her iki veri çerçevesindeki satır sayısının eşit olması koşuluyla iki veri çerçevesini birleştirmek için kullanılmıştır.
X = cbind(c(1,2,3))Y = cbind(c(4,5,6))\[Σ X \cdot Y = 4 + 10 + 18 = 32 \]
olacaktır.
Bu sonucu matrislerle almanın kolay yolu \[ X^T \cdot Y\] bulmak olacaktır.
t(X) %*% Y## [,1]
## [1,] 32Bu sayıyı \(X^2\)’lerin toplamına bölmek istersek, solve fonksiyonunu kullanarak \[(X^T \cdot X)^{-1}\cdot X ^ T \cdot Y \] şeklinde yazabiliriz.
solve(t(X) %*% X)%*%t(X) %*% Y## [,1]
## [1,] 2.285714Yazılan son fonksiyon aslında regresyon katsayıları için olan formüldür.
\[ (X^T \cdot X)^ {-1} \cdot X^T \cdot Y \]
Fonksiyon Yazımı
Girdilerimiz X ve Y olacaktır. X çoklu regresyon yapılırken matris olarak sunulmalıdır.
Dogrusal_reg <- function(X, Y)
solve(t(X) %*% X)%*%t(X) %*% YŞimdi fonksiyonumuz hazır. İlk denememizi oluşturduğumuz X ve Y serileri için kullanalım.
Dogrusal_reg(X=X, Y=Y)## [,1]
## [1,] 2.285714Bu fonksiyonu birde herhangi bir veri seti üzerinden denersek.
bwght verisetini kullanarak cigs ve faminc kullanarak regresyon sonuçlarını tekrarlayalım.
Burada, cigs:sigara faminc: kıtlık anlamında kullanılmaktadır.
Bunun için aşağıdaki işlemleri uygulaya biliriz.
İlk olarak wooldridge paketini çalıştırıyoruz, ardından bwght veri setini data fonksiyonunu kullanarak çalıştırıyoruz.
library(wooldridge)
data("bwght")Fonksiyonumuzu kullanmak için, X verisetimizi oluşturmamız gerek, çünkü X bizim için sadece iki değişkenden oluşuyor. İsmini farklı yazabiliriz.
X yerine Xler ismi yazıla bilir.
head() fonksiyonu, veri setinin ilk altı satırını göstermketedir.
Xler = cbind(bwght$cigs, bwght$faminc)
head(Xler)## [,1] [,2]
## [1,] 0 13.5
## [2,] 0 7.5
## [3,] 0 0.5
## [4,] 0 15.5
## [5,] 0 27.5
## [6,] 0 7.5Ama X veri seti sadece X değişkenlerinden oluşmaz. Kesen diye adlandırılan bir katsayı daha var. Yani iki değil, üç katsayı bulunacak bu yüzden Xler veri setimize birlerden oluşan bir sütun daha eklememiz gerekmektedir.
Xler = cbind(bwght$cigs, bwght$faminc, c(rep(1, nrow(bwght))))
head(Xler)## [,1] [,2] [,3]
## [1,] 0 13.5 1
## [2,] 0 7.5 1
## [3,] 0 0.5 1
## [4,] 0 15.5 1
## [5,] 0 27.5 1
## [6,] 0 7.5 1Xler matrisi tamamlandı.
Y matrisini oluşturursak
Y = cbind(bwght$bwght)
head(Y)## [,1]
## [1,] 109
## [2,] 133
## [3,] 129
## [4,] 126
## [5,] 134
## [6,] 118Artık yeni regresyon kullanılabilir. Bunun için aşağıda yer alan işlemi uygulayabiliriz.
Dogrusal_reg(X=Xler, Y=Y)## [,1]
## [1,] -0.46340754
## [2,] 0.09276474
## [3,] 116.97413048Karşılaştırabileceğimiz gibi sonuçlar aynı.
lm() fonksiyonu, R’de doğrusal bir regresyon modeli oluşturur. Bu işlev, Y’nin sonuç değişkeni ve X’in öngörücü değişken olduğu bir Y ~ X R formülünü alır.
lm(bwght$bwght ~ bwght$cigs + bwght$faminc)##
## Call:
## lm(formula = bwght$bwght ~ bwght$cigs + bwght$faminc)
##
## Coefficients:
## (Intercept) bwght$cigs bwght$faminc
## 116.97413 -0.46341 0.09276Ama lm fonksiyonunun özetini (summary) çıkarınca, katsayılarımızın anlamlılık düzeylerinide gözlemleyebiliyorduk.
summary(lm(bwght$bwght ~ bwght$cigs + bwght$faminc))##
## Call:
## lm(formula = bwght$bwght ~ bwght$cigs + bwght$faminc)
##
## Residuals:
## Min 1Q Median 3Q Max
## -96.061 -11.543 0.638 13.126 150.083
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 116.97413 1.04898 111.512 < 2e-16 ***
## bwght$cigs -0.46341 0.09158 -5.060 4.75e-07 ***
## bwght$faminc 0.09276 0.02919 3.178 0.00151 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 20.06 on 1385 degrees of freedom
## Multiple R-squared: 0.0298, Adjusted R-squared: 0.0284
## F-statistic: 21.27 on 2 and 1385 DF, p-value: 7.942e-10Fonksiyonumuzun yetersiz olduğu görülmektedir. Dolayısıyla geliştirebiliriz.
İstersek standart hatalar için yeni bir fonksiyon oluşturabiliriz. Regresyonun standart hatalar formülü aşağıdaki gibi yazılabilir.
Dogrusal_reg_se <- function(X, Y)
sqrt(as.numeric(var(Y - X%*%solve(t(X) %*% X)%*%t(X) %*% Y))*solve(t(X) %*% X))Dogrusal_reg_se(X=Xler, Y=Y)## Warning in sqrt(as.numeric(var(Y - X %*% solve(t(X) %*% X) %*% t(X) %*% : NaNs
## üretimi## [,1] [,2] [,3]
## [1,] 0.09151077 0.02149117 NaN
## [2,] 0.02149117 0.02916682 NaN
## [3,] NaN NaN 1.048228Bu formül bize varyans kovaryans matrisi dediğimiz bir matris verdi. Aslında standart hatalar için bu matrisin diyagonal elemetlerini almamız yeterli olur. Eğer ikisini birden göstereceğim yeni bir fonksiyon yazarsak daha pratik bir fonksiyon elde edebiliriz. Hatta varyansı bulmak için var komutunu bile kullanmak zorunda değiliz. Onun içinde ayrı bir fonksiyon oluşturmak yerine, varyansı bulabileceğimiz fonsksiyonu, kendi fonksiyonumuz içine yazabiliriz.
dr <- function(X, Y) {
A <- solve(t(X) %*% X)%*%t(X) %*% Y
B <- diag(sqrt(as.numeric(t(Y-as.numeric(t(Y)%*%rep(1, nrow(Y))/nrow(Y)))%*%(Y-as.numeric(t(Y)%*%rep(1, nrow(Y))/nrow(Y)))/(nrow(Y)-1))*solve(t(X) %*% X)))
result <- list(t(A), B)
names(result) <- c("katsayılar", "Standart Hatalar")
suppressMessages(return(result))
}Yeni fonksiyonumuzu test edelim.
dr(X=Xler, Y=Y)## Warning in sqrt(as.numeric(t(Y - as.numeric(t(Y) %*% rep(1, nrow(Y))/nrow(Y)))
## %*% : NaNs üretimi## $katsayılar
## [,1] [,2] [,3]
## [1,] -0.4634075 0.09276474 116.9741
##
## $`Standart Hatalar`
## [1] 0.09290577 0.02961144 1.06420678Artık standart hatalarını da gözlemleyebileceğim bir regresyon fonksiyonumuz var.
Fonksiyona t ve pr değerlerini ekleyelim.
nrow() fonksiyonu, bir veri çerçevesindeki satır sayısını döndürür.
dr <- function(X, Y) {
A <- solve(t(X) %*% X)%*%t(X) %*% Y
B <- diag(sqrt(as.numeric(t(Y-as.numeric(t(Y)%*%rep(1, nrow(Y))/nrow(Y)))%*%(Y-as.numeric(t(Y)%*%rep(1, nrow(Y))/nrow(Y)))/(nrow(Y)-1))*solve(t(X) %*% X)))
C <- t(A)/B
D <- 1-pt(abs(C), nrow(X) - ncol(X) - 1)
result <- list(t(A), B, C, D)
names(result) <- c("katsayılar", "Standart Hatalar", "t değerleri", "pr (-t < & > t)")
suppressMessages(return(result))
}dr(X=Xler, Y=Y)## Warning in sqrt(as.numeric(t(Y - as.numeric(t(Y) %*% rep(1, nrow(Y))/nrow(Y)))
## %*% : NaNs üretimi## $katsayılar
## [,1] [,2] [,3]
## [1,] -0.4634075 0.09276474 116.9741
##
## $`Standart Hatalar`
## [1] 0.09290577 0.02961144 1.06420678
##
## $`t değerleri`
## [,1] [,2] [,3]
## [1,] -4.987931 3.132733 109.9167
##
## $`pr (-t < & > t)`
## [,1] [,2] [,3]
## [1,] 3.437984e-07 0.000884092 0Havuzlanmış Yatay Kesitler: Basit Panel Veri Yöntemleri:
Panel veri hem yatay kesit, hem de zaman serisi boyutlarına sahiptir. Örnek:Zaman İçinde Kadınların Doğurganlık oranları Sander (1992) tarafından kullanılana benzeyen fertil1 veri seti, Milli Düşünce Araştırma Merkezi’nin Genel Sosyal Araştırması’nın 1972’den 1984’e kadar çift sayılı yıllardından oluşmuştur.
Wooldridge, bu verileri, bir kadının doğurduğu toplam çocuk sayısını (kids) açıklamak için kullanıyor.
Soru: Gözlemlenebilen faktörleri kontrol ettikten sonra zaman içinde doğurganlık nasıl değişmiştir?
Kontrol altında tutulan değişkenler: eğitim görülen yıl sayısı, yaş, ırk, 16 yaşında ülkenin hangi bölgesinde yaşadığı ve 16 yaşındaki yaşam çevresidir.
Her zaman yaptığımız gibi önce veri setini indirelim, daha sonra veriyi tablo halinde gösterelim
wooldridge: data komutuyla fertil1 verisetini indirebiliriz
rmarkdown: page_table komutuyla veri setini rmarkdown üzerinde gösterebiliriz
library(wooldridge)
library(rmarkdown)
data("fertil1")
paged_table(fertil1)tabloda year değişkenini görüyoruz. Her yıl için kids değişkeninin ortalama olarak nasıl değiştiğini gösterebiliriz.
dplyr, veri işleme dilbilgisi sağlayan ve veri bilimi analistlerinin en yaygın veri işlemeyi çözmelerine yardımcı olan en çok kullanılan fiilleri sağlayan bir pakettir.
require(dplyr)## Zorunlu paket yükleniyor: dplyr##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, uniongroup_by(year) fonksiyonu, yılları tek tek almak için kullanılır
fertil1 %>%
group_by(year) %>%
summarise("ortalama çocuk sayısı (kids)" = mean(kids))## # A tibble: 7 × 2
## year `ortalama çocuk sayısı (kids)`
## <int> <dbl>
## 1 72 3.03
## 2 74 3.21
## 3 76 2.80
## 4 78 2.80
## 5 80 2.82
## 6 82 2.40
## 7 84 2.24Ortalama çocuk sayısı beklediğimiz gibi düşüyor gibi görünüyor.
require komutu, library komutu yerine kullanacağınız bir komuttur. Eğer kullandığınız paketi, daha önce yükleyip yüklemediğinizi hatırlamıyorsanız, require parentez içine yazılan paketi eğer yüklenmemişse önce yükler daha sonra library yapar.
dplyr paketi, verilerimizi düzenlemek ve temizlemek için sıklıkla kullanacağımız paketlerden biridir. %>% operatörü sayesinde sırasıyla işlem yapabiliriz.
%>% kullanarak cümle yazar gibi komut yazabilirsiniz. Yukarda bulunan chunk’da da görebileceğiniz gibi önce fertil1 veri setini yazdık, %>% yazdıktan sonra bu veri setini senelere (year) göre gruplandırdık. Bunun için group_by komutunu kullandık.
üçüncü satırda senelere göre gruplandırdığımız fertil1 veri setinde kids değişkeninin ortalamasını aldıkç Bunu yapmak için mean komutunu kullandık.
Sonuçlarımızı göstermek için, summarise komutunu kullandık ve ortalamalara ne isim vermek istiyorsak tırnak içinde onu yazdık. (“ortalama çocuk sayısı (kids)”)
diyelim ki bütün değişkenlerin yıllara göre ortalamasını görmek istiyorsunuz. across(everything() bu işlemi sizin için kolaylıkla yapacaktır.
require(dplyr)
fertil1 %>%
group_by(year) %>%
summarise(across(everything(), mean))## # A tibble: 7 × 27
## year educ meduc feduc age kids black east northcen west farm
## <int> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 72 12.2 8.33 8.90 44.9 3.03 0.0833 0.333 0.231 0.128 0.173
## 2 74 12.3 8.94 9.29 44.1 3.21 0.0578 0.237 0.353 0.110 0.208
## 3 76 12.2 8.25 8.99 43.5 2.80 0.0461 0.263 0.316 0.0855 0.237
## 4 78 12.6 9.07 9.80 43.4 2.80 0.0490 0.273 0.329 0.105 0.203
## 5 80 12.9 9.40 9.95 43.7 2.82 0.0704 0.141 0.394 0.155 0.218
## 6 82 13.2 9.56 10.2 43.2 2.40 0.199 0.231 0.290 0.0806 0.167
## 7 84 13.3 10.2 10.7 41.8 2.24 0.0678 0.260 0.333 0.102 0.192
## # ℹ 16 more variables: othrural <dbl>, town <dbl>, smcity <dbl>, y74 <dbl>,
## # y76 <dbl>, y78 <dbl>, y80 <dbl>, y82 <dbl>, y84 <dbl>, agesq <dbl>,
## # y74educ <dbl>, y76educ <dbl>, y78educ <dbl>, y80educ <dbl>, y82educ <dbl>,
## # y84educ <dbl>Diyelim ki, kids değişkeni için sadece ortalamaları değil yıllara göre standart sapmaları da görmek istiyorsunuz, bu durumda sd komutunu kullanabilirsiniz.
require(dplyr)
fertil1 %>%
group_by(year) %>%
summarise("ortalama çocuk sayısı (kids)" = mean(kids), "stadandart sapma çocuk sayısı (kids)" = sd(kids) )## # A tibble: 7 × 3
## year `ortalama çocuk sayısı (kids)` `stadandart sapma çocuk sayısı (kids)`
## <int> <dbl> <dbl>
## 1 72 3.03 1.83
## 2 74 3.21 1.50
## 3 76 2.80 1.66
## 4 78 2.80 1.58
## 5 80 2.82 1.58
## 6 82 2.40 1.70
## 7 84 2.24 1.51İlk olarak ele aldığımız örneğe dönersek kadın doğurganlığının belirleyicileri ile ilgili regresyonu yapalım.
Panel verilerle ilgili başka bir paket kullanarak regresyon analizini gerçekleştireceğiz. Ama bu ilk bir kaç örnek için bildiğimiz lm fonksiyonunu kullanacağız.
summary(lm(kids ~ educ + age + I(age^2) + black + east + northcen + west + farm + othrural + town + smcity + y74 + y76 + y78 + y80 + y82 + y84, data = fertil1))##
## Call:
## lm(formula = kids ~ educ + age + I(age^2) + black + east + northcen +
## west + farm + othrural + town + smcity + y74 + y76 + y78 +
## y80 + y82 + y84, data = fertil1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.9878 -1.0086 -0.0767 0.9331 4.6548
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -7.742457 3.051767 -2.537 0.011315 *
## educ -0.128427 0.018349 -6.999 4.44e-12 ***
## age 0.532135 0.138386 3.845 0.000127 ***
## I(age^2) -0.005804 0.001564 -3.710 0.000217 ***
## black 1.075658 0.173536 6.198 8.02e-10 ***
## east 0.217324 0.132788 1.637 0.101992
## northcen 0.363114 0.120897 3.004 0.002729 **
## west 0.197603 0.166913 1.184 0.236719
## farm -0.052557 0.147190 -0.357 0.721105
## othrural -0.162854 0.175442 -0.928 0.353481
## town 0.084353 0.124531 0.677 0.498314
## smcity 0.211879 0.160296 1.322 0.186507
## y74 0.268183 0.172716 1.553 0.120771
## y76 -0.097379 0.179046 -0.544 0.586633
## y78 -0.068666 0.181684 -0.378 0.705544
## y80 -0.071305 0.182771 -0.390 0.696511
## y82 -0.522484 0.172436 -3.030 0.002502 **
## y84 -0.545166 0.174516 -3.124 0.001831 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.555 on 1111 degrees of freedom
## Multiple R-squared: 0.1295, Adjusted R-squared: 0.1162
## F-statistic: 9.723 on 17 and 1111 DF, p-value: < 2.2e-16Temel yıl 1972 olarak alınmıştır. Doğurganlık için 82 ve 84’de keskin bir düşüş gözleniyor. Örneğin, diğer faktörler sabitken bir kadının 1982’de 1972’ye göre ortalama 0.52 daha az çocuğa sahip olduğu anlamına gelmektedir.
Bu çok büyük bir azalıştır. 1982’deki 100 kadının 1972’deki karşılaştırılabilir 100 kadından 52 tane daha az çocuğu olacağı öngörülmektedir. Örneğin bu azalış, eğitim kontrol edildiği için ortalama eğitim düzeyindeki artışdan kaynaklanan doğurganlıktaki azalıştan ayrıdır.
Bu kukla değişkenlerin katsayıları, regresyona koyduğumuz diğer değişkenlerin yakalayamadığı nedenlerden meydana gelen azalışları temsil eder.
Cinsiyet Ayrımı ve Eğitimin Getirisindeki Değişiklikler
Regresyon modeli, 1975 (temel yıl) ve 1985 yılları arasında havuzlanmış bir log(wage) denklemiyle başlıyor. (wage: saatlik ücret)
\[ log(wage) = \beta_0 + \delta_0 y85 + \beta_1 educ + \delta_1 y85 \cdot educ + \beta_2 exper + \beta_3 exper^2 + \beta_4 union + \beta_5 female + \delta_2 y85 \cdot female + u \]
union, kişi bir sendikaya üyeyse bir diğer durumlarda sıfıra eşit olan bir kukla değişkendir.
y85, gözlem 1985 yılına aitse bire, 1978’e aitse sıfıra eşit olan bir kukla değişkendir.
Veriye göz atalım, cps78_85 verisetini kullanıyoruz.
data("cps78_85")
paged_table(cps78_85)Bu verisetinde kaçtane 1978’e ait gözlem var, kaç tane 1985’e ait gözlem var. dplyr paketinin komutlarını kullanarak bulabilirim. Paketi yukarıda bir kere library() yaptığım için tekrar yapmama gerek yok.
cps78_85 %>%
group_by(year) %>%
summarise(n = n())## # A tibble: 2 × 2
## year n
## <int> <int>
## 1 78 550
## 2 85 53478 yılı için 550 kişi, 85 yılı için 534 kişilik farklı grup vardır.
Yazdığımız modelin sonuçlarına bakalım.
Bunun için summary() fonksiyonu kullanabiliriz.
summary(lm(lwage ~ y85*(educ + female) + exper + I(exper^2) + union, data = cps78_85))##
## Call:
## lm(formula = lwage ~ y85 * (educ + female) + exper + I(exper^2) +
## union, data = cps78_85)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.56098 -0.25828 0.00864 0.26571 2.11669
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.589e-01 9.345e-02 4.911 1.05e-06 ***
## y85 1.178e-01 1.238e-01 0.952 0.3415
## educ 7.472e-02 6.676e-03 11.192 < 2e-16 ***
## female -3.167e-01 3.662e-02 -8.648 < 2e-16 ***
## exper 2.958e-02 3.567e-03 8.293 3.27e-16 ***
## I(exper^2) -3.994e-04 7.754e-05 -5.151 3.08e-07 ***
## union 2.021e-01 3.029e-02 6.672 4.03e-11 ***
## y85:educ 1.846e-02 9.354e-03 1.974 0.0487 *
## y85:female 8.505e-02 5.131e-02 1.658 0.0977 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4127 on 1075 degrees of freedom
## Multiple R-squared: 0.4262, Adjusted R-squared: 0.4219
## F-statistic: 99.8 on 8 and 1075 DF, p-value: < 2.2e-16Sonuçları modelimizde yerine yazalım.
log(wage) = 0.459 + 0.118 y85 + 0.0747 educ + 0.0185 y85 educ + 0.0296 exper - 0.0004 exper^2 + 0.202 union - 0.317 female + 0.085 y85 female + u
Regresyon sonuçlarını yorumlayalım.
1978 için kesim parametresi \(\beta\)’dır. (0.459)
1985 için kesim parametresi \(\beta_0\) + \(\delta_0\)’dır. (0.459 + 0.118 = 0.577) olur.
1978 için eğitimin getirisi \(\beta_1\)’dir. (0.0747)
1985 için eğitimin getirisi \(\beta_1\) + \(\delta_1\)’dir. (0.0747 + 0.0185 = 0.0932)
Burada \(\delta_1\)=0,0185 yedi yıllık dönemde ekstra bir yıllık eğitimin getirisinin nasıl değiştiğini ölçmektedir.
1978 yılında kadınlar ve erkekler arasındaki maaş farkı \(\beta_5\)’dir. (-0.317)
1978 yılında kadınlar ve erkekler arasındaki maaş farkı \(\beta_5\) + \(\delta_2\) ’dir. (-0.317 + 0.085 = -0.232)
Cinsiyet farkına yedi yıllık dönemde bir iyileşme olup olmadığı hipotezini (H0:\(\delta_2\)=0) test edebiliriz. Örneğimizde bu katsayı için t değeri 1.658’dir ve anlamlı değildir. H0 hipotezini reddedemeyiz. Bu katsayı bizim için 0’dan farksızdır, dolayısıyla 7 yıl içinde kadın maaşlarında bir gelişme olmamıştır.
Modelde deneyim (exper) ve sendika (union) 85 yılıyla çarpılmamış. Bu iki zaman döneminde de bu değişkenlerin maaşlar üzerindeki etkisinin farksız olduğunun varsayılması sonucunda yapılmıştır.
Sayfa 448’de nominal ücretlerin, reel ücretlere çevrilmesinin neden bu çalışma için gerekli olmadığı anlatılmıştır. Lütfen bu sayfayı okuyun ve anlamaya çalışın. Bu durumun bilgisayar alıştırmalarında sorulduğu söylenmiş. Gerekirse bilgisayar alıştırması B13.2 de bu durumu doğrulayın.
Bu katsayıların nasıl yorumlanacağından daha önce bahsetmiştik. Örneğin, 1978 yılında eğitimin getirisi yaklaşık %7.5 olduğu tahmin edilmiştir. 1985 yılında eğitimin getirisi yaklaşık 1.85 yüzdelik puan artarak %9.35’e çıkmıştır. Bu katsayı t istatistiğine (1.97) göre %5 seviyesinde istatistiksel olarak anlamlıdır.
1978’de, diğer herşey sabitken bir kadın, bir erkekten yaklaşık %31.7 daha az kazanmıştır. 1985’te ücretlerdeki fark yaklaşık 8.5 yüzdelik puan azalmıştır.