library(wooldridge)
library(car)
## Loading required package: carData
library(lmtest)
## Loading required package: zoo
##
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
data("discrim")
(i) OLS estimation of the initial model
model1 <- lm(log(psoda) ~ prpblck + log(income) + prppov, data = discrim)
summary(model1)
##
## 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-08
# Get two-sided p-value for β1 (prpblck coefficient)
coef_test1 <- summary(model1)$coefficients["prpblck", ]
p_value_prpblck <- coef_test1[4] # p-value for prpblck
(ii) Correlation analysis
cor_test <- cor.test(log(discrim$income), discrim$prppov)
print(cor_test)
##
## Pearson's product-moment correlation
##
## data: log(discrim$income) and discrim$prppov
## t = -31.04, df = 407, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.8650980 -0.8071224
## sample estimates:
## cor
## -0.838467
(iii) Add log(hseval) to the model
model2 <- lm(log(psoda) ~ prpblck + log(income) + prppov + log(hseval), data = discrim)
summary(model2)
##
## Call:
## lm(formula = log(psoda) ~ prpblck + log(income) + prppov + log(hseval),
## data = discrim)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.30652 -0.04380 0.00701 0.04332 0.35272
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.84151 0.29243 -2.878 0.004224 **
## prpblck 0.09755 0.02926 3.334 0.000937 ***
## log(income) -0.05299 0.03753 -1.412 0.158706
## prppov 0.05212 0.13450 0.388 0.698571
## log(hseval) 0.12131 0.01768 6.860 2.67e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.07702 on 396 degrees of freedom
## (9 observations deleted due to missingness)
## Multiple R-squared: 0.1839, Adjusted R-squared: 0.1757
## F-statistic: 22.31 on 4 and 396 DF, p-value: < 2.2e-16
(iv) Joint significance test for log(income) and prppov
# Restricted model (removing both variables)
model3 <- lm(log(psoda) ~ prpblck + log(hseval), data = discrim)
# F-test for joint significance
anova_test <- anova(model3, model2)
print(anova_test)
## Analysis of Variance Table
##
## Model 1: log(psoda) ~ prpblck + log(hseval)
## Model 2: log(psoda) ~ prpblck + log(income) + prppov + log(hseval)
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 398 2.3911
## 2 396 2.3493 2 0.041797 3.5227 0.03045 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#(V)
cat("\nResults Summary:\n")
##
## Results Summary:
cat("\n1. Initial Model Results:")
##
## 1. Initial Model Results:
cat("\n Coefficient for prpblck:", coef(model1)["prpblck"])
##
## Coefficient for prpblck: 0.07280726
cat("\n p-value for prpblck:", p_value_prpblck)
##
## p-value for prpblck: 0.0180976
cat("\n\n2. Correlation between log(income) and prppov:")
##
##
## 2. Correlation between log(income) and prppov:
cat("\n Correlation coefficient:", cor_test$estimate)
##
## Correlation coefficient: -0.838467
cat("\n p-value:", cor_test$p.value)
##
## p-value: 2.349013e-109
cat("\n\n3. Model with log(hseval):")
##
##
## 3. Model with log(hseval):
cat("\n Coefficient for log(hseval):", coef(model2)["log(hseval)"])
##
## Coefficient for log(hseval): 0.1213057
cat("\n p-value for log(hseval):", summary(model2)$coefficients["log(hseval)", 4])
##
## p-value for log(hseval): 2.668125e-11
cat("\n\n4. Joint significance test p-value for log(income) and prppov:")
##
##
## 4. Joint significance test p-value for log(income) and prppov:
cat
## function (..., file = "", sep = " ", fill = FALSE, labels = NULL,
## append = FALSE)
## {
## if (is.character(file))
## if (file == "")
## file <- stdout()
## else if (startsWith(file, "|")) {
## file <- pipe(substring(file, 2L), "w")
## on.exit(close(file))
## }
## else {
## file <- file(file, ifelse(append, "a", "w"))
## on.exit(close(file))
## }
## .Internal(cat(list(...), file, sep, fill, labels, append))
## }
## <bytecode: 0x5ca85a3989e8>
## <environment: namespace:base>