Bu soruyu yanıtlamak için DISCRIM verilerini kullanın. Bunlar, New Jersey ve Pennsylvania’daki fast-food restoranlarındaki çeşitli ürünlerin fiyatlarına ilişkin posta kodu düzeyinde veriler ve posta kodu popülasyonunun özellikleridir. Buradaki fikir, fast-food restoranlarının siyahların daha yoğun olduğu bölgelerde daha yüksek fiyatlar talep edip etmediğini öğrenmektir. Modelimiz
\[ psoda = \beta_0 + \beta_1prpblck + \beta_2income + u \]
Modelin değişkenlerinin ve prppov değişkeninin ne anlama geldiğini yazın.
Ortalama prpblck ve income değerlerini standart sapmalarıyla birlikte bulun. prpblck ve income ölçü birimleri nelerdir?
Bu modeli OLS ile tahmin edin ve sonuçları, n ve R-kare dahil olmak üzere denklem biçiminde rapor edin. (Tahminleri raporlarken bilimsel gösterimi kullanmayın.) prpblck üzerindeki katsayıyı yorumlayın. Sizce ekonomik olarak büyük mü?
Basit regresyon
\[ psoda = \beta_0 + \beta_1prpblck + u \]
Modelini kullanarak basit regresyonu tahmin edin. Ayrımcılık etkisi income’ı kontrol ettiğiniz modele göre daha mı büyük daha mı küçük?
Gelire göre sabit fiyat esnekliğine sahip bir model daha uygun olabilir.
\[ log(psoda) = \beta_0 + \beta_1prpblck + \beta_2(income) + u \]
Modelin tahmin edin ve tahminlerini raporlayın. Eğer prpblck .20 (20 yüzde puanı) artarsa, psoda’nın tahmini yüzde değişimi ne olur? (İpucu: Cevap 2.xx’dir, burada “xx”i doldurursunuz)
Şimdi prppov değişkenini kısım e’deki regresyona ekleyin. _1’e ne olur?
log(income) ve prppov arasındaki ilişkiyi bulun. Kabaca beklediğiniz gibi mi?
Aşağıdaki ifadeyi değerlendirin: “log(income) ve prppov çok yüksek oranda ilişkili olduğundan, aynı regresyonda olmalarına gerek yoktur.”
F-testi yapın ve yorum yapın
options(repos = c(CRAN = "https://cran.r-project.org"))
install.packages("wooldridge")## package 'wooldridge' successfully unpacked and MD5 sums checked
##
## The downloaded binary packages are in
## C:\Users\Hamıdou Souleymanou\AppData\Local\Temp\Rtmp2hCZDy\downloaded_packageslibrary(wooldridge)
data(package = "wooldridge")
data("discrim", package = "wooldridge")
head(discrim)## psoda pfries pentree wagest nmgrs nregs hrsopen emp psoda2 pfries2 pentree2
## 1 1.12 1.06 1.02 4.25 3 5 16.0 27.5 1.11 1.11 1.05
## 2 1.06 0.91 0.95 4.75 3 3 16.5 21.5 1.05 0.89 0.95
## 3 1.06 0.91 0.98 4.25 3 5 18.0 30.0 1.05 0.94 0.98
## 4 1.12 1.02 1.06 5.00 4 5 16.0 27.5 1.15 1.05 1.05
## 5 1.12 NA 0.49 5.00 3 3 16.0 5.0 1.04 1.01 0.58
## 6 1.06 0.95 1.01 4.25 4 4 15.0 17.5 1.05 0.94 1.00
## wagest2 nmgrs2 nregs2 hrsopen2 emp2 compown chain density crmrte state
## 1 5.05 5 5 15.0 27.0 1 3 4030 0.0528866 1
## 2 5.05 4 3 17.5 24.5 0 1 4030 0.0528866 1
## 3 5.05 4 5 17.5 25.0 0 1 11400 0.0360003 1
## 4 5.05 4 5 16.0 NA 0 3 8345 0.0484232 1
## 5 5.05 3 3 16.0 12.0 0 1 720 0.0615890 1
## 6 5.05 3 4 15.0 28.0 0 1 4424 0.0334823 1
## prpblck prppov prpncar hseval nstores income county lpsoda
## 1 0.1711542 0.0365789 0.0788428 148300 3 44534 18 0.11332869
## 2 0.1711542 0.0365789 0.0788428 148300 3 44534 18 0.05826885
## 3 0.0473602 0.0879072 0.2694298 169200 3 41164 12 0.05826885
## 4 0.0528394 0.0591227 0.1366903 171600 3 50366 10 0.11332869
## 5 0.0344800 0.0254145 0.0738020 249100 1 72287 10 0.11332869
## 6 0.0591327 0.0835001 0.1151341 148000 2 44515 18 0.05826885
## lpfries lhseval lincome ldensity NJ BK KFC RR
## 1 0.05826885 11.90699 10.70401 8.301521 1 0 0 1
## 2 -0.09431065 11.90699 10.70401 8.301521 1 1 0 0
## 3 -0.09431065 12.03884 10.62532 9.341369 1 1 0 0
## 4 0.01980261 12.05292 10.82707 9.029418 1 0 0 1
## 5 NA 12.42561 11.18840 6.579251 1 1 0 0
## 6 -0.05129331 11.90497 10.70358 8.394799 1 1 0 0help(discrim)Format
A data.frame with 410 observations on 37 variables:
psoda: price of medium soda, 1st wave
pfries: price of small fries, 1st wave
pentree: price entree (burger or chicken), 1st wave
wagest: starting wage, 1st wave
nmgrs: number of managers, 1st wave
nregs: number of registers, 1st wave
hrsopen: hours open, 1st wave
emp: number of employees, 1st wave
psoda2: price of medium soday, 2nd wave
pfries2: price of small fries, 2nd wave
pentree2: price entree, 2nd wave
wagest2: starting wage, 2nd wave
nmgrs2: number of managers, 2nd wave
nregs2: number of registers, 2nd wave
hrsopen2: hours open, 2nd wave
emp2: number of employees, 2nd wave
compown: =1 if company owned
chain: BK = 1, KFC = 2, Roy Rogers = 3, Wendy’s = 4
density: population density, town
crmrte: crime rate, town
state: NJ = 1, PA = 2
prpblck: proportion black, zipcode
prppov: proportion in poverty, zipcode
prpncar: proportion no car, zipcode
hseval: median housing value, zipcode
nstores: number of stores, zipcode
income: median family income, zipcode
county: county label
lpsoda: log(psoda)
lpfries: log(pfries)
lhseval: log(hseval)
lincome: log(income)
ldensity: log(density)
NJ: =1 for New Jersey
BK: =1 if Burger King
KFC: =1 if Kentucky Fried Chicken
RR: =1 if Roy Rogers
mean(discrim$prpblck)## [1] NAsd(discrim$prpblck)## [1] NAmean(discrim$income)## [1] NAsd(discrim$income)## [1] NAElimizdeki mean ve sd fonksiyonlarını kullanarak ortalama ve standart sapma değerlerini hesaplamayı denedik, ancak sonuç NA çıktı. Bu durum, ilgili değişkenlerin içerisinde bazı eksik (mevcut olmayan) gözlemlerin olabileceğini gösteriyor. discrim veri seti 410 gözlemden oluştuğundan, her bir gözlemi tek tek kontrol etmemiz mümkün değil. Dolayısıyla, bu değişkenlerde kaç adet eksik gözlem bulunduğunu belirleyemiyoruz.
Bu noktada, R bize is.na fonksiyonu ile yardımcı olabilir. is.na, aslında İngilizce bir soru gibidir: “Bu gözlem eksik mi?” diye sorarız. R ise her bir gözlem için, eksik olup olmadığını TRUE (eksik) veya FALSE (eksik değil) olarak geri döner.
Bir örnek vermek gerekirse, is.na(discrim\(income)** ifadesi, **income** değişkenindeki her bir gözlemi kontrol eder. Ardından, **sum(is.na(discrim\)income)) ifadesi, eksik olan gözlemlerin sayısını toplar ve sonucu bize verir. Bu şekilde, eksik gözlem sayısını kolayca bulabiliriz.
sum(is.na(discrim$prpblck))## [1] 1sum(is.na(discrim$income))## [1] 1Gördüğünüz gibi, hem prbblck hem de income değişkenlerinin birer gözlemi eksik (boş) değere sahip. Bu nedenle, mean ve sd fonksiyonlarının bu eksik gözlemlerden etkilendiğini ve sonuçlarının NA çıktığını belirtmemiz gerekiyor.
mean(discrim$prpblck,na.rm = TRUE)## [1] 0.1134864mean(discrim$prpblck,na.rm = TRUE)## [1] 0.1134864sd(discrim$prpblck,na.rm = TRUE)## [1] 0.1824165mean(discrim$income, na.rm = TRUE)## [1] 47053.78sd(discrim$income, na.rm = TRUE)## [1] 13179.29Fonksiyonlara eklediğimiz na.rm (NA Remove, yani “eksik değerleri yok say”) argümanı, hesaplamalarda eksik değerlerin dikkate alınmamasını sağlar. Bu durumda:
prbblck değişkeni için ortalama 0.11, standart sapma ise 0.18 olarak bulunur.
income değişkeni için ortalama 47053, standart sapma ise 13179 olacaktır.
Eğer tüm değişkenler için eksik değerlerin sayısını, geçerli gözlem sayılarını, ortalamalarını ve standart sapmalarını topluca görmek isterseniz, bu konuda size vtable paketi yardımcı olabilir.
install.packages("vtable")## package 'vtable' successfully unpacked and MD5 sums checked
##
## The downloaded binary packages are in
## C:\Users\Hamıdou Souleymanou\AppData\Local\Temp\Rtmp2hCZDy\downloaded_packageslibrary(vtable)
sumtable(discrim, summ=c('notNA(x)', 'countNA(x)', 'mean(x)','sd(x)'),out='return')## Variable NotNA CountNA Mean Sd
## 1 psoda 402 8 1 0.089
## 2 pfries 393 17 0.92 0.11
## 3 pentree 398 12 1.3 0.64
## 4 wagest 390 20 4.6 0.35
## 5 nmgrs 404 6 3.4 1
## 6 nregs 388 22 3.6 1.2
## 7 hrsopen 410 0 14 2.8
## 8 emp 404 6 18 9.4
## 9 psoda2 388 22 1 0.094
## 10 pfries2 382 28 0.94 0.11
## 11 pentree2 386 24 1.4 0.65
## 12 wagest2 389 21 5 0.25
## 13 nmgrs2 404 6 3.5 1.1
## 14 nregs2 388 22 3.6 1.2
## 15 hrsopen2 399 11 14 2.8
## 16 emp2 397 13 18 8.6
## 17 compown 410 0 0.34 0.48
## 18 chain 410 0 2.1 1.1
## 19 density 409 1 4562 5132
## 20 crmrte 409 1 0.053 0.047
## 21 state 410 0 1.2 0.39
## 22 prpblck 409 1 0.11 0.18
## 23 prppov 409 1 0.071 0.067
## 24 prpncar 409 1 0.11 0.12
## 25 hseval 409 1 147399 56070
## 26 nstores 410 0 3.1 1.8
## 27 income 409 1 47054 13179
## 28 county 410 0 14 8
## 29 lpsoda 402 8 0.04 0.085
## 30 lpfries 393 17 -0.088 0.12
## 31 lhseval 409 1 12 0.39
## 32 lincome 409 1 11 0.28
## 33 ldensity 409 1 8 1
## 34 NJ 410 0 0.81 0.39
## 35 BK 410 0 0.42 0.49
## 36 KFC 410 0 0.2 0.4
## 37 RR 410 0 0.24 0.43discrimreg <- lm(psoda~prpblck+income, data = discrim)
summary(discrimreg)##
## Call:
## lm(formula = psoda ~ prpblck + income, data = discrim)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.29401 -0.05242 0.00333 0.04231 0.44322
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 9.563e-01 1.899e-02 50.354 < 2e-16 ***
## prpblck 1.150e-01 2.600e-02 4.423 1.26e-05 ***
## income 1.603e-06 3.618e-07 4.430 1.22e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.08611 on 398 degrees of freedom
## (9 observations deleted due to missingness)
## Multiple R-squared: 0.06422, Adjusted R-squared: 0.05952
## F-statistic: 13.66 on 2 and 398 DF, p-value: 1.835e-06\[ psoda = 0.956 + 0.115prpblck + 0.0000016income + u \]
Örnek büyüklüğü 399 gözlemden oluşmaktadır (398 serbestlik derecesi ve 9 eksik gözlem bulunmaktadır). Ayarlanmış R² değeri ise 0.595’tir. prpblck katsayısı, tüm diğer değişkenler sabit tutulduğunda, prpblck değişkeninde %10’luk bir artış olması durumunda soda fiyatının ekonomik açıdan anlamlı olmayan bir düzeyde yaklaşık 1,2 sent artacağını ifade etmektedir.
basitdiscrimreg <- lm(psoda~prpblck, data = discrim)
summary(basitdiscrimreg)##
## Call:
## lm(formula = psoda ~ prpblck, data = discrim)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.30884 -0.05963 0.01135 0.03206 0.44840
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.03740 0.00519 199.87 < 2e-16 ***
## prpblck 0.06493 0.02396 2.71 0.00702 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.0881 on 399 degrees of freedom
## (9 observations deleted due to missingness)
## Multiple R-squared: 0.01808, Adjusted R-squared: 0.01561
## F-statistic: 7.345 on 1 and 399 DF, p-value: 0.007015Basit regresyon analizi sonucunda prpblck değişkeninin katsayısının tahmini 0.065 olarak bulunmuştur. Bu değer, önceki tahminden daha düşüktür. Bu durum, gelir değişkeni modele dahil edilmediğinde, ayrımcılık etkisinin azaldığını göstermektedir.
logdiscrimreg <- lm(log(psoda)~prpblck+log(income), data = discrim)
summary(logdiscrimreg)##
## Call:
## lm(formula = log(psoda) ~ prpblck + log(income), data = discrim)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.33563 -0.04695 0.00658 0.04334 0.35413
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.79377 0.17943 -4.424 1.25e-05 ***
## prpblck 0.12158 0.02575 4.722 3.24e-06 ***
## log(income) 0.07651 0.01660 4.610 5.43e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.0821 on 398 degrees of freedom
## (9 observations deleted due to missingness)
## Multiple R-squared: 0.06809, Adjusted R-squared: 0.06341
## F-statistic: 14.54 on 2 and 398 DF, p-value: 8.039e-07paste( (0.2*100)*0.122, "yüzdelik artış")## [1] "2.44 yüzdelik artış"Eğer prpblck %20 artarsa, psoda tahmini olarak %2,44 oranında artış gösterecektir.
logdiscrimregprpov <- lm(log(psoda)~prpblck+log(income)+prppov, data = discrim)
summary(logdiscrimregprpov)##
## Call:
## lm(formula = log(psoda) ~ prpblck + log(income) + prppov, data = discrim)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.32218 -0.04648 0.00651 0.04272 0.35622
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.46333 0.29371 -4.982 9.4e-07 ***
## prpblck 0.07281 0.03068 2.373 0.0181 *
## log(income) 0.13696 0.02676 5.119 4.8e-07 ***
## prppov 0.38036 0.13279 2.864 0.0044 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.08137 on 397 degrees of freedom
## (9 observations deleted due to missingness)
## Multiple R-squared: 0.08696, Adjusted R-squared: 0.08006
## F-statistic: 12.6 on 3 and 397 DF, p-value: 6.917e-08Prppov değişkeninin modele eklenmesi, prpblck katsayısının 0,0738’e düşmesine neden olur.
cor(log(discrim$income), discrim$prppov, use = "complete.obs")## [1] -0.838467Korelasyon yaklaşık olarak -0.838’dir. Bu sonuç mantıklıdır, çünkü gelirdeki düşüşlerin genellikle daha yüksek yoksulluk oranlarıyla ilişkilendirilebileceği beklenebilir. Yani, gelir azaldıkça yoksulluk oranlarının artması olasıdır..
Yüksek düzeyde ilişkili olmalarına rağmen, prpblck ve prppov değişkenlerinin her ikisinin de modele dahil edilmesi mükemmel bir doğrusallık ile sonuçlanmaz. Bunun yerine, ayırt edici etkiyi izole edebilmek için başka bir kontrol değişkeni ekleyerek modeli tamamlamak daha etkili olacaktır. Bu ek kontrol değişkeni, her iki faktörün etkilerini daha net bir şekilde ayırmaya yardımcı olur.
install.packages("wooldridge")
library(wooldridge)
data("discrim")
model <- lm(psoda ~ prpblck + income, data = discrim)
r_squared <- summary(model)$r.squared
n <- nrow(discrim)
k <- 2
F_statistic <- (r_squared / k) / ((1 - r_squared) / (n - k - 1))
cat("The calculated F-statistic is:", F_statistic, "\n")## The calculated F-statistic is: 13.96573alpha <- 0.05
critical_F <- qf(1 - alpha, df1 = k, df2 = n - k - 1)
cat("The critical F-value at 5% significance level is:", critical_F, "\n")## The critical F-value at 5% significance level is: 3.017891if (F_statistic > critical_F) {
cat("Result: Reject the null hypothesis (H0). The model is jointly significant.\n")
} else {
cat("Result: Fail to reject the null hypothesis (H0). The model is not jointly significant.\n")
}## Result: Reject the null hypothesis (H0). The model is jointly significant.Eğer hesaplanan F-istatistiği kritik F-değerini aşarsa, null hipotezini reddederiz (𝐻0:𝛽1=𝛽2=0H 0 :β 1 =β 2 =0). Bu, bağımsız değişkenler prpblck ve income’un psoda değişkenindeki varyansı açıklamada ortak olarak anlamlı olduğunu gösterir.
Eğer hesaplanan F-istatistiği, kritik F-değeriyle eşit ya da küçükse, null hipotezini reddedemeyiz, yani prpblck ve income değişkenlerinin psoda üzerinde ortak bir anlamlı etkisi olmadığını söyleyebiliriz.