Psych 250C - Problem Set 2 - 305946624

---

Question 1

---

One researcher (Dr. Boom) is convinced that kids born in the early 80s are more successful than those born in the mid-80s. He believes that this effect is “generational” where adults born in the early 80s will perform better in life than those born in the mid-80s. Based on this data, Dr. Boom examines the relationship between year born and income.

---

1.1) Create a plot with YearBorn on the X-axis and Income2008 on the Y-axis.

# Plot

library(ggplot2)

NLSdata = NLSdata

ggplot(NLSdata, aes(x = YearBorn, y = Income2008)) +
  geom_point() +
  geom_smooth(method = "lm", se = TRUE, color = "red") +
  theme_minimal(base_family = "Times New Roman") +
  theme(
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    axis.text = element_text(size = 12, face = "plain", family = "Times New Roman"),
    axis.title = element_text(size = 13, face = "plain", family = "Times New Roman"),
    plot.title = element_text(size = 16, face = "bold", family = "Times New Roman", hjust = 0.5)
  ) +
  labs(
    x = "Year of Birth",
    y = "Income 2008",
    title = "Income v. Birth Year"
  )

---

1.2) Use SPSS to fit a simple linear regression predicting Income2008 from YearBorn. Include code and output.

---

library(stats)

# Fit Linear Regression Model
modelIncomeBorn <- lm(Income2008 ~ YearBorn, data = NLSdata)

summary(modelIncomeBorn)

Call:
lm(formula = Income2008 ~ YearBorn, data = NLSdata)

Residuals:
   Min     1Q Median     3Q    Max 
-35556 -14406  -3214   9444  95779 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 4659987.3   475120.5   9.808   <2e-16 ***
YearBorn      -2335.6      239.7  -9.743   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 21620 on 4222 degrees of freedom
Multiple R-squared:  0.02199,   Adjusted R-squared:  0.02176 
F-statistic: 94.93 on 1 and 4222 DF,  p-value: < 2.2e-16

---

Question 3

Dr. Boom is convinced by his result and wants to write a paper published in Nature right
away! His colleague Dr. Mill has some doubts. Dr. Mill notices that the dataset also has the
ASVAB, which is the Armed Services Vocational Aptitude Battery, which is used by the US
military to determine qualifications for different positions. The ASVAB measures many
ability constructs such as general science, arithmetic, vocabulary, reading comprehension,
mathematics, etc. Higher scores on the ASVAB should suggest higher intellectual and
practical acumen (think of it as similar to an IQ test). Dr. Mill suggests they consider a
multiple linear regression where they predict Income2008 using both YearBorn and
ASVAB. Dr. Boom agrees. Dr. Boom says “Well because ASVAB is a measure of intellect and
skill, I think that when it is introduced to the model, the relationship between birth year and
income will go away. All of the variance that birth year explains in income should be shared
variance with the ASVAB. So birth year will not uniquely predict income.”

---

3.2) Fit a multiple linear regression model predicting Income2008 using both YearBorn
and ASVAB. Include your code and output.

---

# Fit Multiple Linear Regression Model
modelmultipleregression <- lm(Income2008 ~ YearBorn + ASVAB, data = NLSdata)

summary(modelmultipleregression)

Call:
lm(formula = Income2008 ~ YearBorn + ASVAB, data = NLSdata)

Residuals:
   Min     1Q Median     3Q    Max 
-44998 -13702  -2832   9390 104029 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  4.689e+06  4.582e+05   10.23   <2e-16 ***
YearBorn    -2.355e+03  2.312e+02  -10.19   <2e-16 ***
ASVAB        1.997e-01  1.120e-02   17.83   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 20860 on 4221 degrees of freedom
Multiple R-squared:  0.09049,   Adjusted R-squared:  0.09006 
F-statistic:   210 on 2 and 4221 DF,  p-value: < 2.2e-16

---

Question 6

---

Dr. Boom says “It’s hard to tell how good the model is overall! I wish there was just a simple
correlation coefficient to tell me how good this model is!” Dr. Mill says, “There is!” Calculate
R from the data, provide an interpretation of the estimate, and provide a visualization of
the correlation between the observed and predicted incomes.

---

library (stats)

# Predicted Values From the Model
predicted_values_ps2 <- predict(modelmultipleregression)

# Calculate the correlation coefficient (R)
R_ps2 <- cor(NLSdata$Income2008, predicted_values_ps2)

print(R_ps2)

[1] 0.3008149

# Plot of Predicted Values

library(ggplot2)

ggplot(NLSdata, aes(x = Income2008, y = predicted_values_ps2)) + 
  geom_point() + 
  geom_abline(intercept = 0, slope = 1, color = "red", size = 1.5) +
  theme_minimal(base_family = "Times New Roman")+
  theme(
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    axis.text = element_text(size = 12, face = "plain", family = "Times New Roman"),
    axis.title = element_text(size = 13, face = "plain", family = "Times New Roman"),
    plot.title = element_text(size = 16, face = "bold", family = "Times New Roman", hjust = 0.5)
  ) +
  labs(title = "Correlation between Observed and Predicted Incomes", x = "Observed Income in 2008", y = "Predicted Income in 2008") 

---

Question 7

---

Consider the predicted values from the model. What would the model predicted Income in
2008 be for an individual born in 1984, with an ASVAB score of 75000. Show your work in
how you calculated this estimate.

---

# Dataframe for Prediction
dataframe_ps2 <- data.frame(YearBorn = 1984, ASVAB = 75000)

# Get the Predicted income
predicted_income <- predict(modelmultipleregression, newdata = dataframe_ps2)

print(predicted_income)

      1 
31238.8 

---

Question 8

---

Dr. Boom says “I like correlation coefficients very much. I find them very interpretable. Is
there a way to quantify the unique contributions of each variable to the prediction of
income using a correlation metric?” Dr. Mill says “Yes, in fact, there are two! Let’s look at
the partial and semi-partial correlations.”

---

8.1) Calculate the partial correlation between Income2008 and YearBorn, partialing out
ASVAB. Provide an estimate, an interpretation of the estimate, and a scatterplot that
aligns with the estimate (represent the two partialed variables which are used to
calculate the correlation).

---

# Using the PPCor Package

library(ppcor)

# Partial correlation
partial_correlation_income_year <- pcor.test(NLSdata$Income2008, NLSdata$YearBorn, NLSdata$ASVAB, method = "pearson")

print(partial_correlation_income_year)

# Hardcoded R 

# Regress Income2008 on ASVAB
model_income_on_asvab <- lm(Income2008 ~ ASVAB, data = NLSdata)
residuals_income <- residuals(model_income_on_asvab)

# Regress YearBorn on ASVAB
model_yearborn_on_asvab <- lm(YearBorn ~ ASVAB, data = NLSdata)
residuals_yearborn <- residuals(model_yearborn_on_asvab)

# Correlation of Residuals
partial_cor_income2008yearborn <- cor(residuals_income, residuals_yearborn)

print(partial_cor_income2008yearborn)

[1] -0.154896

partial_residuals_df <- data.frame(residuals_income = residuals_income, residuals_yearborn = residuals_yearborn)

# Plots Residuals

library(ggplot2)

ggplot(partial_residuals_df, aes(x = residuals_yearborn, y = residuals_income)) +
  geom_point() + 
  geom_smooth(method = "lm", color = "red", se = FALSE) +
  theme_minimal(base_family = "Times New Roman") +
  theme(
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    axis.text = element_text(size = 11, face = "plain", family = "Times New Roman"),
    axis.title = element_text(size = 12, face = "plain", family = "Times New Roman"),
    plot.title = element_text(size = 14, face = "bold", family = "Times New Roman", hjust = 0.5)
  ) +
  labs(x = "Residuals for Year of Birth (controlling for ASVAB)",
       y = "Residuals for Income in 2008 (controlling for ASVAB)",
       title = "Partial Correlation between Income in 2008 and Year of Birth")

---

8.2) Calculate the semipartial correlation between Income2008 and YearBorn,
partialling out ASVAB. Provide an estimate, an interpretation of the estimate, and a
scatterplot that aligns with the estimate (represent the one partialed variables which
are used to calculate the correlation).

---

# Using PPcor Package
library(ppcor)

# Semi-partial correlation
semi_partial_correlation_income_year <- spcor.test(NLSdata$Income2008, NLSdata$YearBorn, NLSdata$ASVAB, method = "pearson")

print(semi_partial_correlation_income_year)

NA

# Hardcoded R 

# Regress Income2008 on ASVAB alone
model_income_on_asvab_only <- lm(Income2008 ~ ASVAB, data = NLSdata)
residuals_from_asvab_only <- residuals(model_income_on_asvab_only)

semipartial_cor_income2008yearborn <- cor(NLSdata$YearBorn, residuals_from_asvab_only)

print(semipartial_cor_income2008yearborn)

[1] -0.1548942

semi_residuals_df <- data.frame(YearBorn = NLSdata$YearBorn, residuals_income = residuals_from_asvab_only)

# Plots Semipartial

library(ggplot2)

ggplot(semi_residuals_df, aes(x = YearBorn, y = residuals_income)) +
  geom_point() + 
  geom_smooth(method = "lm", color = "red", se = FALSE) +
  theme_minimal(base_family = "Times New Roman") +
  theme(
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    axis.text = element_text(size = 11, face = "plain", family = "Times New Roman"),
    axis.title = element_text(size = 12, face = "plain", family = "Times New Roman"),
    plot.title = element_text(size = 14, face = "bold", family = "Times New Roman", hjust = 0.5)
  ) +
  labs(x = "Year of Birth",
       y = "Residuals for Income in 2008 (controlling for ASVAB)",
       title = "Semi-Partial Correlation")

---

8.4) Calculate the partial correlation between Income2008 and ASVAB, partialing out
YearBorn. Provide an estimate, an interpretation of the estimate, and a scatterplot
that aligns with the estimate (represent the two partialed variables which are used to
calculate the correlation).

---

# Using the PPCor Package

library(ppcor)

# Partial correlation
partial_correlation_income_asvab <- pcor.test(NLSdata$Income2008, NLSdata$ASVAB, NLSdata$YearBorn, method = "pearson")

print(partial_correlation_income_asvab)

NA

# Hardcoded R 

# Regress Income2008 on YearBorn
model_income2008_on_yearborn <- lm(Income2008 ~ YearBorn, data = NLSdata)
residuals_income2008 <- residuals(model_income2008_on_yearborn)

# Regress ASVAB on YearBorn
model_ASVAB_on_yearborn <- lm(ASVAB ~ YearBorn, data = NLSdata)
residuals_ASVAB <- residuals(model_ASVAB_on_yearborn)

# Correlations on Residuals
partial_cor_income2008asvab <- cor(residuals_income2008, residuals_ASVAB)

print(partial_cor_income2008asvab)

[1] 0.26465

partial_residuals_df_2 <- data.frame(residuals_income = residuals_income2008, residuals_ASVAB = residuals_ASVAB)

# Plots Residuals

library(ggplot2)

ggplot(partial_residuals_df_2, aes(x = residuals_ASVAB, y = residuals_income2008)) +
  geom_point() + 
  geom_smooth(method = "lm", color = "red", se = FALSE) +
  theme_minimal(base_family = "Times New Roman") +
  theme(
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    axis.text = element_text(size = 11, face = "plain", family = "Times New Roman"),
    axis.title = element_text(size = 12, face = "plain", family = "Times New Roman"),
    plot.title = element_text(size = 14, face = "bold", family = "Times New Roman", hjust = 0.5)
  ) +
  labs(x = "Residuals of ASVAB (controlling for Year of Birth)",
       y = "Residuals for Income in 2008 (controlling for Year of Birth)",
       title = "Partial Correlation")

---

8.5) Calculate the semipartial correlation between Income2008 and ASVAB, partialing
out YearBorn. Provide an estimate, an interpretation of the estimate, and a
scatterplot that aligns with the estimate (represent the two partialed variables which
are used to calculate the correlation).

---

# Using the PPCor Package

library(ppcor)

# Semi-partial correlation
semi_partial_correlation_income_asvab <- spcor.test(NLSdata$Income2008, NLSdata$ASVAB, NLSdata$YearBorn, method = "pearson")

print(semi_partial_correlation_income_asvab)

NA

# Hardcoded R 

# Regress Income2008 on ASVAB alone
model_income_on_yearborn_only <- lm(Income2008 ~ YearBorn, data = NLSdata)
residuals_from_yearborn_only <- residuals(model_income_on_yearborn_only)

semipartial_cor_income2008asvab <- cor(NLSdata$ASVAB, residuals_from_yearborn_only)

print(semipartial_cor_income2008asvab)

[1] 0.2646471

semi_residuals_df_2 <- data.frame(ASVAB = NLSdata$ASVAB, residuals_yearborn = residuals_from_yearborn_only)

# Plots Semipartial

library(ggplot2)

ggplot(semi_residuals_df_2, aes(x = ASVAB, y = residuals_yearborn)) +
  geom_point() + 
  geom_smooth(method = "lm", color = "red", se = FALSE) +
  theme_minimal(base_family = "Times New Roman") +
  theme(
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    axis.text = element_text(size = 11, face = "plain", family = "Times New Roman"),
    axis.title = element_text(size = 12, face = "plain", family = "Times New Roman"),
    plot.title = element_text(size = 14, face = "bold", family = "Times New Roman", hjust = 0.5)
  ) +
  labs(x = "ASVAB Score",
       y = "Residuals for Income in 2008 (controlling for Year of Birth)",
       title = "Semi-Partial Correlation")