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

<b>Letný semester 2024 | Cvičenie 8 | 16.04.2024</b><br>
  HTML Markdown -- <a href = "cvicenie8.Rmd">Rmd source file</a> (UTF encoding)
  
  
**************
**************
    
  
<b>Brief outline of the work for the lab session no.8:</b><br><br>
    <ul>
    <li>Homoscedastic vs. heteroscedastic data/model</li>
     <li>Some tools to deal with non-normally distributed data or heteroscedastic data</li>
    <li>Statistical inference in a model without normality (mean vs. variance structure)</li>
    </ul>
    
<br><br>
    
    
We will use a small dataset (47 observations and 8 different covariates) containing some pices of information about fires that occurred in Chicago in 1970 while distinguishing for some information provided in the data---particularly, we will focus on the locality in Chicago (north vs. south) and the proportion of minority citizens. The dataset is loaded from the website and some brief insight is provided below. 

    
 ```{R, warning = F, message = F}  
library("mffSM")

chicago <- read.csv("http://www.karlin.mff.cuni.cz/~maciak/NMFM334/chicago.csv", 
                    header=T, stringsAsFactors = TRUE)
        
chicago <- transform(chicago, fside = factor(side, levels = 0:1, 
                                      labels=c("North", "South")))          
                    
head(chicago)
summary(chicago)
```

The data contain information from Chicago insurance redlining in 47 districts in 1970 where
<i>minor</i> stands for the percentage of minority (0 - 100) in a given district, 
<i>fire</i> stands for the  number of fires per 100 households during the reference period, 
<i>theft</i> denotes the  number of reported thefts per 1000 inhabitants,
<i>old</i>  gives percentage of residential units built before 1939,
<i>insur</i> states the    number of policies per 100 household,
<i>income</i> is the  median income (USD 1000) per household,
and, finally,  <i>side</i>  gives the location of the district within Chicago (0 = North, 1 = South)


<br><br>


### <font color = "#880000">Indidividual work</font> 
How can we manually reconstruct the plots above? 
<ul>
<li>Use the data on fires in Chicago (the dataset loaded above) and perform some basic explanatory analysis while focusing on the number of fires (the dependent variable <i>fire</i> and two explanatory variables---the proportion of the minority (<i>minor</i>) and the location in Chicago---<i>side</i>). 
</li>

<li>What are you expectations about the form of the regression function when being interested in modeling of the conditional number of fires $E[Y | minor, side]$, while conditioning on the minority proportion and the locality information?</li>

<li>Are there some issues that you should be concern with? What about the asssumptions on normal regression model or an ordinary regression model (without the assumption of normality)?  </li>

<li>What particular graphical tools can be used to visualy verify the assumptions imposed on the model? </li>
</ul>
<hr size="3" width="100%" noshade style="color:#440000" align="left" />
<br><br>

### <b>1.Expected number of fires given minority and locality</b> 

We start by looking for a model of dependence of the number of fires (dependent variable <i>fire</i>)
on the percentage of minority (<i>minor</i>). We are also aware of the fact that 
the form of this dependence may differ in north and south parts of Chicago.



 ```{R, warning = F, message = F, fig.width = 10}   
XLAB <- "Minority (%)"
YLAB <- "Fires (#/100 households)"
COLS <- c(North = "darkgreen", South = "red")
BGS <- c(North = "aquamarine", South = "orange")
PCHS <- c(North = 21, South = 23)
plot(fire ~ minor, col = COLS[fside], bg = BGS[fside], pch = PCHS[fside], 
                   xlab = XLAB, ylab = YLAB, data = chicago)
legend(5, 40, legend = c("North", "South"), pch = PCHS, 
              col = COLS, pt.bg = BGS, bty = "n")
lines(with(subset(chicago, fside == "North"), lowess(minor, fire)), 
                                              col = COLS["North"], lwd = 2)
lines(with(subset(chicago, fside == "South"), lowess(minor, fire)), 
                                              col = COLS["South"], lwd = 2)

```

How would you interpret the plot? Do we need interactions in the model?
Does the side of Chicago modify the relationship between the number of fires
and the percentage of minority?

<br><br>

Lets compare the following three models. Based on the output describe the effect of the percentage of minority on the number of fires. Describe also the effect of the side.  

 ```{R, warning = F, message = F}   
m <- list()
m[["Line"]]  <- lm(fire ~ minor * fside, data = chicago)
m[["Log2"]]   <- lm(fire ~ log2(minor) * fside, data = chicago)
m[["Parab"]] <- lm(fire ~ (minor + I(minor^2)) * fside, data = chicago)

summary(m[["Line"]])
summary(m[["Log2"]])
summary(m[["Parab"]])
```

And the corresponding diagnostics: 
 ```{R, warning = F, message = F, fig.width = 10}   
plotLM(m[["Line"]])
plotLM(m[["Log2"]])
plotLM(m[["Parab"]])
```
 
<br><br>

For all three models there can be some issues spotted in the diagnostic plots above (heteroscedasticity, conditional mean specification, maybe normality issues). 
In addition, there is also a bunch of points scattered around the origin and relatively only a few observations are given for larger proportion of the minority or the larger amount of fires. This suggets that log transformation (of the response and the predictor as well) could help.  

<br><br>

### <b>2. Response (logarithmic) transformation</b> 

Typically, when the residual variance increases with the response
expectation, log-transformation of the response often
ensures a homoscedastic model (Assumption A3a from the lecture).

Firstly, compare the plot below (on the log scale on $y$ axes) with the previous plot on the original scale on the $y$ axes: 

 ```{R, warning = F, message = F, fig.width = 10}  
LYLAB <- "Log(Fires)"
par(mfrow = c(1, 1))
plot(log(fire) ~ minor, col = COLS[fside], bg = BGS[fside], pch = PCHS[fside], 
                        xlab = XLAB, ylab = LYLAB, data = chicago)
legend(5, 3.5, legend = c("North", "South"), pch = PCHS, 
               col = COLS, pt.bg = BGS, bty = "n")
lines(with(subset(chicago, fside == "North"), lowess(minor, log(fire))), 
                                              col = COLS["North"], lwd = 2)
lines(with(subset(chicago, fside == "South"), lowess(minor, log(fire))), 
                                              col = COLS["South"], lwd = 2)
```

Using the proposed transformation of the response we also change the underlying data. 
 This means, that instead of the model 
 $$
E[Y_i | \boldsymbol{X}_i] = \boldsymbol{X}_i^\top \boldsymbol{\beta} 
$$
fitted based on the random sample $\{(Y_i, \boldsymbol{X}_i);±i = 1, \dots, n\}$ 
we fit another modelof the form 
$$
E[\log(Y_i) | \boldsymbol{X}_i] = \boldsymbol{X}_i^\top \boldsymbol{\beta} 
$$
 which is, however, based on the different random sample $\{(\log(Y_i), \boldsymbol{X}_i);±i = 1, \dots, n\}$. This also has some crucial consequences. The logarithmic transformation will help to deal with some issues regarding heteroscedastic data but, on the other hand, it will introduce new problemes regarding the model interpretability. 
 
 The idea is that the final model should be (always) interpreted with respect to the original data---the information about the number of fires (<i>fire</i>)---not the logarithm of the number of fires (<i>log(fire)</i>). 

For some better interpretation we will also use the log transformation of the independent variable -- the information about the proportion of minority -- <i>minor</i>. What is the interpretation of the following model: 

 ```{R, warning = F, message = F, fig.width = 10}  
summary(m <- lm(log(fire) ~ log2(minor) * fside, data = chicago))
```

Recall the Jensen inequality and the fact that 
$$
E[\log(Y_i) | \boldsymbol{X}_i] \neq \log(E[Y_i | \boldsymbol{X}_i]).
$$
What are the consequences with respect to the interpretation of the parameter estimates provided in the summary output above? What would be a suitable interpretaion that we can use in this situation? 


<br><br>

It is easy to visualize the fitted regression function. In addition, we can also  add  the confidence bands AROUND the regression function and the prediction bands. However, given the previous problems related to the Jensen inequality, only the  prediction bands make sense for this model and the confidence (AROUND) bands does not have any good interpretation here.

 ```{R, warning = F, message = F, fig.width = 10}  
lpred <- 500
pminor <- seq(1, 100, length = lpred)
pdata <- data.frame(minor = rep(pminor, 2), 
                    fside = factor(rep(c("North", "South"), 
                    each = lpred)))
                    
cifit <-predict(m, newdata = pdata, interval = "confidence")
pfit <-predict(m, newdata = pdata, interval = "prediction")

YLIM <- range(log(chicago[, "fire"]), pfit)
XLIM <- range(chicago[, "minor"])

par(mfrow = c(1, 2))
for (ss in levels(chicago[, "fside"])){
  plot(log(fire) ~ minor, col = COLS[fside], bg = BGS[fside], pch = PCHS[fside],
       xlab = XLAB, ylab = LYLAB, xlim = XLIM, ylim = YLIM,
       data = subset(chicago, fside == ss), main = ss)
  ind <- if (ss == "North") 1:lpred else (lpred + 1):(2 * lpred)
  lines(pminor, cifit[ind, "fit"], col = COLS[ss], lwd = 3)
  lines(pminor, cifit[ind, "lwr"], col = COLS[ss], lwd = 2, lty = 2)
  lines(pminor, cifit[ind, "upr"], col = COLS[ss], lwd = 2, lty = 2)
  lines(pminor, pfit[ind, "lwr"], col = COLS[ss], lwd = 2, lty = 3)
  lines(pminor, pfit[ind, "upr"], col = COLS[ss], lwd = 2, lty = 3)
	legend(x=25, y=1, lwd=c(3,2,2), lty=c(1,2,3), col=COLS[ss], 
	             legend=c("estim. E[log(Y)|X=x]",
	                      "conf. band around E[log(Y)|X=x]", 
	                      "prediction band for log(Y)"))
}  
par(mfrow = c(1, 1))
```

What is  the interpretation of the estimated parameters and the specific numbers 
providsed in the summary output of the model with log-transformed dependent variable? 

<br><br>

### <b>3. Model improvement (mean vs. variance structure)</b>

Now, we will try to improve the model by adding the information about insurance policies (the independent variable <i>insur</i>). We will fit the following model: 

 ```{R, warning = F, message = F, fig.width = 10} 
m2 <- lm(log(fire) ~ (log(minor) + insur) * fside, data = chicago)
```

Note, that the fitted regression planes describing the dependence 
of log(fire) on log(minor) and insur based on the m2 model are the same 
as if the model is fitted separately in the North and South side of Chicago.

Hence if we are interested in just a description of those
planes---the conditional mean structure (and are lazy to extract the appropriate intercepts
and slopes from the m2 model), we can fit separate North/South models
and take the coefficients from there.

Compare the output above with the following two models: 
 ```{R, warning = F, message = F, fig.width = 10} 
m2North <- lm(log(fire) ~ log(minor) + insur, data = chicago, 
                                             subset = (fside == "North"))
m2South <- lm(log(fire) ~ log(minor) + insur, data = chicago, 
                                             subset = (fside == "South")) 
```

<br><br>

### <font color = "#880000">Indidividual work</font> 
Consider the model above and try to answer the following questions:  
<ul>
<li>Can you reconstruct the (mean) parameter estimates from the `m2` model obove using the two separate models for the north and south sides of Chicago? Can you do it vice versa? </li>
<li>Try some different parametrization of the locality. For instance 
 ```{R, warning = F, message = F, fig.width = 10} 
summary(lm(log(fire) ~ log(minor) * fside, data = chicago, 
                          contrasts = list(fside = contr.sum)))
                          
contr.sum(2)
```
What is the interpretation of the estimated parameters in the summary output above? 
</li>
</ul>
<hr size="3" width="100%" noshade style="color:#440000" align="left" />
<br><br>

<br><br>

