Lab 6

options(scipen=999)

library(socviz)
load("Catholic_viz.RData")

DATA PARTITIONING

library(caret)
set.seed(1)
index <- createDataPartition(Catholic_viz$read12, p=0.6, list=FALSE)
train <- Catholic_viz[index,]
test  <- Catholic_viz[-index,]

Descriptive Statistics

library(psych)
library(tidyverse)
library(flextable)
psych::describe(train, fast=TRUE)

         vars    n       mean         sd       min        max      range
id          1 3583 4642399.08 2719725.86 124915.00 7979086.00 7854171.00
read12      2 3583      51.56       9.39     29.15      68.09      38.94
math12      3 3583      52.00       9.59     30.14      71.37      41.23
female      4 3583       0.53       0.50      0.00       1.00       1.00
asian       5 3583       0.06       0.25      0.00       1.00       1.00
hispan      6 3583       0.11       0.31      0.00       1.00       1.00
black       7 3583       0.08       0.26      0.00       1.00       1.00
motheduc    8 3583      13.24       2.03      8.00      18.00      10.00
fatheduc    9 3583      13.55       2.26      8.00      18.00      10.00
lfaminc    10 3583      10.31       0.80      6.21      12.35       6.13
hsgrad     11 3583       0.93       0.26      0.00       1.00       1.00
cathhs     12 3583       0.07       0.26      0.00       1.00       1.00
parcath    13 3583       0.36       0.48      0.00       1.00       1.00
               se
id       45436.17
read12       0.16
math12       0.16
female       0.01
asian        0.00
hispan       0.01
black        0.00
motheduc     0.03
fatheduc     0.04
lfaminc      0.01
hsgrad       0.00
cathhs       0.00
parcath      0.01

INTERPRETATIONS

A significant majority of our observations were female while few were asian, hispanic, or black. The average years of education for both fathers and mothers was around 13 years. Majority of our observations were at least high school graduates but very few attended a catholic high school. About 30 percent of the parents reported being catholic.

Boxplots for Variables

data_long <- train %>%
  select(where(is.numeric)) %>%
  pivot_longer(cols = everything(), names_to = "Variable", values_to = "Value")

ggplot(data_long, aes(x = Variable, y = Value)) +
  geom_boxplot() +
  coord_flip() + 
  facet_wrap(~ Variable, scales = "free", ncol=3) +
  theme_minimal() +
  theme(axis.text.x = element_blank()) + 
  labs(title = "Horizontal Boxplot for Each Variable")

INTERPRETATIONS

Few outliers in the data set but majority are in the family income variable.

Histograms

library(DataExplorer)
plot_histogram(train)

COMMENTS

Both the average reading and math scores were above 50. But the average reading scores were higher than math scores.

Scatterplot

plot_scatterplot(test, by="read12")

COMMENTS

Many of these scatterplots are unable to be interpreted due to their categorical nature.

Correlation Matrix

plot_correlation(train, type="continuous")

INTERPRETATIONS

Only 3 variables have a strong correlation with another. Math scores and reading scores have a strong positive correlation with 0.73. Father and mother years of education also have a strong positive correlation with 0.58. Income with both math and reading scores has a moderate positive correlation.

Estimate Coefficients and Coefficients Table

model <- lm(log(read12) ~ female + asian + hispan +  black + motheduc + fatheduc + 
                                   lfaminc + hsgrad + cathhs + parcath + math12, data=train)
model %>% as_flextable

	Estimate	Standard Error	t value	Pr(>|t|)
(Intercept)	3.016	0.030	99.123	0.0000	***
female	0.057	0.004	12.864	0.0000	***
asian	-0.021	0.009	-2.255	0.0242	*
hispan	-0.001	0.008	-0.168	0.8663
black	-0.029	0.009	-3.341	0.0008	***
motheduc	0.001	0.001	0.383	0.7019
fatheduc	0.004	0.001	2.849	0.0044	**
lfaminc	0.009	0.003	2.644	0.0082	**
hsgrad	0.006	0.009	0.731	0.4645
cathhs	0.005	0.009	0.594	0.5528
parcath	-0.007	0.005	-1.327	0.1847
math12	0.014	0.000	52.979	0.0000	***
Signif. codes: 0 <= '***' < 0.001 < '**' < 0.01 < '*' < 0.05
Residual standard error: 0.1315 on 3571 degrees of freedom
Multiple R-squared: 0.5421, Adjusted R-squared: 0.5407
F-statistic: 384.3 on 3571 and 11 DF, p-value: 0.0000

COMMENTS

The variables that appear to be statistically significant given by their p values are Female, Asian, Black, fathers years of education, family income, and math scores.

Coefficient Magnitude Plot

library(coefplot)
coefplot(model, soft="magnitude", intercept=FALSE)

Variance Inflation Factors

library(car)
vif(model)

  female    asian   hispan    black motheduc fatheduc  lfaminc   hsgrad 
1.018173 1.037221 1.216667 1.117319 1.620753 1.736261 1.389016 1.104961 
  cathhs  parcath   math12 
1.157689 1.280675 1.353957 

Residual Analysis

Residual Range

quantile(round(residuals(model, type = "deviance"),2))

   0%   25%   50%   75%  100% 
-0.62 -0.08  0.01  0.08  0.44 

Residual Plots

plot(model)

INTPRETATIONS

Observations appear to be randomly distributed which indicates that our model is a good fit and is trustworthy.

library(broom)
train2 <- augment(model, data=train)  # this appends predicted to original dataset to build plots on your own

p <- ggplot(train2, mapping = aes(y=.fitted, x=read12))
p  + geom_point()

p <- ggplot(train2, mapping = aes(y=.resid, x=.fitted))
p  + geom_point()

Performance Evaluation

pred <- predict(model, newdata = test, interval="predict")  
test_w_predict <- cbind(test, pred)  

test_w_predict %>%
  head(10) %>%
    select(math12, read12, motheduc, fatheduc, lfaminc) %>% 
  as_flextable(show_coltype = FALSE) 

math12	read12	motheduc	fatheduc	lfaminc
49.8	61.4	14	12	10.3
50.4	59.3	14	11	10.3
54.3	57.6	12	12	10.7
56.7	52.5	12	11	11.0
54.0	42.1	14	14	10.3
63.9	64.0	14	12	10.3
61.4	38.5	14	11	9.1
62.5	60.3	16	14	10.3
58.8	63.6	12	12	9.8
59.7	52.6	14	16	10.3
n: 10

Plot Actual vs Fitted (test)

plot(test_w_predict$read12, test_w_predict$math12, pch=16, cex=1)

INTERPRETATIONS

Based on our plot, we can predict that reading scores will increase when math scores increase and vice versa.

Performance Metrics

library(Metrics)
metric_label <- c("MAE","RMSE", "MAPE")
metrics <- c(round(mae(test_w_predict$read12, test_w_predict$math12),4),
                         round(rmse(test_w_predict$read12, test_w_predict$math12),4),
                         round(mape(test_w_predict$read12, test_w_predict$math12),4))
pmtable <- data.frame(Metric=metric_label, Value = metrics)
flextable(pmtable)

Metric	Value
MAE	5.7632
RMSE	7.3525
MAPE	0.1186

Model Fit by High School Grad

p <- ggplot(data = subset(test_w_predict, hsgrad %in% c("1", "0")),
            aes(x = read12,
                y = motheduc, ymin = fatheduc, ymax = lfaminc,
                color = hsgrad,
                fill = hsgrad,
                group = hsgrad))
p +  geom_point(alpha = 0.5) + 
       geom_line() + geom_ribbon(alpha = 0.2, color = FALSE) +
       labs(title="Actual vs Fitted with Upper and Lower CI",
              subtitle="Graduated and Not Graduated",
              caption="catholic{wooldridge}") +
       xlab(NULL) + ylab("fitted") +
       theme(legend.position = "bottom")

pip <- lm(log(read12) ~ female + asian + black + fatheduc + lfaminc + hsgrad, data=train)
summary(pip)

Call:
lm(formula = log(read12) ~ female + asian + black + fatheduc + 
    lfaminc + hsgrad, data = train)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.59892 -0.11355  0.03404  0.13050  0.47214 

Coefficients:
             Estimate Std. Error t value             Pr(>|t|)    
(Intercept)  3.197836   0.039045  81.901 < 0.0000000000000002 ***
female       0.033752   0.005921   5.701 0.000000012890255207 ***
asian        0.020560   0.012138   1.694               0.0904 .  
black       -0.097986   0.011370  -8.618 < 0.0000000000000002 ***
fatheduc     0.019343   0.001471  13.153 < 0.0000000000000002 ***
lfaminc      0.034011   0.004189   8.120 0.000000000000000635 ***
hsgrad       0.110746   0.011526   9.609 < 0.0000000000000002 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.1765 on 3576 degrees of freedom
Multiple R-squared:  0.1738,    Adjusted R-squared:  0.1725 
F-statistic: 125.4 on 6 and 3576 DF,  p-value: < 0.00000000000000022

library(ggplot2)
p <- ggplot(data=train, mapping=aes(y=log(read12), x=female + asian + black + fatheduc + lfaminc + hsgrad))

p +  geom_point(alpha = 0.2) + 
       geom_smooth(method = "lm", aes(color = "OLS", fill = "OLS"))

p + geom_point(alpha=0.1) +
    geom_smooth(color = "tomato", fill="tomato", method = MASS::rlm) +
    geom_smooth(color = "steelblue", fill="steelblue", method = "lm")

p + geom_point(alpha=0.1) +
    geom_smooth(color = "tomato", method = "lm", linewidth = 1.2, 
                formula = y ~ splines::bs(x, 3), se = FALSE)

p + geom_point(alpha=0.1) +
    geom_quantile(color = "tomato", size = 1.2, method = "rqss",
                  lambda = 1, quantiles = c(0.20, 0.5, 0.85))

FINAL ASSESSMENT

A significant majority of our observations were female while few were Asian, Hispanic, or black. The average years of education for both fathers and mothers was around 13 years. Majority of our observations were at least high school graduates but very few attended a catholic high school. About 30 percent of the parents reported being catholic.
Only 3 variables had a strong correlation with another. Math scores and reading scores had a strong positive correlation. Father and mother years of education also had a strong positive correlation. Income with both math and reading scores had a moderate positive correlation.
The variables that appeared to be statistically significant given by our linear regression model were Female, Asian, Black, fathers years of education, family income, and math scores.
After testing the validity of the model, it is indicated that our model is a moderately good fit with some areas that could be improved. Residual values were randomly distributed which indicated our model was a good fit. But our performance metrics indicated some errors in our model. Visualization for this model has proved valuable has it has provided insights on not only the data that is being presented but also any possible errors that could help us to improve the visual representation of the graph.

END