a.

What is the expected value of 28-day strength when accelerated strength = 2500?
- Plugging 2500 into our regression formula gives \(y = 1800 + 1.3(2500) = 5050\).

b.

By how much can we expect 28-day strength to change when accelerated strength increases by 1 psi?
- Plugging 1 into our regression formula will show us an increase for accelerated strength of 1.3

c.

Answer part (b) for an increase of 100 psi.
- Following the logic of part (b) we will see an increase of 130 for accelerated strength.

d.

Answer part (b) for a decrease of 100 psi.
- This, much like part (c), will see a decrease of 130 for accelerated strength. In other words, a change of -130.

a.

What is the probability that the observed value of 28-day strength will exceed 5000 psi when the value of accelerated strength is 2000?
- We are looking for \(P(Y > 5000 | x^* = 2000)\). We know that \(\sigma = 350\) and that \(f(x) = 1800 + 1.3x\). The following calculation will be shown and then represented with code.

\[P(Y > 5000 | x^* = 2000) = P \left( Z > \frac{f(2000) - 5000}{350}\right)\] \[= 1 - \text{normalcdf(-100, -1.714)}\] \[\approx .957\]

func_ex8 <- function(x) {
    ret_val <- 1800 + (1.3 * x)
    return(ret_val)
}

z <- (func_ex8(2000) - 5000)/350
cat("The probability of observing this is approximately", round(1 -
    pnorm(z), 3))

## The probability of observing this is approximately 0.957

b.

Repeat part (a) with 2500 in place of 2000.
- The same code will be utilized again.

z <- (func_ex8(2500) - 5000)/350
cat("The probability of observing this is approximately", round(1 -
    pnorm(z), 3))

## The probability of observing this is approximately 0.443

c.

Consider making two independent observations on 28-day strength, the first for \(x = 2000\) and the second for \(x = 2500\). What is the probability that the second observation will exceed the first by more than 1000psi?
- For this we need to calculate things a little bit differently. We have the following formula:

num <- func_ex8(2500) - func_ex8(2000) - 1000
den <- sqrt(2 * 350^2)
quotient <- num/den
round(pnorm(quotient, lower.tail = FALSE), 3)

## [1] 0.76

d.

Let \(Y_1\) and \(Y_2\) denote observations on 28-day strength when \(x = x_1\) and \(x = x_2\), respectively. By how much would \(x_2\) have to exceed \(x_1\) in order that \(P(Y_2 > Y_1) = .95\)?
- This isn’t as complex as it sounds. We simply flip around our unknown variables. Instead of knowing that \(\Delta_0 = 1000\) and not knowing \(Z\), we now know that \(Z = \text{invNorm}(.95) \approx 1.64\) and do not know \(\Delta_0\). We can rewrite the formula the following way:

This shows that \(x_2\) would have to be at least 811 larger than \(x_1\) for \(P(Y_2 > Y_1) = .95\).

a.

Determine the equation of the estimated regression line.

temp <- c(170, 172, 173, 174, 174, 175, 176, 177, 180, 180, 180,
    180, 180, 181, 181, 182, 182, 182, 182, 184, 184, 185, 186,
    188)
ratio <- c(0.84, 1.31, 1.42, 1.03, 1.07, 1.08, 1.04, 1.8, 1.45,
    1.6, 1.61, 2.13, 2.15, 0.84, 1.43, 0.9, 1.81, 1.94, 2.68,
    1.49, 2.52, 3, 1.87, 3.08)

df <- data.frame(temp, ratio)

ggplot(data = df, aes(x = temp, y = ratio)) + geom_point() +
    theme_gray() + stat_smooth(method = lm, se = FALSE)

## `geom_smooth()` using formula 'y ~ x'

lm(ratio ~ temp)

## 
## Call:
## lm(formula = ratio ~ temp)
## 
## Coefficients:
## (Intercept)         temp  
##   -15.24497      0.09424

lm() is equivalent to LinReg(a + bx) on my Ti-84 Calculator. (ratio~temp) indicates that ratio here is the response variable and that temp is a predictor.

The output of my code resulted in \(\beta_0 = -15.24\) and \(\beta_1 = 0.09\). This gives the equation: \[\boxed{f(x) = -15.24 + 0.09x}\]

b.

Calculate a point estimate for true average efficiency ratio when tank temperature is 182. - Plugging 182 into our equation gets \[f(182) = -15.24 + 0.094*182 \approx \boxed{1.868}\]

c.

Calculate the values of the residuals from the least squares line for the four observations for which temperature is 182. Why do they not all have the same sign?

What we know is that we need to compare our actual observations at 182 to what 182s output would be given our equation. We calculated the second part already, \(f(182) \approx 1.868\), so we just need to pull out our observed ratios and subtract 1.14 from them.

# Below code creates a list of ratios at the same index
# where the temp is 182 We then just subtract 1.868 from
# those values.
x <- df$ratio[which(df$temp == 182)]
x - 1.868

## [1] -0.968 -0.058  0.072  0.812

What we see is a fairly decent range of values. This is unsurprising as our observations still have randomness to them. Some of those ratios will fall above the line of best fit, others will fall below it. If it falls above the line the residual will be positive, if it’s below the residual will be negative.

d.

What proportion of the observed variation in efficiency ratio can be attributed to the simple linear regression relationship between the two variables? - This is not as intimidating as it sounds. What this problem wants is \(r^2\) which is the proportion of variability in y that is explained by the x value. We can calculate this easily. First we calculate the correlation, ‘r’, using cor() and then we square it.

cat("r^2 is approximately", round(stats::cor(temp, ratio)^2,
    3))

## r^2 is approximately 0.451

a.

Construct a stem-and-leaf display of the MOE values, and comment on any interesting features.

moe <- c(29.8, 33.2, 33.7, 35.3, 35.5, 36.1, 36.2, 36.3, 37.5,
    37.7, 38.7, 38.8, 39.6, 41, 42.8, 42.8, 43.5, 45.6, 46, 46.9,
    48, 49.3, 51.7, 62.6, 69.8, 79.5, 80)
stem(moe, scale = 2)

## 
##   The decimal point is 1 digit(s) to the right of the |
## 
##   2 | 
##   3 | 034
##   3 | 566668899
##   4 | 01334
##   4 | 66789
##   5 | 2
##   5 | 
##   6 | 3
##   6 | 
##   7 | 0
##   7 | 
##   8 | 00

Judging by the plot, it seems there is a huge amount of concentration of values between 35 and 40. There is an immediate drop off in frequency once 50 has been reached, with only 5 moe values at or above 50.

b.

Is the value of strength completely and uniquely determined by the value of MOE? Explain. Without looking at a plot I would say no. Looking at the data in the book, there is a wide spread in strength values whereas the moe values are pretty tightly grouped. I highly doubt that the value of strength is only determined by the value of moe.

strength <- c(5.9, 7.2, 7.3, 6.3, 8.1, 6.8, 7, 7.6, 6.8, 6.5,
    7, 6.3, 7.9, 9, 8.2, 8.7, 7.8, 9.7, 7.4, 7.7, 9.7, 7.8, 7.7,
    11.6, 11.3, 11.8, 10.7)
df <- data.frame(moe, strength)
ggplot(data = df, aes(x = moe, y = strength)) + geom_point() +
    theme_gray()

Judging from the scatterplot we can definitely see there is a chance of a relationship here. It’s still too scattered to say that strength is only being impacted by moe though.

###c. Use the accompanying minitab output (now shown here) to obtain: - The equation of the least squares line - prediction of strength given a moe of 40.

Would you feel comfortable using the least squares line to predict strength when moe is 100? Explain.

The minitab output would give us an equation of \[f(x) = 3.2925 + 0.10748x\]

Plugging 40 into our formula gives us \(f(40) \approx 7.5917\). As for if I would feel comfortable using this equation for predictions of strength given a moe of 100, I would not. We can fairly well predict values around 30 to 50 as we have many samples in that zone. We do not know for sure though if this relationship continues on like this as shown by the higher strength values. We only have 5 above 50 moe, and only 2 near 80. Predcictions in that region would already be unreliable, going beyond that would likely give us inaccurate predictions.

a.

Does a scatterplot of the data support the use of the simple linear regression model? - I would say yes. Based on the scatterplot there seems to be a very strong linear relationship here so it seems to be an appropriate model.

b.

Calculate the point estimates of the slope and intercept of the population regression line.

model <- lm(y ~ x, data = df)
model

## 
## Call:
## lm(formula = y ~ x, data = df)
## 
## Coefficients:
## (Intercept)            x  
##      -1.128        0.827

Based on the lm() functions output, our intercept is -1.128 and our slope is .827. This gives the equation \[f(x) = -1.128 + .827x\]

c.

Calculate a point estimate of the true average runoff volume when rainfall volume is 50.

Plugging 50 into our formula gives \(f(50) = -1.128 + .827 \cdot 50 \approx \boxed{40.22}\)

d.

Calculate a point estimate of the standard deviation \(\sigma\).

sse <- sum((stats::predict(model) - df$y)^2)
sigma <- sqrt(sse/(length(y) - 2))
cat("The sum of squares of residuals is", round(sse, 2), "\nThe point estimate of the standard deviation is",
    round(sigma, 2))

## The sum of squares of residuals is 357.01 
## The point estimate of the standard deviation is 5.24

e.

What proportion of the observed variation in runoff volume can be attributed to the simple linear regression relationship between runoff and rainfall?

cat("r^2 is approximately", round(cor(x, y)^2, 3))

## r^2 is approximately 0.975

a.

Obtain the equation of the estimated regression line. Then create a scatterplot of the data and graph the estimated line. Does it appear that the model relationship will explain a great deal of the observed variation in y?

model = lm(y ~ x, data = df)
model

## 
## Call:
## lm(formula = y ~ x, data = df)
## 
## Coefficients:
## (Intercept)            x  
##    118.9099      -0.9047

Based on the output of the function we get the following equation: \[f(x) = 118.9 - 0.9047x\]

ggplot(data = df, aes(x = x, y = y)) + geom_point() + theme_gray() +
    stat_smooth(method = lm, se = FALSE)

## `geom_smooth()` using formula 'y ~ x'

It appears that a lot, if not most, of the variation in y can be explained by x here.

b.

Interpret the slope of the least squares line.

The slope indicates that, if there is a relationship, that it’s negative. In other words, as the predictor grows the dependent variable decreases. In this case, it decreases by quite a lot!

c.

What happens if the estimated line is used to predict porosity when unit weight is 135? Why is this not a good idea.

new <- data.frame(x = c(135))
my_func <- function(x) {
    118.9 - 0.9048 * x
}
cat("The prediction for porosity when unit weight is 135 is",
    round(predict(model, newdata = new), 3), "\nNote: Above uses the predict function",
    "\n\nThe same prediction done manually gives", my_func(135))

## The prediction for porosity when unit weight is 135 is -3.229 
## Note: Above uses the predict function 
## 
## The same prediction done manually gives -3.248

Two methods for predicting are shown to be roughly equivalent here. The manual calculation utilized some rounding earlier on that explains the discrepancy. I did both methods just to sate my own curiosity.

As for why this prediction is a bad idea, the explanation is similar to an earlier question. We don’t have any x values that go that high, as such we don’t know if that behavior continues on indefinitely or if it shifts at some point.

d.

Calculate the residuals corresponding to the first two observations.

z <- df$y[1:2] - my_func(df$x[1:2])
cat("The first two residuals are", z[1], "and", z[2])

## The first two residuals are -0.5248 and 0.47528

e

Calculate and interpret a point estimate of \(\sigma\).

sse <- sum((stats::predict(model) - df$y)^2)
sigma <- sqrt(sse/(length(y) - 2))
cat("The sum of squares of residuals is", round(sse, 2), "\nThe point estimate of the standard deviation is",
    round(sigma, 2))

## The sum of squares of residuals is 11.44 
## The point estimate of the standard deviation is 0.94

The standard deviation is actually quite close to that of a standardized normal distribution. That’s pretty interesting.

f.

What proportion of observed variation in porosity can be attributed to the approximate linear relationship between unit weight and porosity?

cat("r^2 is approximately", round(cor(x, y)^2, 3))

## r^2 is approximately 0.974

a.

Fit the simple linear regression model to this data, predict colony density when surface area = 70 and 71, and calculate the corresponding residuals. How do they compare?

model

## 
## Call:
## lm(formula = y ~ x, data = df)
## 
## Coefficients:
## (Intercept)            x  
##    -305.881        9.963

my_generalized_func <- function(a, b, x) {
    a + b * x
}
prediction <- my_generalized_func(-305.881, 9.963, 70:71)
actual_value1 <- df$y[which(df$x == 70)]
actual_value2 <- df$y[which(df$x == 71)]
actual_value <- c(actual_value1, actual_value2)
resid <- c(actual_value[1:2] - prediction[1:2])

cat("The predicted colony density when surface area", "= 70 is",
    prediction[1], "\nThe predicted colony density when surface area",
    "= 71 is", prediction[2], "\nThe calculated residuals are",
    resid[1], "and", resid[2])

## The predicted colony density when surface area = 70 is 391.529 
## The predicted colony density when surface area = 71 is 401.492 
## The calculated residuals are -378.529 and 1527.508

The equation we get from this model is: \[f(x) = -305.881 + 9.963x\]

b.

Calculate the coefficient of determination.

cat("The r^2 value for this data is", round(cor(x, y)^2, 3))

## The r^2 value for this data is 0.124

c. 

The second observation has a very extreme y value. This observation may have a substantial impact on the fit of the model and subsequent calculations. Eliminate it and recalculate the equation of the estimated regression line. Does it appear to differ substantially from the equation before the deletion? What is the impact on \(r^2\) and \(s\)?

df <- df[-c(2), ]

model = lm(y ~ x, df)
model

## 
## Call:
## lm(formula = y ~ x, data = df)
## 
## Coefficients:
## (Intercept)            x  
##     34.3734       0.7792

The new equation we get from this is: \[f(x) = 34.3734 + 0.7792x\] This is, of course, very different from the previous equation.

cat("The r^2 value is", round(cor(df$x, df$y)^2, 3))

## The r^2 value is 0.024

The new \(r^2\) value is far less than previously! There doesn’t seem to be much of a linear relationship at all. For curiosity’s sake let us generate a new scatterplot with the new fitted line.

ggplot(data = df, aes(x = x, y = y)) + geom_point() + theme_gray() +
    stat_smooth(method = lm, se = FALSE)

## `geom_smooth()` using formula 'y ~ x'

As we can see from this plot, now that the outlier is removed the values seem far more scattered than it appeared when it was zoomed out.