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

<b>Letný semester 2024 | Cvičenie 6 | 02.04.2024</b><br>
  HTML Markdown -- <a href = "cvicenie6.Rmd">Rmd source file</a> (UTF encoding)
  
  
**************
**************
    
  
<b>Outline of the lab session no.5:</b><br><br>
    <ul>
    <li>Model assumptions and diagnostic tools</li>
     <li>Graphical diagnostic tools</li>
    </ul>
    
<br><br>
    
    
Firstly, some necessary R libraries and the working dataset: 

 ```{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)
```

<br><br>

### <b>1. Basic model diagnostics</b> 

We will use the following model at the beginning: 

 ```{R, warning = F, message = F}   
aLogMg <- lm(yield ~ log2(Mg), data = Dris)

summary(aLogMg) 
```

We will use the model above and we will closely inspect the model residuals. 
Three specific plots will be of particular interest: 

<ul>
<li><b>Raw residuals</b> $u_1, \dots, u_n$, that are supposed to satisfy: 
<ul>
<li>$E[u_i | \boldsymbol{X}_i] = 0$</li>
<li>$Var(u_i|\boldsymbol{X}_i) = \sigma^2 m_{ii}$, where $m_ii$ is a diagonal element 
of the matrix $\mathbb{M} = (\mathbb{i} - \mathbb{H})$</li>
<li>in a normal linear model, it also holds that $u_i$ follow the normal distribution, however, 
with different variance for different $i \in \{1, \dots, n\}$</li>
</ul>
</li>

<li><b>Standardized residuals</b> $v_i = u_i / \sqrt{MSe \times m_Pii{}}$, for $i = 1, \dots, n$
that satify
<ul>
<li>$E[v_i | \boldsymbol{X}_i] = 0$</li>
<li>$Var(v_i | \boldsymbol{X}_i) = 1$</li>
<li>however, the standardized residuals are not normal -- not even under the normal linear model</li>
</ul>

</li>
<ul>


<br><br>

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

<br><br>

### <font color = "#880000">Indidividual work</font> 
How can we manually reconstruct the plots above? 
<ul>
<li>First plot: 
 ```{R, warning = F, message = F, fig.width = 10}  
plot(aLogMg, which = 1)

### manual computation
points(fitted(aLogMg),residuals(aLogMg),lwd=2,col="blue4",bg="skyblue",pch=21)
lines(lowess(fitted(aLogMg),residuals(aLogMg)),col=2,lwd=3)
```
</li>

<li>Second plot: 
 ```{R, warning = F, message = F, fig.width = 10}  
plot(aLogMg, which = 2, ask = FALSE)

### standardized residuals
H = qr.Q(aLogMg[["qr"]]) %*% t(qr.Q(aLogMg[["qr"]]))
M = diag(nrow(H))-H
stdres = residuals(aLogMg)/sqrt(deviance(aLogMg)/aLogMg$df.residual*diag(M))   
qq = qqnorm(stdres,plot.it=FALSE)
points(qq$x, qq$y,lwd=2,col="blue4",bg="skyblue",pch=21)
all.equal(qq$y,stdres)
qqline(stdres,lwd=2,lty=3)
```
</li>

<li>Third plot
 ```{R, warning = F, message = F, fig.width = 10}  
plot(aLogMg, which = 3, ask = FALSE)  

### manual computation
points(fitted(aLogMg),sqrt(abs(stdres)),lwd=2,col="blue4",bg="skyblue",pch=21)
lines(lowess(fitted(aLogMg),sqrt(abs(stdres))),col=2,lwd=3)
```
</li>
</ul>
<hr size="3" width="100%" noshade style="color:#440000" align="left" />
<br><br>


Note,  that there are actually infinitely many possibilities how to check that 
$E(u_i| \boldsymbol{X}_i) = 0$ (a scatterplot of fitted values vs. residuals is only one 
of them, other common possibility is to plot the covariate values vs. the residuals. 

 ```{R, warning = F, message = F, fig.width = 10}  
plot(log2(Dris$Mg), residuals(aLogMg), xlab=expression(log[2](Mg)), 
     ylab="Residuals", main="Residuals vs. covariate", pch = 21, col = "blue4", 
     bg = "skyblue");
lines(lowess(residuals(aLogMg)~log2(Dris$Mg)), col="red")
```
Also check the folowing and explain the results: 

 ```{R, warning = F, message = F, fig.width = 10} 
sum(aLogMg$residuals)
sum(aLogMg$residuals*log2(Dris$Mg))
```

<br><br>


### <b>2. Polynomials and splines as regressors</b> 

For some example of possible model improvements, we will consider a different dataset 
now and we will fit a series of different model that we will mutually compare. 


A dataset providing a series of (indepeendent) measurements of head acceleration
in a simulated motorcycle accident (used to test crash helmets) is given in the datafile 
<i>Motocycle</i> (within the R packages <i>mffSM</i>). 

 ```{R, warning = F, message = F, fig.width = 10} 
data(Motorcycle, package = "mffSM")
help(Motorcycle)

head(Motorcycle)
dim(Motorcycle)
summary(Motorcycle)
```

We are primarily interested in how the head acceleration (covariate `haccel`) depends on the time 
after the accident (given in miliseconds -- the timecovaraite).


 ```{R, warning = F, message = F, fig.width = 10} 
BTY <- "o";
par(mfrow = c(1, 1), bty = BTY, mar = c(4, 4, 1, 1) + 0.1)
plot(haccel ~ time, data = Motorcycle, pch = 21, bg = "orange", col = "red3",
     xlab = "Time [ms]", ylab = "Head acceleration [g]")

fit0 = lm(haccel~time, data=Motorcycle)
abline(fit0)
```

Using the diagnostic plots, we can  clearly see that the model denoted as <i>fit0</i> 
is hardly appropriate: 

 ```{R, warning = F, message = F, fig.width = 10} 
plotLM(fit0)
```

As far as there is a quite flexible shape observed in the plot we need some more 
flexibile model to capture the structure of the underlying conditional expectation of the 
acceleraton given the time. We will use polynomials for fiting the model 
(while still staying within the linear regression framework). 

 ```{R, warning = F, message = F, fig.width = 10} 
### Fitting polynomial of various degrees
fit6 <- lm(haccel~poly(time, 6, raw=TRUE), data=Motorcycle)  
fit12 <- lm(haccel~poly(time, 12, raw=TRUE), data=Motorcycle)
```

The model can be visualized as follows: 

 ```{R, warning = F, message = F, fig.width = 10} 
xgrid <- seq(0, 60, length = 500)

par(mfrow = c(1, 1), bty = BTY, mar = c(4, 4, 1, 1) + 0.1)
plot(haccel ~ time, data = Motorcycle, pch = 21, bg = "orange", col = "red3",
     xlab = "Time [ms]", ylab = "Head acceleration [g]")
lines(xgrid, predict(fit6, newdata=data.frame(time=xgrid)), lwd=2, col="blue")
lines(xgrid, predict(fit12, newdata=data.frame(time=xgrid)), lwd=2, 
      col="darkgreen")
legend(37,-70, lty=rep(1, 2), col=c("blue", "darkgreen"), 
       legend=paste("p=", c(6,12),sep=""), lwd=c(2,2), cex=1.2)
```


Classical polynomials are, however, not computationally very effective. Especially for 
higher order of the polynomial approximation the solution tends to be unstable. 

One possible solution how to introduce some reasonable stability into the model 
where some higher order polynomials are used, are <b>splines</b>.

 ```{R, warning = F, message = F, fig.width = 10}
### Choice of knots 
knots <- c(0, 11, 12, 13, 20, 30, 32, 34, 40, 50, 60) 
(inner <- knots[-c(1, length(knots))])  # inner knots 
(bound <- knots[c(1, length(knots))])   # boundary knots 

### B-splines of degree 3
DEGREE <- 3; # piecewise cubic 
Bgrid <- bs(xgrid, knots = inner, Boundary.knots = bound, degree = DEGREE, 
            intercept = TRUE)
dim(Bgrid)   # 13-dimensional linear space of piecewise polynomials

COL <- rainbow_hcl(ncol(Bgrid), c = 80)
par(mfrow = c(1, 1), bty = BTY, mar = c(4, 4, 1, 1) + 0.1)
plot(xgrid, Bgrid[,1], type = "l", col = COL[1], lwd = 2, xlab = "x", 
     ylab = "B(x)", xlim = range(knots), ylim = range(Bgrid, na.rm = TRUE))
for (j in 2:ncol(Bgrid)) lines(xgrid, Bgrid[,j], col = COL[j], lwd = 2)
points(knots, rep(0, length(knots)), pch = 21, bg = "blue", col = "skyblue", 
       cex = 1.5)
```

Thus, we can used B-splines to fit the model. <br>
More details about the B-splines can be found, for instance, here 
<a href = "http://www.chebfun.org/examples/approx/BSplineConv.html">http://www.chebfun.org/examples/approx/BSplineConv.html</a><br><br>

 ```{R, warning = F, message = F, fig.width = 10}
m1 <- lm(haccel ~ bs(time, knots = inner, Boundary.knots = bound, 
                     degree = DEGREE, intercept = TRUE) - 1, data = Motorcycle)
names(m1$coefficients) <- paste("B", 1:13)
summary(m1)
```

Note, that as the intercept is not formally included in in the model (because of
the -1 in the model formula), neither 'Multiple R-squared'
nor the 'F-statistic' quantity in the output are meaningful numbers. 

<br><br>
### <font color = "#880000">Indidividual work</font> 
Can you explain why the values reported for $R^2$ and $R_{adj}^2$ are not 
appropriate? What is the right meaning of this numbers? Can you manually reconstruct them? 

<ul>
<li> The $R^2$ value from the model is
 ```{R, warning = F, message = F, fig.width = 10}
summary(m1)$r.squared
```
and the manual calculations will give
 ```{R, warning = F, message = F, fig.width = 10}
1 - deviance(m1)/sum(Motorcycle$haccel^2)
sum(fitted(m1)^2)/sum(Motorcycle$haccel^2)
```
</li>
<li>Interpretation the F-test from the summary output above is a test of a submodel
 ```{R, warning = F, message = F, fig.width = 10}
summary(m1)$fstatistic

# the usual F-test is the following
m2 <- lm(haccel ~ 1, data = Motorcycle);
anova(m2,m1)

# in a model without intercept it is this one
(m0 <- lm(haccel ~ -1, data=Motorcycle))
anova(m0,m1)
```
</li>

<li>Nevertheless, as the intercept paramter is in the linear span of the columns of \mathbb{X}$,
it makes sense to calculcate the (usual) multiple $R^2$ value.
This can do it 'manually' as follows: 
 ```{R, warning = F, message = F, fig.width = 10}
var(fitted(m1))/var(Motorcycle$haccel)
cor(Motorcycle$haccel, fitted(m1))^2
1 - deviance(m1)/deviance(m2)
```
</li> 

<li>Alternatively one can use B-splines (the basis without intercept) 
with an intercept in the model to get some meaningful numbers in the output 
 ```{R, warning = F, message = F, fig.width = 10}
m0 <- lm(haccel ~ bs(time, knots = inner, Boundary.knots = bound, 
                     degree = DEGREE, intercept = FALSE), data = Motorcycle)
names(m0$coefficients) <- c("(Intercept)", paste("B", 1:12))
summary(m0)

anova(m2,m0)
```



</li>

</ul>



What about the interpretation of the classical diagnostic plots 
for the model fitted with B-slines? 
 ```{R, warning = F, message = F, fig.width = 10}
plotLM(m1)
```


<br><br>