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

<b>Letný semester 2024 | Cvičenie 5 | 26.03.2024</b><br>
HTML Markdown -- <a href = "cvicenie5.Rmd">Rmd source file</a> (UTF encoding)


**************
**************

<b>Outline of the lab session no.5:</b><br><br>
<ul>
<li>Multiple regression model with interactions</li> 
<li>Transformations of the explanatory variables</li> 
<li>Hierarchically well formulated model vs. non-hierarchical model</li>
</ul>

<br><br>


For the purpose of this lab session we will consider the dataset `Dris` from the R library <i>mffSM</i>. 
We will be  interested  in the dependence of the yield on the amount of  magnesium concentration (yield ~ Mg) 
and the amount of the calcium concentration (yield ~ Ca). 

All together, there are 368 observations avaiable in the dataset and 6 different covariates (although, we will not use all of them). A brief insight about the data can be taken from the following code: 



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

data(Dris, package = "mffSM")

dim(Dris)
head(Dris)
summary(Dris)
```

Firstly, we visualize marginal dependence of yield on th magnesium concentration and 
the calcium concentration. 

```{R, warning = F, message = F, fig.width = 10}
par(mfrow = c(1,2))

plot(yield ~ Mg, data = Dris, pch = 22, bg = "gray", xlab = "Magnesium concentration", ylab = "Yield")
plot(yield ~ Ca, data = Dris, pch = 22, bg = "gray", xlab = "Calcium concentration", ylab = "Yield")
```

Note, that the points are not equaly distributed along the $x$ axis. There are way more points at the left corners of the scatterplots than in the right parts. Ths suggest, that a logarithmic transformation (of the magnesium and calcium concentration) could help to align the data more uniformly along the $x$ axis,

For more intuitive interpretation we use the logarithm with the base two (meaning that unit increase of the logaritmically transformed covariates, i.e.,  $\log_2(x) + 1 = \log_2 (x) \cdot \log_2 2  = \log_2 (2 x)$ equals double amount of the original covariate). 

```{R, warning = F, message = F, fig.width = 10}
### Add transformations of Mg, Ca, N to the dataset 
Dris <- transform(Dris, lMg = log2(Mg), lMg2 = (log2(Mg))^2, lCa = log2(Ca), 
                  lCa2 = (log2(Ca))^2, lN = log2(N), lN2=(log2(N))^2)

par(mfrow = c(1,2))
plot(yield ~ lMg, data = Dris, pch = 22, bg = "gray", xlab = "Magnesium concentration [log scale]", ylab = "Yield")
plot(yield ~ lCa, data = Dris, pch = 22, bg = "gray", xlab = "Calcium concentration [log scale]", ylab = "Yield")
```

<br><br>

For the following, we will use the tranformed data. <br><br>

### <b>1. Simple and multiple linear regression models</b> 

Let us start with two simple regression model (one explanatory variable only) and we compare these 
models with a multiple regression model (with both explanatory variables being used at the same time).

```{R}
summary(m10 <- lm(yield ~ lMg, data = Dris))
summary(m01 <- lm(yield ~ lCa, data = Dris))
```

Both marginal model can be now compared with a larger model, containing both predictor variables: 

```{R}
### Additive model
summary(m11 <- lm(yield ~ lMg + lCa, data = Dris))
```

Note, there quite different values for the estimated paramters when being interested in the effect of magnesium or the effect of calcium. Which model should be taken as a reference one? 

<br><br>

### <font color = "#880000">Indidividual work</font> 
Interpret the models fitted above -- try to explain the meaning of the estimated parameters. 
In particular, focus on the following: 
<ul>
<li>Interpret the estimated regression coefficients in the models.</li>
<li>What is being  tested on the rows started with 'lMg' and 'lCa'?</li>
<li>When comparing all three models, which one is, according your our opinion, the most suitable one?</li>
</ul>
<hr size="3" width="100%" noshade style="color:#440000" align="left" />
<br><br>


In the following, we will try slightly more complex (multiple) regression models:
```{R}
summary(m20 <- lm(yield ~ lMg + lMg2, data = Dris))
summary(m02 <- lm(yield ~ lCa + lCa2, data = Dris))
```

Is there some statistical improvements in the two models above? 

And what about the following models? 

```{R}
summary(m21 <- lm(yield ~ lMg + lMg2 + lCa, data = Dris)) 
summary(m22 <- lm(yield ~ lMg + lMg2 + lCa + lCa2, data = Dris)) 
```

What is now the effect of the magnesium concentration on the expected (estimated) amount of yield? 
How this effect can be quantified using different models? <br><br>

We can actually visualize the effect of the magnesium concentration on the yeld (using different models fitted above). Specifically, we are interested in the effect of the change of lMg by additional log2(1.1) which corresponds with the increase of Mg by 10%). What will be the corresponding change of the expected abount of yield, i.e. $E[yield | lMg, lCa]$?

```{R, fig.width = 10, warning = F}
eps <- 0.1;
delta <- log2(1+eps);
lMg.grid <- seq(min(Dris$lMg), max(Dris$lMg)-delta, length = 1000); 
effect.lMg <- function(b1,b2) b1*delta + b2*(2*delta*lMg.grid + delta^2); 

par(mfrow = c(1,2))

plot(lMg.grid, effect.lMg(coef(m20)['lMg'], coef(m20)['lMg2']), col = "blue", 
               xlab = "log2(Mg)", ylab = "Change of expected yield", 
               main = "The effect of 10%-increase of Mg on the yield", 
               type="l", lwd=2, ylim=range(effect.lMg(coef(m20)['lMg'], 
               coef(m20)['lMg2']), effect.lMg(coef(m21)['lMg'], 
               coef(m21)['lMg2']), effect.lMg(coef(m22)['lMg'], 
               coef(m22)['lMg2'])));
lines(lMg.grid, effect.lMg(coef(m21)['lMg'], coef(m21)['lMg2']), 
                lwd=2, col="red")
lines(lMg.grid, effect.lMg(coef(m22)['lMg'], coef(m22)['lMg2']), 
                lwd=2, col="yellow", lty="dotted") 
abline(h=0, lty="dotted")
legend(3.25, -0.1, col=c("blue", "red", "yellow"), lwd=c(2,2,2), 
                lty=c("solid", "solid", "dotted"), 
                legend=c("m20", "m21", "m22"), cex=1.3)


Mg.grid <- seq(min(Dris$Mg), max(Dris$Mg)-1, length = 1000); 
effect.Mg <- function(b1,b2){
  return(b1*(log2(Mg.grid+1) - log2(Mg.grid)) + 
           b2*((log2(1+Mg.grid))^2 - (log2(Mg.grid))^2))
}

plot(Mg.grid, effect.Mg(coef(m20)['lMg'], coef(m20)['lMg2']), col = "blue",
              xlab = "Mg", ylab = "Change of expected yield",
              main = "The effect of 10%-increase of Mg on the yield",
              type="l", lwd=2, ylim=range(effect.Mg(coef(m20)['lMg'], 
              coef(m20)['lMg2']), effect.Mg(coef(m21)['lMg'], 
              coef(m21)['lMg2']), effect.Mg(coef(m22)['lMg'], 
              coef(m22)['lMg2'])));
lines(Mg.grid, effect.Mg(coef(m21)['lMg'], coef(m21)['lMg2']), 
              lwd=2, col="red")
lines(Mg.grid, effect.Mg(coef(m22)['lMg'], coef(m22)['lMg2']), 
              lwd=2, col="yellow", lty="dotted")
abline(h=0, lty="dotted")
legend(15, 0.3, col=c("blue", "red", "yellow"), lwd=c(2,2,2), 
            lty=c("solid", "solid", "dotted"), legend=c("m20", "m21", "m22"), 
            cex=1.3)
rug(Dris$Mg+runif(nrow(Dris))-0.5);
```


Considering the models above, how would you (statistically) answer a question whether that given the information about the 
magnesium ('Mg') the amount of yield is independent of the concentration of the calcium?

How would you test that given Mg (with the quadratic dependence) the amount of  'yield' is independent of Ca?

Can we say that the amount of yield is independent of Ca?  Can we say that given Mg and Ca the amount of yield is independent 
of N?


<br><br>

### <b>2. Regression models with simple interations</b> 
The interactions in the model allow for modelling a flexible effect of some covariate given the value of some other covariate (or multiple covariates). 
We consider a simple model with one interaction term (denoted as ``lMg:lCa):

```{R}
summary(m1_inter <- lm(yield ~ lMg + lCa + lMg:lCa, data = Dris))
```


What is the interpretation of the estimated parameters now? What is the effect of the magnesium concentration and the calcium concentration on the amount of yield? Compare the interaction model above with the mode below (multiple model with three explanatory variables, where the last variable $Z$ is created artificiall from the magnesium concentration and the calcium concentration).

```{R}
Dris$Z <- Dris$lMg * Dris$lCa
summary(m1_Z <- lm(yield ~ lMg + lCa + Z, data = Dris))
```

Again, we can quantify the effect of the magnesium concentration on the yield however, this time 
we have to properly take into account the underlying concentration of calcium. How can be this done? 


```{R, fig.width = 10}
Mg.grid <- seq(min(Dris$Mg), max(Dris$Mg), length = 1000)
Ca.grid <- seq(min(Dris$Ca), max(Dris$Ca), length = 1000)
lCa.grid <- seq(min(Dris$lCa), max(Dris$lCa), length = 1000)
eps <- 0.1;
delta <- log2(1+eps)


par(mfrow = c(1, 2))
plot(lCa.grid, delta*coef(m1_inter)['lMg'] + 
               delta*coef(m1_inter)['lMg:lCa'] * lCa.grid, col = "blue", 
         xlab = "log2(Ca)", ylab = "Change of expected yield", 
         main = "The effect of 10%-increase of Mg on the yield", type="l", lwd=2)

plot(Ca.grid, coef(m1_inter)['lMg'] * (log2(Mg.grid+1)-log2(Mg.grid)) + 
              coef(m1_inter)['lMg:lCa'] * (log2(Mg.grid+1)-
              log2(Mg.grid)) * log2(Ca.grid), col = "blue", 
              xlab = "Ca", ylab = "Change of expected yield", 
              main = "The effect of increase of Mg by +1 on the yield", 
              type="l", lwd=2)
```


<br><br>

### <b>3. Interations with categorical covariates</b> 

Now, we will consider another dataset from the <i>mffSM</i> packages.  Data are provided  from the American Association of University Professors (AAUP) on annual faculty salary survey of American colleges and universities, 1995.

```{R, fig.width = 10}
data(AAUP, package = "mffSM")

dim(AAUP)
head(AAUP)
summary(AAUP)


### A data subset containing only relevant variables
aaup <- subset(AAUP, select = c("FICE", "name", "state", "type", "n.prof",
                                "n.assoc", "n.assist", "salary.assoc"))

dim(aaup)
head(aaup)
summary(aaup)

### considering only a subset of covariates
Data <-subset(aaup, complete.cases(aaup[, c("salary.assoc", "n.prof", "n.assoc",
                                            "n.assist", "type")]))
```
<br><br>

### <font color = "#880000">Indidividual work</font> 
<ul>
<li>Consider the dataset on annual income of university professors in USA and perfom simple exploratory analysis.</li>
<li>Focus on average yearly salaries for different types of the universities.</li>
<li>In our exploratory analysis try to take into account also an information about different numbers of univeristy professors (of different ranks)
at the given university.</li>
</ul>
<hr size="3" width="100%" noshade style="color:#440000" align="left" />
<br><br>


```{R, fig.width = 10}
COL3 <- rainbow_hcl(3)
plot(salary.assoc ~ type, data = aaup, col = COL3, xlab = "Type of institution", 
                          ylab = "Salary associate professor [USD 100]")
```



We will fit a few different models that we will compare later. Firstly, a simple additive model: 

```{R, fig.width = 10}
summary(m1orig <- lm(salary.assoc ~ type + n.prof + n.assoc + n.assist,  data = Data))
```

Interpret the estimated parameters in the model and compare it with the interpretation of the estimated parameters 
in the following model: 

```{R, fig.width = 10}

Data <- transform(Data, n.prof40 = n.prof - 40,
                  n.assoc40 = n.assoc - 40,
                  n.assist40 = n.assist - 40)

### Simple additive model with transformed variables 
summary(m1 <- lm(salary.assoc ~ type + n.prof40 + n.assoc40 + n.assist40, data = Data))
```

And, once more the same model, however, with a different parametrization used for the categorical covariate 
of the university type: 


```{R, fig.width = 10}
###  contrast sum ("contr.sum") parametrization for "type"
m1B <- lm(salary.assoc ~ type +n.prof40 + n.assoc40 + n.assist40, data = Data,  contrasts = list(type = "contr.sum"))
```


Finally, one more complicated model with interations (in two different parametrization of the 
continuous covariates): 

```{R, fig.width = 10}
summary(m2orig <- lm(salary.assoc ~ type + n.prof*n.assoc + n.assist,  data = Data)) 
summary(m2 <- lm(salary.assoc ~ type + n.prof40*n.assoc40 + n.assist40, data = Data)) 
```

### <font color = "#880000">Indidividual work</font> 
<ul>
<li>Compare both models and interpret the estimated parameters. </li>
<li>Focus on the estimated salaries for different types of university while controlling for specific 
amounts of the academic personel at the given university.</li>
</ul>
<hr size="3" width="100%" noshade style="color:#440000" align="left" />
<br><br>

### <b>4. Hierarchical structure of the model</b> 
Finally, we will briefly compare a hierarchically well formulated model and a model that is non-hierarchical, Above, we already considered the model for the expected anual salary of the university professors in the USA and the proportion of the professors, associate professors and assistant professors were lowered by 40 (in order to ensure a better interpretation of the intercept parameter and better interpretation of the interaction terms). 

Consider firstly a hierarchically well formulated model with and without the corresponding linear transformation of the number of professors. 

```{R, fig.width = 10}
summary(m3orig <- lm(salary.assoc ~ type + n.assoc * n.prof,  data = Data)) 
summary(m3 <- lm(salary.assoc ~ type + n.assoc40 * n.prof40, data = Data)) 
```

<br><br>
### <font color = "#880000">Indidividual work</font> 
Both models above a hierarchically well formulated models. 
<ul>
<li>Compare the estimated paramters and try to reconstract the parameter estimates from 
one model using the estimates from the second model..</li>
<li>Compare summary statistics for both models (i.e., degrees of freedom, residual sum of squares, $F$-statistic, etc. )</li>
</ul>
<hr size="3" width="100%" noshade style="color:#440000" align="left" />
<br><br>


Now, we will fit two wery similar models but both will be non-hierarchical. Compare both models 
and try to understand the effect of the transformation used. 

```{R, fig.width = 10}
summary(m3orig <- lm(salary.assoc ~  - 1 + type + n.assoc * n.prof - n.prof,  data = Data)) 
summary(m3 <- lm(salary.assoc ~  - 1 + type + n.assoc40 * n.prof40 - n.prof40, data = Data)) 
```


<br><br>
