---
title: "V.\\ Statistical analysis in R (solution)"
author: "Data Science Lab, University of Copenhagen"
date: "August 2025"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

```{r, include=FALSE}
library(emmeans)
```


### Data

1. We load the MASS package to get access to the `birthwt` dataset and make the relevant manipulations of the dataset. In the end, notice how levels 3-6 of `ftvFac` are collapsed into one level denoted "MoreThanOnce". 
 
    ```{r, message=FALSE}
# Load MASS package
library(MASS)

### Open help page for birthwt
?birthwt
    ```

2. The `tidyverse` package is needed to make tibbles, so it is loaded first. Then, we construct the tibble. 
   
    ```{r, message=FALSE}
### Make a tibble with the data
library(tidyverse)
birthData <- as_tibble(birthwt)
birthData
    ```

3. We recode the `smoke` variable from a numerical variable to a categorical one.

    ```{r}
### Make smoke into a factor
birthData <- mutate(birthData, smoke = factor(smoke))
    ```

4. We examine the `ftv` variable, then construct a categorical version of the `ftv` variable, then make a version where levels 2 and above are collapsed into just one level. 

    ```{r}
# Check the ftv variable, make it a factor and collapse some of the levels
table(birthData$ftv)
birthData <- mutate(birthData, ftvFac = factor(ftv))
birthData <- mutate(birthData, 
                    visits = fct_collapse(ftvFac, Never="0", Once="1", other_level="MoreThanOnce"))
table(birthData$visits, birthData$ftvFac)
birthData %>% group_by(ftvFac, visits) %>% count()
    ```

5. We make boxplots for `bwt` for different levels of one or more variables. Birth weight appears to be smaller for infants born by smokers compared to non-smokers, while birth weight does not appear to depend on the frequency of vistis to the doctor. This is just preliminary consideratiions; formal analyses are carried out later. 

    ```{r, message=FALSE, fig.height=3, fig.width=4}
### Parallel boxplots
ggplot(birthData, aes(x=visits, y=bwt)) + geom_boxplot() 
ggplot(birthData, aes(x=smoke, y=bwt)) + geom_boxplot() 
ggplot(birthData, aes(x=visits, y=bwt, color=smoke)) + geom_boxplot() + facet_wrap(~smoke) 
    ```

6. Scatterplot with `lwt` on the $x$-axis and `bwt` on the $y$-axis. The color and shape of the points can be controlled by the options `color` and `shape`. There is no clear association between the two variables (mother's weight before pregnancy and infant weight), but let us see in the analysis. 
 
    ```{r, message=FALSE, fig.height=3, fig.width=5}
### Scatter plot
ggplot(birthData, aes(x=lwt, y=bwt)) + geom_point() 

### Points coloured after visits
ggplot(birthData, aes(x=lwt, y=bwt, color=visits)) + geom_point() 

### Points coloured after visits, point types after smoke status
ggplot(birthData, aes(x=lwt, y=bwt, color=visits, shape=smoke)) + geom_point() 
    ```

### Regression

7. We use the `lm` function to fit the simple linear regression model with infant's birth weight as response and mothers weight before pregnancy as predictor. Then, we extract estimates and confidence intervals (CIs). Specifically we get a 95% CI 1.05-7.81 for the slope parameter. Notice that zero is not included in the CI, in line with the p-value of 0.01 which appears in the summary. So, despite the impression from the scatter plot above, there is a significant association between the two variables. 

    ```{r}
### Simple linear regression
reg1 <- lm(bwt ~ lwt, data=birthData) 
summary(reg1)
confint(reg1)
    ```

8. We make the usual model validation plots. None of them gives us any reasons to worry about the validity of the model assumtions. 

    ```{r}
### Model validation
par(mfrow=c(2,2))
plot(reg1)
    ```

9. We include extra predictors in the model, so we now get a multiple linear regression. The regression coefficients are estimated to 4.12 for mother's weight (`lwt`), 7.49 for age and 15.4 for number of visits to the doctor (`ftv`). The interpretation of the last one, say, is that the birth weight is expected to increase by 15.4 grams per extra visit by the doctor, all other predictors kept unchanged. Notice that the estimate is not significantly different from zero (p=0.76) though. 

    ```{r}
### Include age and ftv as covariate
reg2 <- lm(bwt ~ lwt + age + ftv, data=birthData)
summary(reg2)$coefficients
    ```

10. We make a new (boring dataset) with just one data line and values of each of the three predictors from the multiple linear regression from the previous question. The `predict` function computes the corresponding predicted value for the birth weight. It needs to know the predictor values as well as which model to use for the prediction (you could try using `reg1` instead of `reg2`; then you would get a slightly different prediction). The `interval` option makes it possible to supplement the point prediction with a 95% prediction interval, i.e., an interval that contains the birth weight for 95% of infants with the given predictor characteristics. 

    ```{r}
### Prediction
newData <- data.frame(lwt=100, age=25, ftv=0)
newData
predict(reg2, newData)
predict(reg2, newData, interval="prediction")
    ```

11. We use a different dataset for the regression, namely a subset constructed with the `filter` function. There is just for the sake of the example, so we refrain from looking at the results in any detail. 

    ```{r}
### Use only data from mothers with weight below 160
reg3 <- lm(bwt ~ lwt + age + ftv, data=filter(birthData, lwt<160))
summary(reg3)$coefficients
    ```


### ANOVA

12. We use `lm` to fit the oneway ANOVA. The non-smokers (`smoke=0`) are used as reference, so the expected birth weight for that group is given by the intercept, 3056 grams. For the estimate for smokers, we must add the contrast/difference caused by the smoke variable (which is negative). We get 3056-284 = 2772. With `emmeans` we get both numbers automatically (but then not the difference). The p-value for the effect of smoking is 0.0087, which can be found both from the `summary` and the `pairs` output, so there is a statistically significant effect of smoking.
  
    ```{r}
### Oneway ANOVA against smoke
oneway1 <- lm(bwt ~ smoke, data=birthData)
summary(oneway1)$coefficients
emmeans(oneway1,~smoke)
pairs(emmeans(oneway1,~smoke))
    ```

13. We now fit a another oneway ANOVA where `visits` from question 4 is used as the explanatory variable. The groupwise estimated birth weights are most easily found with `emmeans` and are 2865 g, 3108 g, and 2951 g, respectively. The output from `pairs` shows that none of the groups differ statistically significantly from each other (p-valuer are 0.14 or larger). Finally, the overall test gives a p-value of 0.17 (`drop1`), confirming that groups do not differ significantly. 
 
    ```{r}
### Oneway ANOVA against visits
oneway2 <- lm(bwt ~ visits, data=birthData)
emmeans(oneway2,~visits)
pairs(emmeans(oneway2,~visits))
drop1(oneway2,test="F")
    ```

14. A twoway ANOVA with both `smoke` and `visits` is fitted. It is the model _without_ interaction because there is a plus (+) between the two variables. Reference groups are chosen for each variable (`smoke=0` corresponding to non-smokers and `vists=Never`), and the intercept 2981 g is the estimated birth weight for the reference combination. For the contrast, the birth weight is estimated to decrease by 257 g if the mother is a smoker, and to increase by 193 g and 74 g if the mother visits the doctor once or more than once, respectively. We also get p-values for each of these contrasts, and conclusions are unchanged from the oneway ANOVAs: There is a significant effect of smoking but not a significant effect of visiting the doctor. 

    ```{r}
### Twoway ANOVA without interaction
twoway1 <- lm(bwt ~ visits + smoke, data=birthData)
summary(twoway1)$coefficients    
    ```

15. The twoway ANOVA _with_ interaction is fitted by using a multiplication sign or star (*) between the two variables. This allows for the effect os smoking to differ between visits groups, or vice versa. The hypothesis of no interaction can be carried out with the `anova` function, where the models to be compared (the models with and without interaction) are specified explicitly, or with `drop1`. The same test is carried out, so we get the same p-value, 0.14. So, the interaction effect is not significant. By the way, the function name `anova` is a bit annoying since it can be used for many model types, not only ANOVA models. 

    ```{r}
### Twoway ANOVA with interaction, test for interaction in two ways
twoway2 <- lm(bwt ~ visits * smoke, data=birthData)
anova(twoway2, twoway1)
drop1(twoway2,test="F")
    ```

16. The `emmeans` function can be used to compute group-specific estimates, on average over the other variables in the model. We get the estimated birth weights of 3066 g and 2772 g, respectively. This is not exactly the same as in question 12: The estimates in question 12 were simple averages in the two groups, while the new estimates are averages of estimates. This is not the same when the dataset is unbalanced wrt. to the two variables.  

    ```{r}
emmeans(twoway2, ~smoke)
    ```

### Models with numerical as well as categorical predictors

17. `model1` includes both a numerical and a categorical explanatory variable. Birth weight is estimated 
to increase 4.24 g for each extra pound of mother's weight (no matter the smoking status) and to differ 272 g for smoking and non-smoking mothers. The intercept is the estimated birth weight for a non-smoking mother with a weight of zero pounds! This is of course nonsense - as is often the case with regression type models when zero values of a numerical predictor cannot occur in practice. 

    ```{r}
### Model with linear (lwt,bwt) association. 
### Intercept differ between smokers and non-smokers, one common slope.
model1 <- lm(bwt ~ lwt + smoke, data=birthData)
summary(model1)
    ```

18. Interaction between a numerical and categorical variable is interpreted as different slopes (wrt. the numerical variable) for the different groups. Here, this means two straight lines (`x=lwt` and `y=bwt`), one for smokers and one for non-smokers, with different intercepts and different slopes. The model from question 17 had different intercepts, but the same slope. Both the intercept and the slope are reported in terms of a reference groups (non-smokers) and a difference between smokers and non-smokers, and the two estimated regression lines are therefore given by: bwt = 2351 + 5.39 * lwt for non-smokers and bwt = (2351+41) + (5.39-2.42) * lwt = 2392 - 2.97 * lwt for smokers. The interaction term, as reported the last line in the coefficients table, estimates the difference in the weight effect between non-smokers and smokers, so this is where we should look to answer the last question. The p-value is 0.48, so there is no evidence that the effect of mother's weight affects the birth weight is different for smokers and non-smokers. The same test can be carried out by comparing to the additive model from question 17, using the `anova` function. 

    ```{r}
### Model with linear (lwt,bwt) association. Intercept and slope both differ between visit groups
model2 <- lm(bwt ~ lwt * smoke, data=birthData)
summary(model2)
### Test if slopes differ between visit groups
anova(model2, model1)
    ```

19. We now include more variables in the model. The model validation plots cause no worries. Smoking is estimated to decrease infant birth weight by 244 g (SE 107 g), and the effect is statistically significant (p=0.017).


    ```{r}
### Model with many effects (no interactions)
model3 <- lm(bwt ~ lwt + smoke + age + visits, data=birthData)
par(mfrow=c(2,2))
plot(model3)
summary(model3)
    ```

### Logistic regression

20. We now fit a logistic regression, modelling the probability of a low birthweight defines as <2500 g. The effect of smoking is just not significant (p=0.06), which most likely reflects that we loose information going from the actual birth weights to the binary version. Notice that the estimated coefficient for smoking is now positive (0.62). This means that smoking increases the probability of getting a child with low birth weight; hence the direction of a potential effect is still the same: Smoking increases the likelihood of a smaller birth weight. 

    ```{r}
### Logistic regression with many predictors (no interactions)
logreg1 <- glm(low ~ lwt + smoke + age + visits, data=birthData, family="binomial")
summary(logreg1)
    ```


### Linear mixed models

21. We first make the artificial center variable as described in the exercise. You should of course never make such a construct unless it reflects the data collection. 

    ```{r}
### Make artificial center variable
set.seed(123)
center <- sample(rep(1:19, each=10)[1:189])
birthData <- mutate(birthData,center=factor(center))
    ```

22. We use the new variable as a random factor in a linear mixed model (LMM) and identify the relevant estimate for smoking to be 244 g (SE 107). This is exactly the same as in the model without random effects (question 19). That would typically not be the case, but happens here because the center-to-center standard deviation is essentially estimated to zero (3.109e-06), which is not too surprising since the center variable had nothing to do with actual data collection. Notice that the summary of an `lmer` object does not provide p-values. There are ways to get those p-values, for example by using `emmeans` or the **lmerTest** package instead of the the **lme4** package, but we will not go into details about that here.

    ```{r, message=FALSE}
### Remember to install lme4 before this can run
# install.packages("lme4")
library(lme4)

### Linear mixed model with random effect of center
lmm1 <- lmer(bwt ~ lwt + smoke + age + visits + (1|center), data=birthData)
summary(lmm1)
    ```

**End of solution**
