Homework 4
Robert Papshev
Problem 1
#plot library
library(ggplot2)
#load dataset
library(faraway)
# Calculating r^2 tp explain variation
model_fit = lm(gamble ~., data = teengamb)
summary(model_fit)##
## Call:
## lm(formula = gamble ~ ., data = teengamb)
##
## Residuals:
## Min 1Q Median 3Q Max
## -51.082 -11.320 -1.451 9.452 94.252
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 22.55565 17.19680 1.312 0.1968
## sex -22.11833 8.21111 -2.694 0.0101 *
## status 0.05223 0.28111 0.186 0.8535
## income 4.96198 1.02539 4.839 1.79e-05 ***
## verbal -2.95949 2.17215 -1.362 0.1803
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 22.69 on 42 degrees of freedom
## Multiple R-squared: 0.5267, Adjusted R-squared: 0.4816
## F-statistic: 11.69 on 4 and 42 DF, p-value: 1.815e-06Percentage of variation explained by predictors
\[ r^2 = 0.5267 \]
52.67% of variation in gambling is explained by the predictors sex, status, income, and verbal
#Observation with highest residual
gamble_resid <- resid(model_fit)
which.max(gamble_resid)## 24
## 24#Compute mean and median of residuals
mean(gamble_resid)## [1] 2.645638e-16median(gamble_resid)## [1] -1.451392#Compute the correlation of the residuals with the fitted values
cor(gamble_resid, fitted.values(model_fit))## [1] 4.775848e-17#Compute the correlation of the residuals with the income
cor(gamble_resid, teengamb$income)## [1] -1.988947e-17#For all other predictors held constant, what would be the difference in predicted expenditure on gambling for a male compared to a female?
0 * coef(model_fit)["sex"] - 1 * coef(model_fit)["sex"]## sex
## 22.11833Problem 2
#load dataset
setwd('/Users/robertpapshev/Downloads')
bears=read.table("bears.txt", header=TRUE)
#Select only first observation of bear
bears=bears[bears$Obs.No==1,]
#Bears in dataset
nrow(bears)## [1] 99Why is it important to only include the first observation of a bear for linear regression?
For linear regression, independece is necessary and this allows for the observations to be independent of each other. If we included multiple observations of the same bear, it would violate the assumption of independence and could lead to biased results.
#Set variables
y <- bears$Weight
x1 <- bears$Chest.G
x2 <- bears$Head.W
#Both simple regression models
model1 <- lm(y~x1)
model2 <- lm(y~x2)
#Show results of models
model1##
## Call:
## lm(formula = y ~ x1)
##
## Coefficients:
## (Intercept) x1
## -266.68 12.65model2##
## Call:
## lm(formula = y ~ x2)
##
## Coefficients:
## (Intercept) x2
## -201.70 61.48summary(model1)##
## Call:
## lm(formula = y ~ x1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -77.75 -18.32 -0.63 17.22 97.78
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -266.6770 13.2722 -20.09 <2e-16 ***
## x1 12.6462 0.3586 35.27 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 30.77 on 97 degrees of freedom
## Multiple R-squared: 0.9276, Adjusted R-squared: 0.9269
## F-statistic: 1244 on 1 and 97 DF, p-value: < 2.2e-16summary(model2)##
## Call:
## lm(formula = y ~ x2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -186.60 -40.84 -11.71 26.70 223.84
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -201.697 34.906 -5.778 9.13e-08 ***
## x2 61.482 5.372 11.446 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 74.61 on 97 degrees of freedom
## Multiple R-squared: 0.5746, Adjusted R-squared: 0.5702
## F-statistic: 131 on 1 and 97 DF, p-value: < 2.2e-16Do the estimated regression slopes suggest positive or negative relationships?
The estimated regression slopes for both models are positive, suggesting that as Chest.G and Head.W increase, Weight also increases on average.
\[ model1-slope = 12.65 \] \[ model2-slope = 61.48 \]
What are the estimated variances of the error term for the two models (you can check the result with the summary function)?
\[ sigma^2-for-model1 = 946.79 \] \[ sigma^2-for-model2 = 5566.65 \]
What are the coefficients of determination R2 for both models. What are their interpretations?
\[ r^2-model1 = .9276 \]
92.76% of the variation in Weight is explained by girth around the bear’s chest
\[ r^2-model2 = .5746 \]
57.46% of the variation in Weight is explained by the length of the bear’s head
Between x1=“Chest.G” and x2=“Head.W”, which is the better predictor for y=“Weight”? (Address this comparing the coefficients of determination R2 of the two regressions).
The better predictor for Weight is Chest.G, as it has a higher coefficient of determination (R2 = .9276) compared to Head.W (R2 = .5746). This indicates that Chest.G explains a larger proportion of the variation in Weight than Head.W does.
#Fit a multiple linear regression model with predictors x1=“Chest.G” and x2=“Head.W”.
model3 <- lm(y~x1+x2)
#Summarization of results
summary(model3)##
## Call:
## lm(formula = y ~ x1 + x2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -84.365 -17.478 2.572 18.953 100.887
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -277.6130 14.6721 -18.92 <2e-16 ***
## x1 11.9565 0.5429 22.02 <2e-16 ***
## x2 5.6343 3.3536 1.68 0.0962 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 30.48 on 96 degrees of freedom
## Multiple R-squared: 0.9297, Adjusted R-squared: 0.9283
## F-statistic: 634.9 on 2 and 96 DF, p-value: < 2.2e-16What is the estimated variance sigma 2 of the error term for the new model?
\[ sigma^2 = 929.03 \]
What are the coefficients of determination R2 for the new model. What is its interpretation?
\[ r^2 = .9297 \]
92.97% of the variation in Weight is explained by the predictors chest girth and head length.
Do you think this model is better than the one with only x1? Why?
The new model with both predictors (Chest.G and Head.W) is slightly better than the model with only Chest.G, as it has a slightly higher R2 value (0.9297 vs 0.9276). This indicates that including Head.W as an additional predictor explains a small additional proportion of the variation in Weight. However, the improvement is minimal, suggesting that Chest.G alone is already a strong predictor of Weight.
Problem 3
#plot library
library(ggplot2)
#load dataset
library(faraway)
# Calculating r^2 and residual standard error
model1 = lm(lpsa ~ lcavol, data = prostate)
# r^2 and residual standard error
summary(model1)$sigma## [1] 0.7874994summary(model1)$r.squared## [1] 0.5394319# Adding each predictor in order one at a time and recording r^2 and residual
model2 <- lm(lpsa ~ lcavol + lweight, data = prostate)
summary(model2)$sigma## [1] 0.7506469summary(model2)$r.squared## [1] 0.5859345model3 <- lm(lpsa ~ lcavol + lweight + svi, data = prostate)
summary(model3)$sigma## [1] 0.7168094summary(model3)$r.squared## [1] 0.6264403model4 <- lm(lpsa ~ lcavol + lweight + svi + lbph, data = prostate)
summary(model4)$sigma## [1] 0.7108232summary(model4)$r.squared## [1] 0.6366035model5 <- lm(lpsa ~ lcavol + lweight + svi + lbph + age,
data = prostate)
summary(model5)$sigma## [1] 0.7073054summary(model5)$r.squared## [1] 0.6441024model6 <- lm(lpsa ~ lcavol + lweight + svi + lbph + age + lcp,
data = prostate)
summary(model6)$sigma## [1] 0.7102135summary(model6)$r.squared## [1] 0.645113model7 <- lm(lpsa ~ lcavol + lweight + svi + lbph + age + lcp + pgg45,
data = prostate)
summary(model7)$sigma## [1] 0.7047533summary(model7)$r.squared## [1] 0.6544317model8 <- lm(lpsa ~ lcavol + lweight + svi + lbph + age + lcp + pgg45 + gleason,
data = prostate)
summary(model8)$sigma## [1] 0.7084155summary(model8)$r.squared## [1] 0.6547541# Plot the trends in these two statistics
models <- list(model1, model2, model3, model4,
model5, model6, model7, model8)
RSE <- sapply(models, function(x) summary(x)$sigma)
R2 <- sapply(models, function(x) summary(x)$r.squared)
RSE## [1] 0.7874994 0.7506469 0.7168094 0.7108232 0.7073054 0.7102135 0.7047533
## [8] 0.7084155R2## [1] 0.5394319 0.5859345 0.6264403 0.6366035 0.6441024 0.6451130 0.6544317
## [8] 0.6547541plot(1:8, RSE,
type = "b",
xlab = "Number of Predictors",
ylab = "Residual Standard Error",
main = "Residual Standard Error")
plot(1:8, R2,
type = "b",
xlab = "Number of Predictors",
ylab = "R-squared",
main = "R-squared")
# Plot lpsa against l cavol
plot(prostate$lcavol, prostate$lpsa,
xlab = "lcavol",
ylab = "lpsa",
main = "lpsa vs lcavol")
# Fit both regressions
reg1 <- lm(lpsa ~ lcavol, data = prostate)
abline(reg1, col ="red")
reg2 <- lm(lcavol ~ lpsa, data = prostate)
a <- coef(reg2)[1]
b <- coef(reg2)[2]
abline(a = -a/b, b = 1/b, col = "blue")
legend("topleft", legend = c("lpsa on lcavol", "lcavol on lpsa"), col = c("red", "blue"), lwd = 2)
# Intersection of lines
mean(prostate$lcavol)## [1] 1.35001mean(prostate$lpsa)## [1] 2.478387