# <b>Základy regrese | NMFM 334</b>

<b>Letný semester 2024 | Cvičenie 7 | 09.04.2024</b><br>
  HTML Markdown -- <a href = "cvicenie7.Rmd">Rmd source file</a> (UTF encoding)
  
  
**************
**************
    
  
<b>Outline of the lab session no.6:</b><br><br>
    <ul>
    <li>Statistical inference in a normal linear model</li>
     <li>Confidence intervals, confidence regions, statistical tests</li>
     <li>Prediction interval</li>
    <li>Different bounds for and along the linear regression line</li>
    </ul>
    
<br><br>
    
    
We will consider a simple linear regression model (modelling the amount of yield 
given the information about the concentration of magnesium -- <i>Mg</i>). For some model 
improvements we will use the logarithm transformation of the independent covariate (<i>Mg</i>). 

    
 ```{R, warning = F, message = F}  
library("mffSM")
library("splines")       # B-splines modelling
library("MASS")          # function stdres for standardized residuals

### Working data 
data(Dris, package = "mffSM")
help("Dris")

head(Dris)
summary(Dris)

### a simple linear model fitted in the last session
aLogMg <- lm(yield ~ log2(Mg), data = Dris)

summary(aLogMg) 
```


We will use the model above and we will construct three different types of the confidence intervals/bounds for the model: 

<ol>
<li>Pointwise confidence intervals (region/band) around the regression line</li>

<li>Confidence region (band) FOR the (whole) regression line simultaneously</li>

<li>Prediction interval (for some new given value of $\boldsymbol{X}$)</li>
</ol>

<br><br>


### <b>1. Pointwise confidence intervals (region/band) around the regression line<</b> 

Firstly, use the help session in R to find out, what is the purpose of the predict command? Try the following: 


 ```{R, warning = F, message = F}   
help(predict)
help(predict.lm)
```

### <font color = "#880000">Indidividual work</font> 
How can we manually reconstruct the plots above? 
<ul>
<li>What is the difference between the outputs from "predict" and "fitted"? A few values 
are shown in the output below. 

 ```{R, warning = F, message = F}   
data.frame(Predict = predict(aLogMg), Fitted = fitted(aLogMg))[1:10,]
all.equal(predict(aLogMg), fitted(aLogMg))
```
</li>

<li>Now, we will use the same command, however, for some new value of $X$ (concentration of magnesium)
 ```{R, warning = F, message = F}   
predict(aLogMg, newdata = data.frame(Mg = 10))

(point.fit<-aLogMg$coef[1]+aLogMg$coef[2]*log2(10))
```
</li>

<li>It can be even obtained with the corresponding standard error
 ```{R, warning = F, message = F}   
predict(aLogMg, newdata = data.frame(Mg = 10), se.fit = TRUE)
```
The standard error can be also reconstructed manually
 ```{R, warning = F, message = F}   
lvec <- c(1, log2(10))
(se.fit <- c(sqrt(t(lvec) %*% vcov(aLogMg) %*% lvec)))
all.equal(predict(aLogMg, newdata = data.frame(Mg = 10), se.fit = TRUE)$se.fit,
          se.fit)
          
# Recall - residual degrees of freedom, i.e. $df
aLogMg$df.residual
# sqrt(MS_e) = estimate of sigma, i.e. $residual scale
sqrt(deviance(aLogMg)/aLogMg$df.residual)
summary(aLogMg)$sigma
```
</li>

<li>Alternativelly, the output can be also obtained from the `LSest()` function in the R library <i>mffSM</i> 

 ```{R, warning = F, message = F}   
LSest(aLogMg, L = c(1, log2(10)))
```
</li>
</ul>
<hr size="3" width="100%" noshade style="color:#440000" align="left" />
<br><br>




Whan we are interested in the estimate of the conditional mean $E[Y | \boldsymbol{X} = \boldsymbol{x}]$ we can either provide the point estimate or the interval estimate instead. 
For the interval estimate, we need the confidence level and the estimate of the variance parameter. 

It can be given automatically as follows; 
 ```{R, warning = F, message = F}   
predict(aLogMg, newdata = data.frame(Mg = 10), interval = "confidence", 
        level = 0.95)
```
 
Or, alternatively, we can reconstruct everything manually: 
 ```{R, warning = F, message = F}   

### Manual computation of the confidence interval 
crossprod(c(1, log2(10)), coef(aLogMg))   ## What is this?
as.numeric(crossprod(c(1, log2(10)), coef(aLogMg))) + 
  c(-1, 1) * se.fit * qt(0.975, df = aLogMg[["df.residual"]])

### We can calculate a sequence of confidence intervals for E(Y|x) 
### for different values of x
predict(aLogMg, newdata = data.frame(Mg = 10:20), interval = "confidence", 
        level = 0.95)
predict(aLogMg, newdata = data.frame(Mg = 10:20), interval = "confidence", 
        level = 0.95, se.fit = TRUE)

### A denser grid: preparation of a grid of points for the pointwise confidence 
### band around the regression line
x.Mg <- seq(min(Dris[, "Mg"]) - 1, max(Dris[, "Mg"]) + 1, length = 101)

### Calculate the confidence intervals
ciEY.aLogMg <- predict(aLogMg, data.frame(Mg = x.Mg), interval = "confidence", 
                       level = 0.95)
```


Observations together with the estimated regression line and a sequence of confidence intervals for $E[Y|\boldsymbol{X} = \boldsymbol{x}]$ are visualized as a band -- confidence region (band) AROUND the regression line. 

 ```{R, warning = F, message = F, fig.width = 10}   
plot(yield ~ log2(Mg), pch = 23, col = "orange", bg = "beige", 
     xlab = "Log2(Mg)", ylab = "Yield", data = Dris)
lines(log2(x.Mg), ciEY.aLogMg[, "fit"], lwd = 2, col = "red3")
lines(log2(x.Mg), ciEY.aLogMg[, "lwr"], lwd = 2, col = "red3", lty = 2)
lines(log2(x.Mg), ciEY.aLogMg[, "upr"], lwd = 2, col = "red3", lty = 2)

### Where is the band narrowest? 
lines(rep(mean(log2(Dris$Mg)), 2), c(min(Dris$yield), max(Dris$yield)))
lines(c(min(log2(Dris$Mg)), max(log2(Dris$Mg))), rep(mean(Dris$yield), 2))
```

Everything can be also plotted with respect to the original scale of the magnesium concentration: 

 ```{R, warning = F, message = F, fig.width = 10}  
plot(yield ~ Mg, pch = 23, col = "orange", bg = "beige", xlab = "Mg", 
     ylab = "Yield", data = Dris)
lines(x.Mg, ciEY.aLogMg[, "fit"], lwd = 2, col = "red3")
lines(x.Mg, ciEY.aLogMg[, "lwr"], lwd = 2, col = "red3", lty = 2)
lines(x.Mg, ciEY.aLogMg[, "upr"], lwd = 2, col = "red3", lty = 2)

### What about now?
lines(rep(mean(Dris$Mg), 2), c(min(Dris$yield), max(Dris$yield)))
lines(c(min(Dris$Mg), max(Dris$Mg)), rep(mean(Dris$yield), 2))
```

<br><br>



### <b>2. Confidence region (band) FOR the (whole) regression line simultaneously</b> 

The second option provides a confidence band for the whole regression line simultaneously. 
Note, that there is diffferent interpretation of the point-wise band constructed above (where each 
interval can be interpreted only individually, for some given concentration of magnesium -- generally, for some given value of the explanatory vector $\boldsymbol{X} = \boldsymbol{x}$). <br><br>

Some theoretical details towards the band FOR the whole regression line simultaneously were given during the lecture, nevertheless, the code provided below can be considered as 
purely illustrative. It is only important to understand the difference 
 between point-wise band "AROUND" the regression line constructed in 1). and the simultaneous band "FOR" the whole regression line at once.

 ```{R, warning = F, message = F, fig.width = 10}  
ciEY2.aLogMg <- predict(aLogMg, data.frame(Mg = x.Mg), se.fit = TRUE)
forEY.aLogMg <- data.frame(fit = ciEY2.aLogMg[["fit"]],
          lwr = ciEY2.aLogMg[["fit"]] - 
                   ciEY2.aLogMg[["se.fit"]] * 
                     sqrt(qf(0.95, 2, ciEY2.aLogMg[["df"]])*2),
          upr = ciEY2.aLogMg[["fit"]] + 
                   ciEY2.aLogMg[["se.fit"]] * 
                     sqrt(qf(0.95, 2, ciEY2.aLogMg[["df"]])*2))
head(forEY.aLogMg)
```

The computation of the confidence band is based on the confidence ellipse 
for the regression parameters (beta_0, beta_1)

 ```{R, warning = F, message = F, fig.width = 10}  
confidenceEllipse(aLogMg)
```

It can be visualized on the logarithm scale


 ```{R, warning = F, message = F, fig.width = 10}  
plot(yield ~ log2(Mg), pch = 23, col = "orange", bg = "beige", xlab = "Log(Mg)",
     ylab = "Yield", data = Dris)
lines(log2(x.Mg), ciEY.aLogMg[, "fit"], lwd = 2, col = "darkgreen")
### Band AROUND
lines(log2(x.Mg), ciEY.aLogMg[, "lwr"], lwd = 2, col = "red3", lty = 2)
lines(log2(x.Mg), ciEY.aLogMg[, "upr"], lwd = 2, col = "red3", lty = 2)
### Band FOR
lines(log2(x.Mg), forEY.aLogMg[, "lwr"], lwd = 2, col = "blue2", lty = 4)
lines(log2(x.Mg), forEY.aLogMg[, "upr"], lwd = 2, col = "blue2", lty = 4)
legend(4, 8.5, legend = c("For", "Around"), col = c("blue2", "red3"), 
       lty = c(4, 2), lwd = 2)
```

or the original scale of the logarithm concentration: 

 ```{R, warning = F, message = F, fig.width = 10}  
plot(yield ~ Mg, pch = 23, col = "orange", bg = "beige", xlab = "Mg", 
     ylab = "Yield", data = Dris)
lines(x.Mg, ciEY.aLogMg[, "fit"], lwd = 2, col = "darkgreen")
### Band AROUND
lines(x.Mg, ciEY.aLogMg[, "lwr"], lwd = 2, col = "red3", lty = 2)
lines(x.Mg, ciEY.aLogMg[, "upr"], lwd = 2, col = "red3", lty = 2)
### Band FOR
lines(x.Mg, forEY.aLogMg[, "lwr"], lwd = 2, col = "blue2", lty = 4)
lines(x.Mg, forEY.aLogMg[, "upr"], lwd = 2, col = "blue2", lty = 4)
legend(18, 8.5, legend = c("For", "Around"), col = c("blue2", "red3"), 
       lty = c(4, 2), lwd = 2)
```


### <b>3. Prediction interval (for some new given value of $\boldsymbol{X}$)</b> 

Finally, we will compare both confidence bands constructed above with one more confidence band -- the prediction band (which will be again interpreted in a point-wise manner as the first confidence band around the regression line). <br><br>

The idea is to construct an interval which covers---with probability $1 - \alpha$ 
the NEW value of the response $Y$ obtained for $\boldsymbol{X}$ being set to $\boldsymbol{x}$. 

Compare the following two commands: 


 ```{R, warning = F, message = F, fig.width = 10}  
predict(aLogMg, newdata = data.frame(Mg = 10), interval = "prediction")    
predict(aLogMg, newdata = data.frame(Mg = 10), interval = "confidence") 
```

Manual calculation of the prediction interval can follow the lines below: 

 ```{R, warning = F, message = F, fig.width = 10}  
s = summary(aLogMg)$sigma
se.fit2 = sqrt(s^2+(se.fit)^2)
as.numeric(crossprod(c(1, log2(10)), coef(aLogMg))) + 
  c(-1, 1) * se.fit2 * qt(0.975, df = aLogMg[["df.residual"]])
  
### Prediction band - a sequence of prediction intervals for a dense grid of 
### x-values
pred.aLogMg <- predict(aLogMg, data.frame(Mg = x.Mg), interval = "prediction", 
                       level = 0.95)
head(pred.aLogMg)
```

Finally, we can visualize the result -- using the logarithm scale


 ```{R, warning = F, message = F, fig.width = 10}  
plot(yield ~ log2(Mg), pch = 23, col = "orange", bg = "beige", xlab = "Log(Mg)",
     ylab = "Yield", data = Dris)
lines(log2(x.Mg), ciEY.aLogMg[, "fit"], lwd = 2, col = "red3")
### Band AROUND the regression line
lines(log2(x.Mg), ciEY.aLogMg[, "lwr"], lwd = 2, col = "red3", lty = 2)
lines(log2(x.Mg), ciEY.aLogMg[, "upr"], lwd = 2, col = "red3", lty = 2)
### Prediction band
lines(log2(x.Mg), pred.aLogMg[, "lwr"], lwd = 2, col = "blue2", lty = 4)
lines(log2(x.Mg), pred.aLogMg[, "upr"], lwd = 2, col = "blue2", lty = 4)
legend("topleft", legend = c("Prediction", "Around the regression line"), 
       col = c("blue2", "red3"), lty = c(4, 2), lwd = 2)
```

or using the original scale of the magnesium concentration: 

 ```{R, warning = F, message = F, fig.width = 10}  
plot(yield ~ Mg, pch = 23, col = "orange", bg = "beige", xlab = "Mg", 
     ylab = "Yield", data = Dris)
lines(x.Mg, ciEY.aLogMg[, "fit"], lwd = 2, col = "red3")
### Band AROUND the regression line
lines(x.Mg, ciEY.aLogMg[, "lwr"], lwd = 2, col = "red3", lty = 2)
lines(x.Mg, ciEY.aLogMg[, "upr"], lwd = 2, col = "red3", lty = 2)
### Prediction band
lines(x.Mg, pred.aLogMg[, "lwr"], lwd = 2, col = "blue2", lty = 4)
lines(x.Mg, pred.aLogMg[, "upr"], lwd = 2, col = "blue2", lty = 4)
legend("bottomright", legend = c("Prediction", "Around the regression line"), 
       col = c("blue2", "red3"), lty = c(4, 2), lwd = 2)
```
<br><br>


### <font color = "#880000">Indidividual work</font> 
How can we manually reconstruct the plots above? 
<ul>
<li>Recall the theoretical foundations of all three confidence bands constructed above. 
What are the main principial differences? 
</li>

<li>What are particular utilization purposees for each of the confidence bands constructed above? What type inferential tool do they offer? </li>
</ul>
<hr size="3" width="100%" noshade style="color:#440000" align="left" />
<br><br>

<br><br>

