Fonksiyon çalışıyor
Bu fonksiyonu birde herhangi bir veri seti üzerinden deneyelim. Geçen hafta yapmış olduğumuz, bwght verisetini kullanarak cigs ve faminc kullanarak regresyon sonuçlarını tekrarlayalım.
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. Ama adına X demek zorunda değilsiniz, Xler ismini verebiliriz.
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ırdığımız bir katsayımız daha var. Yani iki değil, üç katsayı bulacağız bu yüzden Xler veri setimize birlerden oluşan bir sütun daha eklememiz gerek.
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 1Şimdi Xler matrisimiz tamamlandı.
Y matrisini de oluşturalım.
Y = cbind(bwght$bwght)
head(Y)## [,1]
## [1,] 109
## [2,] 133
## [3,] 129
## [4,] 126
## [5,] 134
## [6,] 118Artık yeni fonksiyonumuzu kullanabiliriz.
Dogrusal_reg(X=Xler, Y=Y)
Karşılaştırabileceğimiz gibi sonuçlar aynı.
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-10Signif. 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-10 Fonksiyonumuz bu konuda yetersiz. Fonskiyonu istediğiniz kadar geliştirebilirsiniz.
İstersek standart hatalar için yeni bir fonksiyon oluşturabiliriz. Regresyonun standart hatalar formülüü hatırlayın.
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) %*% :
## Production de NaN## [,1] [,2] [,3]
## [1,] 0.09151077 0.02149117 NaN
## [2,] 0.02149117 0.02916682 NaN
## [3,] NaN NaN 1.048228Warning in sqrt(as.numeric(var(Y - X %% solve(t(X) %% X) %% t(X) %% : NaNs produced
Bu formül bize varyans kovaryans matrisi dediğimiz bir matris verdi. Aslında standart hatalar için bu matrisin diyagonal elemetlerini almam yeterli. Eğer ikisini birden göstereceğim yeni bir fonksiyon yazarsam daha pratik bir fonksiyon elde edebilirim. 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)))
## %*% : Production de NaN## $katsayılar
## [,1] [,2] [,3]
## [1,] -0.4634075 0.09276474 116.9741
##
## $`Standart Hatalar`
## [1] 0.09290577 0.02961144 1.06420678Warning in sqrt(as.numeric(t(Y - as.numeric(t(Y) %% rep(1, nrow(Y))/nrow(Y))) %% : NaNs produced
$katsayılar
Artık standart hatalarını da gözlemleyebileceğim bir regresyon fonksiyonum var.
Fonksiyona t ve pr değerlerini ekleyelim.
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)))
## %*% : Production de NaN## $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 0Panel veri hem yatay kesit, hem de zaman serisi boyutlarına sahiptir. Örnek 13.1 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
library(wooldridge) # data komutuyla fertil1 verisetini indirmek içim
library(rmarkdown) # page_table komutuyla veri setini rmarkdown üzerinde gösterebilmek için
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. (Bunu merak ettiğiniz diğer değişkenler için de yapın.)
require(dplyr)## Le chargement a nécessité le package : dplyr##
## Attachement du package : 'dplyr'## Les objets suivants sont masqués depuis 'package:stats':
##
## filter, lag## Les objets suivants sont masqués depuis 'package:base':
##
## intersect, setdiff, setequal, union