```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE,tidy=TRUE,message=FALSE,warning=FALSE,results="hold")
showsol<-FALSE
```

(Latest changes: 08.11.2018 - typos).

(Warning: some changes may occur before the second week.)

# Overview

**Aim**: Present methods for analysing correlated responses in a (normal/Gaussian) regression setting.

We will only consider _two-level models_ and in particular focus on _random intercept and random slope models._

---

## Learning material

* Textbook: Fahrmeir et al (2013): Chapter 2.4, 7.1-7.3, 7.7. In greater detail: pages 349-354 (not "Alternative view on the random intercept model"), 356-365 (not 7.1.5 "Stochastic Covariates""), 368-377 (not "Bayesian Covariance Matrix""), 379-380 (not "Testing Random Effects or Variance Parameters"", only last part on page 383), 383 (middle), 401-409 (orange juice). Note: Bayesian solutions not on the reading list.

* Alternative readings: Zuur et al. (2009): "Mixed Effects Models and Extensions in Ecology with R", chapter 5 (pages 101-142). Available as free ebook from Springer for NTNU students. More explanations and less mathematics than Fahrmeir et al (2013), more focus on understanding. [Link to ebook Chapter 5](https://link.springer.com/chapter/10.1007/978-0-387-87458-6_5)

* [Classnotes 01.11.2018](https://www.math.ntnu.no/emner/TMA4315/2018h/M7w12018.pdf) 
* [Classnotes 08.11.2018](https://www.math.ntnu.no/emner/TMA4315/2018h/M7w22018.pdf) 


#<a id="interactivew1">Interactive session week 1 </a>

Plan:  

* Problem 1abcd (not ef - covered in week 2).
* Problem 2abe (not cd - covered in week 2).
* If time: [Exam 2017, Problem 3](#exam2017)

---

## Exercise 1: Taken from UiO, STK3100, 2011, problem 1

We will in this exercise look at the following model:

$$Y_{ij} = \mathbf{x}_{ij}^T \pmb{\beta} + \gamma_{0i} + {\varepsilon}_{ij}, \ \varepsilon_{ij} \sim N(0, \sigma^2) $$
$$\gamma_{0i} \sim N(0, \tau_0^2)$$

where $i \in \{1, \dots, m\}$ is the group index, and $j \in \{1, \dots, n_i\}$ is the index for repeated measurements in each group (cluster). We assume that all random variables are independent.

<span class="question">a) What is this model called? Discuss where and when such models are useful</span>


<details><summary>Answer</summary>
This is a random intercept model. Such models are useful if you want to model correlations between variables from the same group (clustered data) or the same invidual (longitudinal data).
</details>

<span class="question">b) What is the *marginal* model for $\mathbf{Y}_i = (Y_{i1}, \dots, Y_{in_i})$? Show how you arrive at your answer and write with vectors and matrices (for cluster $i$).</span>

<span class="question">What are the advantages of having an expression for the marginal distribution of $\mathbf{Y}_i$ when doing estimation? </span>

<details><summary>Answer</summary>
We have 

$$\mathbf{Y}_i = \mathbf{X}_i \pmb{\beta} + \mathbf{1}\gamma_{0i} + \pmb{\varepsilon}_i, \ \pmb{\varepsilon}_i \sim N(0, \sigma^2\mathbf{I}) $$

which gives

$$ \mathbf{Y}_i \sim N(\mathbf{X}_i\pmb{\beta}, \tau_0^2 \mathbf{11}^T + \sigma^2 \mathbf{I}) $$

where $\mathbf{1}$ is a column vector of ones. 

Advantage of having the marginal distribution: Then we also have an explicit expression for the likelihood

$$L(\pmb{\beta}, \tau_0^2, \sigma^2) = \prod_{i = 1}^N f(y_i; \pmb{\beta}, \tau_0^2, \sigma^2) $$

which we need to find MLE and REML estimates.

</details>

<!-- New: also write down the formula for ${\bf Y}$ (all clusters together). -->



In the rest of the exercise we will look at a dataset consisting of observations of cod from the year 2000. The dataset is from Havforskningsinstitutten in Bergen, and we use only a subset of a large dataset with information about fish in Barentshavet. The following variables are available on a random sample within each catch (a draft of trawl (trål)):

* `length`: the length of fish (in cm)
* `weight`: the weight of fish (in grams)
* `age`: the age of fish (in years)
* `haulsize`: the size of the total catch (in ton (1000 kg))

Let $i$ be the index for catch, and $j$ be the index for an individual fish in a catch. We will not use the age variable in this exercise.

We start by looking at the model

$$ \log(\texttt{weight}_{ij}) = \beta_0 + \beta_1 \log(\texttt{length}_{ij}) + \beta_2\log(\texttt{haulsize}_i) + \gamma_{0i} + \varepsilon_{ij} $$

where $\gamma_{0i}\sim N(0, \tau_0^2)$, $\varepsilon_{ij} \sim N(0,\sigma^2)$, and all random effects are independent. Below is an excerpt from the output from fitting this model in `R`.

```{r, echo = TRUE,eval=TRUE}

library(lme4)

filepath <- "https://www.math.ntnu.no/emner/TMA4315/2018h/fishdata.dat"
fish <- read.table(filepath, header = TRUE)

fit1 <- lmer(log(weight) ~ log(length) + log(haulsize) + (1 | haul), data = fish)

summary(fit1)

```


<span class="question">c) Write down the estimates for all the parameters in the model. What is the correlation between two weight-variables from the same catch (haul)?</span>

<details><summary>Answer</summary>

```{r, echo = FALSE}

library(lme4)

filepath <- "https://www.math.ntnu.no/emner/TMA4315/2018h/fishdata.dat"
fish <- read.table(filepath, header = TRUE)

fit1 <- lmer(log(weight) ~ log(length) + log(haulsize) + (1 | haul), data = fish)

#summary(fit1)

```

We have the following parameter estimates:

```{r, echo = FALSE}

df <- data.frame(Parameter = c("$\\beta_0$", "$\\beta_1$", "$\\beta_2$", "$\\sigma_b$", "$\\sigma$"),
                 Estimate = c(fixef(fit1), sqrt(summary(fit1)$varcor$haul[1]), sigma(fit1)))

library(knitr)
kable(df, digits = 4)

```

The correlation between two weight-variables from the same catch (haul) is found from inserting parameter estimates into 
$$\text{Corr}(Y_{ij},Y_{il})=\frac{\tau_0^2}{\tau_0^2+\sigma^2} \text{ for }j\neq l$$

```{r}
vc <- VarCorr(fit1)
print(vc)
df <- as.data.frame(vc)
numerator <- df[4][[1]][1]
denominator <- sum(df[4][[1]])
cat("numerator:", numerator, "\n")
cat("denominator:", denominator, "\n")
ICC <- numerator/denominator
ICC
```

Then, $$\widehat{\frac{\tau_0^2}{\tau_0^2 + \sigma^2}} = \frac{ {`r numerator`} }{ {`r denominator`} }={`r ICC`}$$

</details>

<span class="question">d) We are now interested in $\theta = \exp(\beta_0 + \beta_1 \log(66) + \beta_2 \log(0.46))$. The standard deviation of $\hat{\beta}_0 + \hat{\beta}_1 \log(66) + \hat{\beta}_2 \log(0.46)$ is 0.007. Explain how you can calculate this value (you do not have to do the calculation). Then create and calculate a 95 % confidence interval for $\theta$.</span>

<details><summary>Answer</summary>
$\pmb{\hat{\beta}}$ is approximately Gaussian (normal), and then

$$\log(\hat{\theta}) = \hat{\mu} = \hat{\beta}_0 + \hat{\beta}_1 \log(66) + \hat{\beta}_2 \log(0.46)$$

is also Gaussian with variance

$$Var(\hat{\mu}) = Var(\hat{\beta}_0) + \log(66)^2Var(\hat{\beta}_1) + \log(0.46)^2Var(\hat{\beta}_2) + \\ 2\cdot \log(66)Cov(\hat{\beta}_0, \hat{\beta}_1) + 2\cdot \log(0.46)Cov(\hat{\beta}_0, \hat{\beta}_2) + 2\log(66)\log(0.46)Cov(\hat{\beta}_1, \hat{\beta}_2)) = 0.00706^2$$

Then we get a 95 % confidence interval for $\mu$:

$$\hat{\mu} \pm 1.96 \cdot \sqrt{Var(\hat{\mu})} = [7.849683, 7.877365] $$

We take the exponential of this interval to get a confidence interval for $\theta$: $[2564.922, 2636.914]$

```{r}

# it is much simpler to just use vcov(fit1), but here we write down everything so we can compare with the printout
vbeta0 <- summary(fit1)$coefficients[1,2]^2
vbeta1 <- summary(fit1)$coefficients[2,2]^2
vbeta2 <- summary(fit1)$coefficients[3,2]^2
# cov2cor gives correlations from a covariance matrix. This is where we could just use vcov()
cbeta01 <- sqrt(vbeta0)*sqrt(vbeta1)*cov2cor(vcov(fit1))[1,2]
cbeta02 <- sqrt(vbeta0)*sqrt(vbeta2)*cov2cor(vcov(fit1))[1,3]
cbeta12 <- sqrt(vbeta1)*sqrt(vbeta2)*cov2cor(vcov(fit1))[2,3]

variance <- vbeta0 + log(66)^2*vbeta1 + log(0.46)^2*vbeta2 + 2*log(66)*cbeta01 + 2*log(0.46)*cbeta02 + 2*log(66)*log(0.46)*cbeta12
SD <- sqrt(variance)

# OR

SD2 <- as.numeric(sqrt(c(1, log(66), log(0.46)) %*% vcov(fit1) %*% c(1, log(66), log(0.46))))

cat(SD, SD2)

critval <- qnorm(0.975)

CI <- c(-critval*SD, critval*SD) + c(fixef(fit1)%*%c(1, log(66), log(0.46)))
cat(CI)
cat(exp(CI))

```

</details>

**Remark:** skip (e), we will not focus so much on model selection. However, a _top-down-strategy_ (which this exam question points to) will be explained in the end of this module page (week 2)

<span class="question">e) Assume that you want to compare different models with respect to which random effects *and* fixed effects that should be included in the model. Write down a general strategy for doing model evaluation on these kinds of models. </span>

<details><summary>Answer</summary>
1. Start with a model as close to the saturated as possible (all covariates and as many interactions as possible)
2. Find the optimal structure on all random effects (use REML)
3. Find the optimal structure on all fixed effects (use ML)
4. Present the final model with REML estimates
</details>


<span class="question">f) Below is a (horizontal version of a) catepillarplot of the random effects with confidence intervals, and a density plot for the random effects and a QQ-plot for the random effects. Explain what you see.</span>

```{r}
library(ggplot2)
library(sjPlot)
plot_model(fit1,y.offset=0.5,type="re")
ggplot() + geom_density(aes(x = ranef(fit1)$haul[[1]])) + theme_minimal()+
  labs(x = "x", y = "y", title = "Density")
gg <- plot_model(fit1, type = "diag", prnt.plot = FALSE, title = "Quantile plot", geom.size = 1)
gg[[2]]
```

<details><summary>Answer</summary>
From the plots we see that it is reasonable to say that the random intercepts follow a normal distribution with mean 0 and standard deviation around 0.05 (estimated value is 0.05505):

```{r}
ggplot(data = data.frame(x = ranef(fit1)$haul[[1]]), aes(x = x)) + geom_density() + theme_minimal() + 
  labs(x = "x", y = "y", title = "Density") +
  stat_function(fun = dnorm, args = list(mean = 0, sd = attr(VarCorr(fit1)$haul, "stddev")), col = "red")
```


</details>

---

## Exercise 2: Taken from UiO, STK3100, 2013, problem 2

The data in this problem is part of a longitudinal study of income in the US, the Panel Study of Income Dynamics, begun in 1968. The subset consists of 42 heads of household who were aged 25-39 in 1968. The variables included are

* annual nominal income, which is the response variable
* `age`, age in 1968
* `cyear`, coded as -10 in 1968, 0 in 1978 and 10 in 1988
* `educ`, years of education in 1968
* `sex`, M = male, F = female
* `person`, containing information on person id

Below is an excerpt from the output from fitting a linear mixed model (LMM) with the function `lmer` in `R`. The log income (`lincm`) is the response,

$$\texttt{lincm}_{ij} = \beta_0 + \beta_1 \texttt{age}_i + \beta_2 \texttt{cyear}_{ij} + \beta_3 \texttt{educ}_i + \beta_4 \texttt{sex}_i + \gamma_{0i} + \varepsilon_{ij}$$
$$i = 1, \dots, 42, \ j = 1968, 1978, 1988$$

where $\gamma_{0i}$ represent the random effects.

```{r, echo = FALSE}

library(faraway)
library(lme4)

data(psid, package = "faraway")
psid$cyear <- psid$year - 78
psid2 <- subset(psid, cyear %in% c(-10, 0, 10))
psid2 <- psid2[psid2$person %in% names(which(table(psid2$person) >= 3)),]
psid2$lincm <- log(psid2$income)

fit2 <- lmer(lincm ~ age + cyear + educ + sex + (1 | person), data = psid2)

summary(fit2)

```


<span class="question">a) Formulate the model on matrix form and explain the meaning of the different parts. State the model assumptions. </span>

<details><summary>Answer</summary>
On a vector form for responses $\textbf{Y}_i^T = (Y_{i1}, Y_{i2}, Y_{i3})^T$:

$$
\begin{pmatrix} Y_{i1} \\ Y_{i2} \\ Y_{i3} \end{pmatrix} 
= 
\begin{pmatrix} 
  1 & \texttt{age}_{i1} & \texttt{cyear}_i & \texttt{sex}_i \\
  1 & \texttt{age}_{i2} & \texttt{cyear}_i & \texttt{sex}_i \\
  1 & \texttt{age}_{i3} & \texttt{cyear}_i & \texttt{sex}_i \\
\end{pmatrix}
\begin{pmatrix} \beta_0 \\ \beta_1 \\ \beta_2 \\ \beta_3 \\ \beta_4 \end{pmatrix}
+
\begin{pmatrix} 1 \\ 1 \\ 1 \end{pmatrix} \mathbf{\gamma}_{0i}
+
\begin{pmatrix} \varepsilon_{i1} \\ \varepsilon_{i2} \\ \varepsilon_{i3} \end{pmatrix}
$$

Which has the form of a linear mixed model $\mathbf{Y}_i = X_i\beta + U_i\mathbf{\gamma}_{0i} + \varepsilon_i$, where the columns of the $n_i \times p$-matrix $U_i$ is a subset of the columns of the $n_i \times p$-matrix $X_i$. Here, $\pmb{\gamma}_{0i}$ and $\varepsilon_i$ are independent, and

$$\pmb{\gamma}_{0i} \sim N_{q+1}(0, Q) \text{ and } \varepsilon_i = \begin{pmatrix} \varepsilon_{i1} \\ \varepsilon_{i2} \\ \varepsilon_{i3} \end{pmatrix} \sim N_{n_i}(0, \Sigma_i), \ i = 1, \dots, N $$

In our case, $N = 42$, $n_i = 3$, $p = 4$, $q = 0$ (since we only have random intercept), and $\Sigma_i = \sigma^2 I_3$ where $I_3$ is the $3 \times 3$ identity matrix. Here $\beta_0, \dots, \beta_4$ describe the fixed effects and $\textbf{b}_i$ the random effects.

Assumptions: $\pmb{\gamma}_{0i}$ and $\varepsilon_{i}$ must be independent, and we need a linear relationshop between $y$ and $x$. In general, we need things to be like the model says (e.g. Gaussian errors).
</details>

<span class="question">b) Determine an approximate 95 % interval for the coefficient of `cyear`. Do you think the nominal income has been constant in the period covered by the survey?</span>

<details><summary>Answer</summary>

```{r}

library(faraway)
library(lme4)

data(psid, package = "faraway")
psid$cyear <- psid$year - 78
psid2 <- subset(psid, cyear %in% c(-10, 0, 10))
psid2 <- psid2[psid2$person %in% names(which(table(psid2$person) >= 3)),]
psid2$lincm <- log(psid2$income)

fit2 <- lmer(lincm ~ age + cyear + educ + sex + (1 | person), data = psid2)

summary(fit2)

critval <- qnorm(0.025, lower.tail = FALSE)
interval <- summary(fit2)$coefficients[3,1] + critval * summary(fit2)$coefficients[3,2]*c(-1, 1)
Q <- c(summary(fit2)$varcor$person) # q = 1, so this is not a matrix (a 1x1-matrix)
sigma_i <- sigma(fit2)^2*diag(3)

```


From the output we see that the coefficient of `cyear`, $\beta_3$ is estimated to be $\hat{\beta}_3 = 0.08416257$ with estimated standard error 0.008188946. Hence, an approximate interval has limits $0.08416257 \pm 1.96 \times 0.008188946$, so the interval is $(0.06811253, 0.10021261)$. The interval does not contain zero, so a test for constant nominal income would be rejected.

</details>

<span class="question">e) Use the values in the `R`-output to calculate the estimated covariance matrix for the response $(Y_{i1}, Y_{i2}, Y_{i3})^T$ by hand. Hint: compound symmetry due to random intercept model. </span>

<details><summary>Answer</summary>
The elements in the covariance matrix are

$$ \text{Cov}(Y_{ij}, Y_{kl}) = \begin{cases} \tau_0^2 + \sigma^2 = \text{Var}(Y_{ij}) & \text{ for } i=k, j=l \\ \tau_0^2 & \text{ for } i=k, j \neq l \\ 0 & \text{ for } i \neq k, j \neq l \end{cases} $$


From d) it follows that $\hat{\Sigma}_{\mathbf{Y}_i} = U_i\hat{Q}U_i^T + \hat{\Sigma}_i$. In this case $Z_i = (1,1,1)^T$ and $\hat{\Sigma}_i = \hat{\sigma}^2 I_3$. From the output we have that $\hat{Q} = \hat{\tau}_0^2 =$ `r Q` $\approx 0.0018$ and $\hat{\sigma} =$ `r sigma_i[1,1]` $\approx 0.563$. Thus we get

$$\hat{\Sigma}_{\mathbf{Y}_i} = \hat{\tau}_0^2 \begin{pmatrix} 1 & 1 & 1 \\ 1 & 1 & 1 \\ 1 & 1 & 1 \end{pmatrix} + \hat{\sigma}^2 \begin{pmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{pmatrix}$$

$$ = \begin{pmatrix} 0.0018+0.563 & 0.0018 & 0.0018 \\ 0.0018 & 0.0018+0.563 & 0.0018 \\ 0.0018 & 0.0018 & 0.0018+0.563 \end{pmatrix} $$

$$ = \begin{pmatrix} 0.5648 & 0.0018 & 0.0018 \\ 0.0018 & 0.5648 & 0.0018 \\ 0.0018 & 0.0018 & 0.5648 \end{pmatrix} $$

(diagonal elements all equal to $\hat{\sigma}^2 + \hat{\tau}_0^2$ and off-diagonal elements all equal to $\hat{\tau}_0^2$).

</details>

**Remark** Skip (c), which will be lectured in week 2.

<span class="question">c) Explain how one can test the simultaneous significance of two fixed effects (e.g. `age` and `educ`) in LMMs.</span>

<details><summary>Answer</summary>
The question involves the fixed effects. One way to do it is to fit two models:

* one full model containing all the fixed effects, and
* one nested model where the two effects `age` and `educ` are omitted, 

and obtain the maximal value of the likelihood at the two models, call them $L_{max, full}$ and $L_{max, nested}$. The likelihood ratio test consists of comparing 

$$-2 \log(L_{max, nested}/L_{max, full}) =-2(l_{max, nested}-l_{max, full})$$

to a $\chi^2$-distribution where the degrees of freedom are the difference of the number of parameters in the two models, i.e., the number of restrictions which is **two** in this case.

It is important that the ordinary maximum likelihood estimates are used. The restricted maximum likelihood method (REML) consists of basing the estimates on a linear transformation of the data. These transformations are different for the full and the nested model and involve unequal reductions of the data. Therefore it does not make sense to compare the REML likelihoods since they are based on different data.

An alternative is to use the Wald test. The approximate distribution of the estimators of the coefficients have covariance matrix 

$$\Sigma_{\hat{\beta}} = \left( \sum_{i = 1}^N (X_i^T \Sigma_i(\hat{\theta})^{-1}X_i) \right)^{-1}$$

where $\Sigma_i(\theta)$ is the covariance matrix of $\mathbf{Y}_i$ and $\theta$ are the parameters that describe this covariance. The Wald statistic for the null hypothesis

$$H_0 : C\beta = d$$

with $C$ a $r \times s$-matrix ($s$ is the number of parameters in the hypothesis), is 

$$(C\hat{\beta}-d)^T(C\hat{\Sigma}_{\hat{\beta}}C^T)^{-1}(C\hat{\beta}-d)$$

and is approximately $\chi^2$-distributed with $s$ degrees of freedom under $H_0$. In our case, $s = 2$, $C = I_2$, and $d = (0,0)^T$.

```{r}

# lrt:

fit_full <- lmer(lincm ~ age + cyear + educ + sex + (1 | person), data = psid2)
fit_nest <- lmer(lincm ~ cyear + sex + (1 | person), data = psid2)

# un-restricted (regular) likelihood
lrt <- -2*(c(logLik(fit_nest)) - c(logLik(fit_full))); lrt
pchisq(lrt, 2, lower.tail = FALSE) # reject H0


# wald:

C <- diag(2)
beta <- as.vector(summary(fit2)$coefficients[c(2,4),1])
sigma <- as.matrix(vcov(fit2)[c(2,4), c(2,4)])
d <- matrix(c(0,0), ncol = 1)

# since d is a zero-vector and C is a diagonal matrix we can ignore them in the calculations
# wald <- t(C%*%beta-d)%*%solve(C%*%sigma%*%t(C))%*%(C%*%beta-d); wald
wald <- t(beta)%*%solve(sigma)%*%beta; wald

pchisq(wald, 2, lower.tail = FALSE) # small p-value, reject H0

```

Note that you do not need any new information to find the Wald statistic, as you have all you need from the `summary`. To calculate this was not expected on the exam.

```{r}

beta # both under "Fixed effects" in the summary
sigma # the squared roots of the diagonal is under "Fixed effects" in the summary (standard deviations)
sqrt(diag(sigma))

# the covariances are trickier, but can be found by using that 
# the correlation between age and educ is -0.201, found under
# "Correlation of Fixed Effects", and then multiplying with the standard deviations
-0.201*0.015207*0.027582

```

</details>


<span class="question">d) Describe how the random effects $\gamma_{0i}$ can be predicted/estimated. What information is missing for you to calculate these. </span>

<details><summary>Answer</summary>
The estimates of $\pmb{\gamma}_{0i}$ are based on the conditional expectations E$(\pmb{\gamma}_{0i} | \mathbf{Y}_1, \dots, \mathbf{Y}_N) =$ E$(\pmb{\gamma}_{0i} | \mathbf{Y}_i)$ since $\mathbf{Y}_1, \dots, \mathbf{Y}_N$ are independent and $\pmb{\gamma}_{0i}$ only depends on $\mathbf{Y}_i$. Thd simultaneuous distribution of $\pmb{\gamma}_{0i}$ and $\mathbf{Y}_i$ is a $(q + n_i)$-dimensional multinormal with expectation and covariance matrix

$$\begin{pmatrix} 0 \\ X_i\beta \end{pmatrix} \text{ and } \begin{pmatrix} Q & QZ_i^T \\ Z_iQ & Z_iQZ_i^T + \Sigma_i \end{pmatrix} $$

It then follows from the properties if the multinomial distribution that

$$E(\pmb{\gamma}_{0i} | \mathbf{Y}_i) = 0 + QU_i^T(U_iQU_i^T + \Sigma_i)^{-1}(\mathbf{Y}_i - X_i\beta) $$

Hence $\pmb{\gamma}_{0i}$ is estimated by

$$\hat{Q}Z_i^T(U_i\hat{Q}U_i^T + \hat{\Sigma}_i)^{-1}(\mathbf{Y}_i - X_i\hat{\beta}) $$

where all estimates are the REML estimates.

```{r}

U_i <- matrix(c(1,1,1), ncol = 1)
Q <- c(summary(fit2)$varcor$person) # q = 1, so this is not a matrix (a 1x1-matrix)
beta <- as.matrix(fixef(fit2))
sigma_i <- sigma(fit2)^2*diag(3)

# choosing i = 2
i <- 2
Y_i <- as.matrix(subset(psid2, person == i)$lincm)
X_i <- as.matrix(data.frame(int = c(1,1,1), 
                            subset(psid2, person == i)[,c(1, 7, 2)], 
                            sex = ifelse (psid2$sex[psid2$person == i][1] == "M", c(1,1,1), c(0,0,0))))


Q%*%t(U_i)%*%solve(U_i%*%Q%*%t(U_i) + sigma_i)%*%(Y_i - X_i%*%beta)
# can check and see if the estimate coincides with the estimate from lmer
ranef(fit2)$person[which(row.names(ranef(fit2)$person) == i),] # and it does

```


Need $Y_i$ and $X_i$, everything else can be found (with some calculations) from the summary. 

</details>

In part c) and d) it was not necessary to do any numerical calculations at the exam.

<span class="question">f) New: do c) and d) in `R` (after you have done c) and d) by hand!!!, i.e., do numerical calculations (choose $i = 2$ in d)). The rmd-file of this module contains necessary code to download the dataset, modify it, and fit the full model. Alternatively use `purl` to extract the code.</span>

<details><summary>Answer</summary>
This is a random intercept model. Such models are useful if you want to model correlations between variables from the same group (clustered data) or the same invidual (longitudinal data).
</details>

---


