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ünü 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## $katsayilar
## [,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## $katsayilar
## [,1] [,2] [,3]
## [1,] -0.4634075 0.09276474 116.9741
##
## $`Standart Hatalar`
## [1] 0.09290577 0.02961144 1.06420678
##
## $`t degerleri`
## [,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## Warning: le package 'dplyr' a été compilé avec la version R 4.1.3## Error: le chargement du package ou de l'espace de noms a échoué pour 'dplyr' in loadNamespace(i, c(lib.loc, .libPaths()), versionCheck = vI[[i]]) :
## le package 'rlang' 1.0.1 est déjà chargé, mais >= 1.0.2 est requisOrtalama çocuk sayısı beklediğimiz gibi düşüyor gibi görünüyor.
Bir sonraki aşamaya geçmeden önce yazdığımız kodları analamaya çalışalım.
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)”)
Biraz pratik yaptıktan sonra, bu ortalamaları çıkarmak sizin içinde basit hale gelecek.
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.
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. Not: console’a help(“summarise”) yazarsanız, başka neleri bulabileceğiniz hakkında bilgi sahibi olabilirsiniz. R-studio’yu kullandıkça deneyimiziniz artacaktır.