Part 1: Linear Regression Recreation with challenger2
## Example: Space Shuttle Launch Data
launch <- read.csv("challenger2.csv")
# estimate beta manually
b <- cov(launch$temperature, launch$distress_ct) / var(launch$temperature)
b
[1] -0.03364796
# calculating the covariance between two variables
# estimate alpha manually
a <- mean(launch$distress_ct) - b * mean(launch$temperature)
a
[1] 2.814585
# calculates the difference between the mean of distress_ct and the product of b and the mean of temperature
# calculate the correlation of launch data
r <- cov(launch$temperature, launch$distress_ct) /
(sd(launch$temperature) * sd(launch$distress_ct))
# the var (temp) we chose to explain these stress events in 33%, we need more variables
r
[1] -0.3359996
cor(launch$temperature, launch$distress_ct)
[1] -0.3359996
# 33%
# computing the slope using correlation
r * (sd(launch$distress_ct) / sd(launch$temperature))
[1] -0.03364796
# calculates the slope of the regression line that best fits the relationship between temperature and distress_ct
# confirming the regression line using the lm function (not in text)
model <- lm(distress_ct ~ temperature, data = launch)
model
Call:
lm(formula = distress_ct ~ temperature, data = launch)
Coefficients:
(Intercept) temperature
2.81458 -0.03365
# creates a linear regression model between the temperature and distress_ct variables in the launch data frame.
summary(model)
Call:
lm(formula = distress_ct ~ temperature, data = launch)
Residuals:
Min 1Q Median 3Q Max
-1.0649 -0.4929 -0.2573 0.3052 1.7090
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.81458 1.24629 2.258 0.0322 *
temperature -0.03365 0.01815 -1.854 0.0747 .
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.7076 on 27 degrees of freedom
Multiple R-squared: 0.1129, Adjusted R-squared: 0.08004
F-statistic: 3.436 on 1 and 27 DF, p-value: 0.07474
# provides a summary with statistics
# creating a simple multiple regression function
reg <- function(y, x) {
x <- as.matrix(x)
x <- cbind(Intercept = 1, x)
b <- solve(t(x) %*% x) %*% t(x) %*% y
colnames(b) <- "estimate"
print(b)
}
# This function can be used to fit a linear regression model between the x and y variables
# examine the launch data
str(launch)
'data.frame': 29 obs. of 4 variables:
$ distress_ct : int 0 1 0 0 0 0 0 0 1 1 ...
$ temperature : int 66 70 69 68 67 72 73 70 57 63 ...
$ field_check_pressure: int 50 50 50 50 50 50 100 100 200 200 ...
$ flight_num : int 1 2 3 4 5 6 7 8 9 10 ...
# provides the structure of the launch data frame
# test regression model with simple linear regression
reg(y = launch$distress_ct, x = launch[2])
estimate
Intercept 2.81458456
temperature -0.03364796
# use regression model with multiple regression
reg(y = launch$distress_ct, x = launch[2:4])
estimate
Intercept 2.239817e+00
temperature -3.124185e-02
field_check_pressure -2.586765e-05
flight_num 2.762455e-02
# calculates the coefficients of the linear regression model
# confirming the multiple regression result using the lm function (not in text)
model <- lm(distress_ct ~ temperature + field_check_pressure + flight_num, data = launch)
model
Call:
lm(formula = distress_ct ~ temperature + field_check_pressure +
flight_num, data = launch)
Coefficients:
(Intercept) temperature field_check_pressure flight_num
2.240e+00 -3.124e-02 -2.587e-05 2.762e-02
LS0tCnRpdGxlOiAiSW4gY2xhc3MgYWN0aXZpdHkgNyIKb3V0cHV0OiBodG1sX25vdGVib29rCi0tLQoKIyMjIyBQYXJ0IDE6IExpbmVhciBSZWdyZXNzaW9uIFJlY3JlYXRpb24gd2l0aCBjaGFsbGVuZ2VyMgoKYGBge3J9CiMjIEV4YW1wbGU6IFNwYWNlIFNodXR0bGUgTGF1bmNoIERhdGEKbGF1bmNoIDwtIHJlYWQuY3N2KCJjaGFsbGVuZ2VyMi5jc3YiKQpgYGAKCmBgYHtyfQojIGVzdGltYXRlIGJldGEgbWFudWFsbHkKYiA8LSBjb3YobGF1bmNoJHRlbXBlcmF0dXJlLCBsYXVuY2gkZGlzdHJlc3NfY3QpIC8gdmFyKGxhdW5jaCR0ZW1wZXJhdHVyZSkKYgojIGNhbGN1bGF0aW5nIHRoZSBjb3ZhcmlhbmNlIGJldHdlZW4gdHdvIHZhcmlhYmxlcwpgYGAKCmBgYHtyfQojIGVzdGltYXRlIGFscGhhIG1hbnVhbGx5CmEgPC0gbWVhbihsYXVuY2gkZGlzdHJlc3NfY3QpIC0gYiAqIG1lYW4obGF1bmNoJHRlbXBlcmF0dXJlKQphCiMgY2FsY3VsYXRlcyB0aGUgZGlmZmVyZW5jZSBiZXR3ZWVuIHRoZSBtZWFuIG9mIGRpc3RyZXNzX2N0IGFuZCB0aGUgcHJvZHVjdCBvZiBiIGFuZCB0aGUgbWVhbiBvZiB0ZW1wZXJhdHVyZQpgYGAKCmBgYHtyfQojIGNhbGN1bGF0ZSB0aGUgY29ycmVsYXRpb24gb2YgbGF1bmNoIGRhdGEKciA8LSBjb3YobGF1bmNoJHRlbXBlcmF0dXJlLCBsYXVuY2gkZGlzdHJlc3NfY3QpIC8KICAgICAgIChzZChsYXVuY2gkdGVtcGVyYXR1cmUpICogc2QobGF1bmNoJGRpc3RyZXNzX2N0KSkKIyB0aGUgdmFyICh0ZW1wKSB3ZSBjaG9zZSB0byBleHBsYWluIHRoZXNlIHN0cmVzcyBldmVudHMgaW4gMzMlLCB3ZSBuZWVkIG1vcmUgdmFyaWFibGVzCnIKYGBgCgpgYGB7cn0KY29yKGxhdW5jaCR0ZW1wZXJhdHVyZSwgbGF1bmNoJGRpc3RyZXNzX2N0KQojIDMzJSAKYGBgCgoKYGBge3J9CiMgY29tcHV0aW5nIHRoZSBzbG9wZSB1c2luZyBjb3JyZWxhdGlvbgpyICogKHNkKGxhdW5jaCRkaXN0cmVzc19jdCkgLyBzZChsYXVuY2gkdGVtcGVyYXR1cmUpKQojICBjYWxjdWxhdGVzIHRoZSBzbG9wZSBvZiB0aGUgcmVncmVzc2lvbiBsaW5lIHRoYXQgYmVzdCBmaXRzIHRoZSByZWxhdGlvbnNoaXAgYmV0d2VlbiB0ZW1wZXJhdHVyZSBhbmQgZGlzdHJlc3NfY3QKYGBgCgpgYGB7cn0KIyBjb25maXJtaW5nIHRoZSByZWdyZXNzaW9uIGxpbmUgdXNpbmcgdGhlIGxtIGZ1bmN0aW9uIChub3QgaW4gdGV4dCkKbW9kZWwgPC0gbG0oZGlzdHJlc3NfY3QgfiB0ZW1wZXJhdHVyZSwgZGF0YSA9IGxhdW5jaCkKbW9kZWwKIyBjcmVhdGVzIGEgbGluZWFyIHJlZ3Jlc3Npb24gbW9kZWwgYmV0d2VlbiB0aGUgdGVtcGVyYXR1cmUgYW5kIGRpc3RyZXNzX2N0IHZhcmlhYmxlcyBpbiB0aGUgbGF1bmNoIGRhdGEgZnJhbWUuCmBgYAoKYGBge3J9CnN1bW1hcnkobW9kZWwpCiMgcHJvdmlkZXMgYSBzdW1tYXJ5IHdpdGggc3RhdGlzdGljcwpgYGAKCmBgYHtyfQojIGNyZWF0aW5nIGEgc2ltcGxlIG11bHRpcGxlIHJlZ3Jlc3Npb24gZnVuY3Rpb24KcmVnIDwtIGZ1bmN0aW9uKHksIHgpIHsKICB4IDwtIGFzLm1hdHJpeCh4KQogIHggPC0gY2JpbmQoSW50ZXJjZXB0ID0gMSwgeCkKICBiIDwtIHNvbHZlKHQoeCkgJSolIHgpICUqJSB0KHgpICUqJSB5CiAgY29sbmFtZXMoYikgPC0gImVzdGltYXRlIgogIHByaW50KGIpCn0KCiMgVGhpcyBmdW5jdGlvbiBjYW4gYmUgdXNlZCB0byBmaXQgYSBsaW5lYXIgcmVncmVzc2lvbiBtb2RlbCBiZXR3ZWVuIHRoZSB4IGFuZCB5IHZhcmlhYmxlcwpgYGAKCmBgYHtyfQojIGV4YW1pbmUgdGhlIGxhdW5jaCBkYXRhCnN0cihsYXVuY2gpCiMgcHJvdmlkZXMgdGhlIHN0cnVjdHVyZSBvZiB0aGUgbGF1bmNoIGRhdGEgZnJhbWUKYGBgCgoKYGBge3J9CiMgdGVzdCByZWdyZXNzaW9uIG1vZGVsIHdpdGggc2ltcGxlIGxpbmVhciByZWdyZXNzaW9uCnJlZyh5ID0gbGF1bmNoJGRpc3RyZXNzX2N0LCB4ID0gbGF1bmNoWzJdKQpgYGAKCmBgYHtyfQojIHVzZSByZWdyZXNzaW9uIG1vZGVsIHdpdGggbXVsdGlwbGUgcmVncmVzc2lvbgpyZWcoeSA9IGxhdW5jaCRkaXN0cmVzc19jdCwgeCA9IGxhdW5jaFsyOjRdKQojIGNhbGN1bGF0ZXMgdGhlIGNvZWZmaWNpZW50cyBvZiB0aGUgbGluZWFyIHJlZ3Jlc3Npb24gbW9kZWwKYGBgCgpgYGB7cn0KIyBjb25maXJtaW5nIHRoZSBtdWx0aXBsZSByZWdyZXNzaW9uIHJlc3VsdCB1c2luZyB0aGUgbG0gZnVuY3Rpb24gKG5vdCBpbiB0ZXh0KQptb2RlbCA8LSBsbShkaXN0cmVzc19jdCB+IHRlbXBlcmF0dXJlICsgZmllbGRfY2hlY2tfcHJlc3N1cmUgKyBmbGlnaHRfbnVtLCBkYXRhID0gbGF1bmNoKQptb2RlbApgYGAKCg==