---
title: "V.\\ Statistical analysis in R (presentation)"
author: "Data Science Lab, University of Copenhagen"
date: today
format: 
  html:
    toc: true
    toc-depth: 3
    toc-location: left    
    embed-resources: true
    standalone: true
execute:
  warnings: true
  messages: true
  echo: true
---


## Structure of a statistical analysis in R

The very basic structure of an R script doing a classical statistical analysis is as follows:

1. Load packages that you will be using.

1. Read the dataset to be analyzed. Possibly also do some data cleaning and manipulation.

1. Visualize the dataset by graphics and other descriptive statistics.

1. Fit and validate a statistical model.

1. Hypothesis testing. Possibly also post hoc testing.

1. Visualize and/or report: model parameters and/or model predictions.

Of course there are variants of this set-up, and in practice there will often be some iterations of the steps.

In this manuscript we will examplify the proposed steps in the analysis of two simple datasets:

A) One-way ANOVA of the expression of the IGFL4 gene against the skin type in psoriasis patients.

A) Multiple linear regression of volume of cherry trees against their diameter and height.


#### Step 1: Load packages

For both examples we will use `ggplot2` to make plots, and to be prepared for data manipulations we simply load this together with the rest of the `tidyverse`. Moreover, we will also use extra plotting functionalities from the `ggfortify` and `GGally` packages.

The psoriasis data are provided in an Excel sheet, so we also load `readxl`. Finally, we will use the package `emmeans` to make post hoc tests and to report model predictions.

Remember that you should install the wanted packages before they can be used (but you only need to install the packages once!) To take full advantage of the `emmeans` package you should also install the packages `multcomp` and  `multcompView` packages.

Thus,
```{r}
#| warning: false
#| message: false
#install.packages("tidyverse")
#install.packages("ggfortify")
#install.packages("GGally")
#install.packages("readxl")
#install.packages("emmeans")
#install.packages("multcomp")
#install.packages("multcompView")
#install.packages("lme4")
library(tidyverse)
library(ggfortify)
library(GGally)
library(readxl)
library(emmeans)
library(multcomp)
library(multcompView)
library(lme4)
```

Now we have done step 1 for both analyses. Next we will look specifically at the two examples. Finally, we conclude with a brief outlook on other statistical models in R.


## Example A: Analysis of variance

#### Step 2: Data

Psoriasis is an immune-mediated disease that affects the skin. Researchers carried out an micro-array experiment with skin from 37 people in order to examine a potential association between the disease and a certain gene (IGFL4). For each of the 37 samples the gene expression was measured as an intensity. Fifteen skin samples were from psoriasis patients and from a part of the body affected by the disease (`psor`); 15 samples were from psoriasis patients but from a part of the body not affected by the disease (`psne`); and 7 skin samples were from healthy people (`control`).  

The data are saved in the file **psoriasis.xlsx**. At first the variable `type` is stored as a character variable we change it to a factor (and check that indeed there are 15, 15 and 7 patients in the three groups). 

```{r}
psorData <- read_excel("psoriasis.xlsx")
psorData
psorData <- mutate(psorData, type = factor(type)) 
count(psorData, type)
```


#### Step 3: Descriptive plots and statistics

To get an impression of the data we make two plots and compute group-wise means and standard deviations.

```{r}
ggplot(psorData, aes(x=type, y=intensity)) +
  geom_point() + 
  labs(x="Skin type", y="Intensity")

ggplot(psorData, aes(x=type, y=intensity)) +
  geom_boxplot() + 
  labs(x="Skin type", y="Intensity")

psorData %>% 
  group_by(type) %>% 
  summarise(avg=mean(intensity), sd=sd(intensity))
```


#### Step 4: Fit of oneway ANOVA, model validation

The scientific question is whether the gene expression level differs between the three types/groups. Thus, the natural type of analysis is a oneway analysis of variance (ANOVA). The oneway ANOVA is fitted with `lm`. It is a good approach to assign a name (below *oneway*) to the object with the fitted model. This object contains all relevant information and may be used for subsequent analysis. Note that we have logarithmic transformed the response as intensities often live on a multiplicative scale. 

```{r}
oneway <- lm(log(intensity) ~ type, data=psorData)
oneway
```

Mathematically this model makes assumptions about linearity (which, however, will always be satisfied for a oneway ANOVA -- why?), variance homogeneity, and normality. These assumptions are checked visually via four plots generated from the model residuals. In R it is easy to make these plots; simply apply `plot` on the `lm`-object.

```{r}
par(mfrow=c(2,2))   # makes room for 4=2x2 plots!
plot(oneway)
```

The model validation plots do not show any evident problems. 

**Remark:** Using the `autoplot()` function from the `ggfortify` package you can also make a ggplot-version of these plots. This can be aesthetically more pleasing, you avoid calling the `par()` function, and there is more room for the plot on a small laptop screen. But mathematically it is the same plots that will be produced!
```{r}
autoplot(oneway)
```


#### Step 5: Hypothesis test + Post hoc tests

It is standard to carry out an $F$-test for the overall effect of the explanatory variable. To be precise, the hypothesis is that the expected values are the same in _all_ groups. The most easy way to do this test is to use `drop1`. The option `test="F"` is needed to get the $F$-test:

```{r}
drop1(oneway,test="F")
```

Thus, the overall test for homogeneity between the groups show, that the groups are not all the same. But it might be that the gene expression in two of the three groups, say, are not significantly different. To investigate that we do post hoc testing. This is nicely done within the framework of *estimated marginal means* using the `emmeans` package. Here `emmeans` makes the estimated marginal means (that is the predicted gene expression intensity on the log scale), and the `pairs` command provide post hoc pairwise comparisons:

```{r}
emmeans(oneway, ~type)
pairs(emmeans(oneway, ~type))
```

It can be convenient to summarize the pairwise comparisons in the so-called *compact letter display*. Here a grouping is performed such that groups _not sharing_ a letter are significantly different. Previous versions of the **emmeans** package had a function for doing this (called `CLD`). However, the author of the **emmeans** package does not approve of this way of displaying the results, and consequently he has removed the `CLD` function from the package! But instead you can use the `cld` function from the **multcomp** package. Notice that both the **multcomp** and the **multcompView** package must be installed.

Using the syntax `multcomp::cld` we may avoid loading the **multcomp** package. Long story short; here is the code:

```{r}
cld(emmeans(oneway, ~type))
```

We see that the "grouping letters" appearing in the column `.group` actually are numbers. But the interpretation is the same. Namely that `psor` (letter=1) is significantly different from both `healthy` (letter=2) and `psne` (letter=2), whereas that two latter are not significantly different (they share the letter "2"). Note: If you want actual letters instead of digits, then add the option `Letters=letters`.

Or restated in biological terms: The groups *healthy* and *psne* appears to have the same expression of the IGFL4 gene, and this expression is significantly larger than in the *psor* group. Biologically this is consistent with the IGFL4 gene being related to psoriasis as *psne* corresponds to skin not affected by the disease.


#### Step 6: Report of model parameters

In this example it is natural to report the estimated intensities as well as the comparisons between the 3 groups. As the analysis was done using a log transformation it is also advisable to backtransform. Both things can be done automatically by the **emmeans** package, where the option `type="response"` (here `type` is the name of an option inside the R command, and it is a coincidence that it has the same name as the variable *type*) requests the backtransformation:

```{r}
emmeans(oneway, ~type, type="response")
confint(pairs(emmeans(oneway, ~type, type="response")))
```

Although we don't recommend usage of standard errors in the context of backtransformed parameters, we see that a standard error is provided on the backtransformed parameters.

#### Iteration of step 3 to step 6

Actually the used statistical model is wrong! At least there appears to be a strong correlation between the n'th obserservations of *psne* of *psor*. This is evident from the following two descriptive statistics:

```{r}
# correlation
cor(psorData$intensity[1:15], psorData$intensity[16:30])
# scatter plot
qplot(x=psorData$intensity[1:15], y=psorData$intensity[16:30]) + 
  geom_point() + 
  geom_abline(intercept=0,slope=1) + 
  coord_equal(xlim=c(0.6,2.1),ylim=c(0.6,2.1))
```

What is going on here?

The first 30 observations presumably do not come from 30 different people, but are rather two measurements from each of 15 people, one from an affected skin area and one from non-affected area. In that case, the correct analysis should include a random effect of person. First we introduce the apparently missing variable in the dataset:

```{r}
psorData$Id <- factor(c(1:15,1:15,16:22))
```

Then we can do the correct mixed effects model using the `lmer()` function from the **lme4** package:
```{r}
oneway_with_random_effect <- lmer(log(intensity) ~ type + (1|Id), data=psorData)
```

The syntax used in the other steps now change a bit at some places due to the usage of the `lme4` package. We will not do all the steps here, but for comparison let's check that the preferred analysis including the random effect indeed is more powerful:

```{r}
drop1(oneway, test="F")
drop1(oneway_with_random_effect, test="Chisq")
```


## Example B: Simple and multiple linear regression

#### Step 2: Data

For the first part of the presentation we will use data from 31 cherry trees available in the built in dataset `trees`: Girth, height and tree volume have been measured for each tree. Please note, that the help page `?trees` states, that girth is the tree diameter. Interest is in the association between tree volume on the one side and girth and height on the other side. 

```{r}
data(trees)
```

#### Step 3: Visualization of raw data

We start by plotting the data. If you `plot` a data frame, then you get a display of pairwise scatter plots (using Base R graphics):

```{r}
plot(trees)
```

A more fancy version of this is available in **GGally**, where the diagonal is used to visualize the marginal distribution of the three tree variables:

```{r}
ggpairs(trees)
```

#### Step 4: Fitting and validating a simple linear regression model

We start out with a simple linear regression with volume as outcome and girth as covariate. This model is fitted with the `lm` function. It is a good approach to assign a name (below *linreg1*) to the object with the fitted model. This object contains all relevant information and may be used for subsequent analysis.

```{r}
linreg1 <- lm(Volume~Girth, data=trees)
linreg1
```

Mathematically this model makes assumptions about linearity, variance homogeneity, and normality. These assumptions are checked visually via four plots generated from the model residuals. In R it is easy to make these plots; simply apply `plot` on the `lm`-object.

```{r}
par(mfrow=c(2,2))   # makes room for 4=2x2 plots!
plot(linreg1)
```

In this case there appears to be problems with the model. In particular, the upper left plot shows a quadratic tendency suggesting that the linearity assumption is not appropriate. Perhaps a data transformation can help us out.

#### Step 4 iterated: Transformation

The linear regression model above says that an increment of $\Delta_g$ of girth corresponds to an increment of volume by $\beta\Delta_g$ where $\beta$ is the slope parameter, i.e., interpretation is concerned with _absolute changes_. Perhaps it is more reasonable to talk about _relative changes_: An increment of girth by a factor $c_g$ corresponds to an increment of volume of $\gamma c_g$. This corresponds to a linear regression model on the log-transformed variables, which is easily fitted. Notice that `log` in R is the natural logarithm.

```{r}
linreg2 <- lm(log(Volume) ~ log(Girth), data=trees)
par(mfrow=c(2,2))
plot(linreg2)
```

The statistical validity of this model is much better! Which is nice as it also has a much more meaningful biological/physical interpretation!


#### Step 6: Report of the model

Parameter estimates, standard errors are computed and certain hypothesis tests are carried out by the `summary()` function:

```{r}
summary(linreg2)
```

The coefficient table is the most important part of the output. It has two lines --- one for each parameter in the model --- and four columns. The first line is about the intercept, the second line is about the slope. The columns are

* `Estimate`: The estimated value of the parameter (intercept or slope).
* `Std. Error`: The standard error (estimated standard deviation) associated with the estimate in the first column.
* `t value`: The $t$-test statistic for the hypothesis that the corresponding parameter is zero. Computed as the estimate divided by the standard error.
* `Pr(>|t|)`: The $p$-value associated with the hypothesis just mentioned. In particular the $p$-value in the second line is for the hypothesis that there is no effect of girth on volume.

Below the coefficient table you find, among others, the `Residual standard error`, i.e., the estimated standard deviation for the observations.

If you prefer to report confidence intervals for the parameters (intercept and slope) instead of standard error then use `confint()`:

```{r}
confint(linreg2)
```

Finally, don't forget to make the mathematical interpretation of the model and its parameter estimates! Thus, the model postulates the relation $$\ln V \approx \alpha + \beta \ln G$$ with estimates $\hat{\alpha} = -2.353$ and $\hat{\beta} = 2.200$. This corresponds to $V \approx e^\alpha \cdot G^\beta$, and we see that an increment of girth by 10%, say, corresponds to an increment of volume by a factor of $1.1^{2.2} = 1.23$, that is, by 23%.


#### Step 6: Visualization of the model

An excellent synthesis of a statistical analysis is often provided by a plot combining the raw data with the model fit.

The estimated regression line can be added to a scatter plot in different ways (and with/without corresponding standard error curves):

```{r, fig.height=3, fig.width=4}
ggplot(trees, aes(x=log(Girth), y =log(Volume))) + 
  geom_point() +
  geom_abline(intercept=-2.35332, slope=2.19997, col="red")

ggplot(trees, aes(x=Girth, y=Volume)) + 
  geom_point() +
  geom_smooth(method="lm", se=TRUE) +
  scale_x_log10() +
  scale_y_log10()
```



#### Step 4 iterated: Multiple linear regression

Several covariates can be included in a regression model, corresponding to multiple linear regression; simply write a plus between the covariates. Here is fitting of a model with log-girth as well as log-height included as explanatory variables:

```{r}
linreg3 <- lm(log(Volume) ~ log(Girth) + log(Height), data=trees)
summary(linreg3)$coefficients
```

Don't forget to do model validation:
```{r}
par(mfrow=c(2,2))
plot(linreg3)
```


#### Step 5: Hypothesis test

The $F$-test for comparison of two nested models may be carried out by the `anova` function. For example, `linreg3` is nested in `linreg2`, and the comparison between the two corresponds to the hypothesis that there is no extra information in height when girth is included (on log-scale).

```{r}
anova(linreg3, linreg2)
```

However, as an alternative to `anova` it is more simple to use the `drop1` command, which makes all the tests corresponding to removal of one of the explanatory variables. If you want $p$-values, then ask for $F$-tests:
```{r}
drop1(linreg3,test="F")
```

We see that `anova` and `drop1` give the same p-value ($p=7.8 * 10^{-6}$) for the effect of log(Height). And in a Gaussian model (like this) this coincides with the p-value from the $t$-test for the hypothesis that there is no effect of a single numerical covariate. Thus, 

```{r}
summary(linreg3)$coefficients
```

So should you do all three tests? Definitely not! 

* In general it is not recommendable to use `summary` for doing tests: The $F$-tests and $t$-tests only coincide for hypothesis on a single parameter. And you may also be deceived by the model parametrizations when interpreting $t$-tests.
* To use `drop1` you only need to fit the large model (here *linreg3*). Thus, it gives shorter and more easily readable R code.

## Use of pipe operator

Although we did not use it systematically in the above, the pipe operator may increase readability of the code. For example, the two commands below do exactly the same: 

```{r}
confint(pairs(emmeans(oneway, ~type, type="response")))
oneway %>% emmeans(~type, type="response") %>% pairs() %>% confint()
```


## Outlook: Other analyses

The `lm` function is used for linear models, that is, models where data points are assumed to be independent with a Gaussian distribution (and typically also with the same variance). Obviously this class of models is not always appropriate, and there exists functions for many, many more situations and data types. Here we just mention a few functions corresponding to common data types and statistical problems. 

* `glm()`: For independent, but non-Gaussian data. Examples are binary outcomes (logistic regression) and outcomes that are counts (Poisson regression). glm is short for generalized linear models, and the `glm()` function is part of the base installation of R. Remark: SAS users should not confuse this with **PROC GLM**, which does linear normal models like `lm()`.
* `lmer()` and `glmer()`: For data with dependence structures that can be described by random effects, e.g., block designs. lme is short for linear mixed effects  (Gaussian data), glmer is short for generalized linear mixed effects  (binary or count data). Both functions are part of the **lme4** package. Actually, we did touch upon a mixed effects model above.
* `nls()`: For non-linear regression, e.g., dose-response analysis. nls is short for non-linear least squares. The function is included in the base installation of R.

The functions mentioned above are used in a similar way as `lm()`: a model is fitted with the function in question, and the model object subsequently examined with respect to model validation, estimation, computation of confidence limits, hypothesis tests, prediction, etc. with functions `summary()`, `confint()`, `drop1()`, `predict()`, `emmeans()`, `pairs()` as indicated above. 

**End of presentation.**

