---
title: "Check perfect model"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Check perfect model}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include=FALSE}
old <- options(digits = 3)
knitr::opts_chunk$set(echo = TRUE)
```

## Objective

In general, the broken stick model smoothes the observed growth trajectory. What happens of all observations are already aligned to the break ages? Does the model perfectly represent the data? Is the covariance matrix of the random effects ($\Omega)$ equal to the covariance between the measurements? Is $\sigma^2$ equal to zero?

## Data generation

We adapt code from http://www.davekleinschmidt.com/sst-mixed-effects-simulation/simulations_slides.pdf to generate test data:

```{r}
library("plyr")
library("mvtnorm")
make_data_generator <- function(resid_var = 1,
                                ranef_covar = diag(c(1, 1)), n = 100
                                ) {
  ni <- nrow(ranef_covar)
  generate_data <- function() {
    # sample data set under mixed effects model with random slope/intercepts 
    simulated_data <- rdply(n, {
      b <- t(rmvnorm(n = 1, sigma = ranef_covar))
      epsilon <- rnorm(n = length(b), mean = 0, sd = sqrt(resid_var))
      b + epsilon
    })
  data.frame(
    subject = rep(1:n, each = ni),
    age = rep(1:ni, n),
    simulated_data)
  }
}
```

We choose between the perfect situation where $\sigma^2 = 0$ and the noisy case $\sigma^2 = 1$ and where the ages align perfectly.

```{r}
resid_var <- 0
resid_var <- 1
set.seed(77711)
covar <- matrix(c(1, 0.7, 0.5, 0.3,
                  0.7, 1, 0.8, 0.5,
                  0.5, 0.8, 1, 0.6,
                  0.3, 0.5, 0.6, 1), nrow = 4)
gen_dat <- make_data_generator(n = 10000, 
                               ranef_covar = covar,
                               resid_var = resid_var)
data <- gen_dat()
head(data)
```

We wish to reproduce the correlation matrix among the $y$'s from the mixed model estimates. The **target correlation matrix** is:

```{r}
library("tidyr")
library("dplyr")
d <- as_tibble(data[,-3])
broad <- t(spread(d, subject, X1))[-1,]
cor(broad)
```

## Fit model

Fit broken stick model, with knots specified at ages `1:4`.

```{r}
library("brokenstick")
knots <- 1:3
boundary <- c(1, 4)
fit <- brokenstick(X1 ~ age | subject, data, 
                   knots = knots, boundary = boundary,
                   method = "lmer")
omega <- get_omega(fit, hide = "no")
beta <- coef(fit, hide = "no")
sigma2 <- fit$sigma2
round(beta, 2)
round(sigma2, 4)

# correlation random effects
round(covar, 3)
round(omega, 2)

# covariances measured data
round(omega + diag(sigma2, 4), 3)
round(cov(broad), 3)

# convert to time-to-time correlation matrix
round(cov2cor(omega + diag(sigma2, 4)), 3)
round(cor(broad), 3)

z <- predict(fit, x = "knots", include_data = FALSE, shape = "wide")[, -1]
# off-diagonal elements of covariance of broken stick estimates approach correlation
# not enough variance in the diagonal because of smoothing
cov(z)
# correlations of broken stick estimates are inflated because of smoothing
cor(z)
```

## Conclusions

1. If $\sigma^2=0$, then the off-diagonal elements of $\Omega$ reproduce the correlations among the $y$'s. The estimate of $\sigma^2$ is too high (about 0.13 instead of 0).
2. If $\sigma^2 > 0$, then $\hat C = \Omega + \hat\sigma^2 I(n_i)$ reproduces the sample covariance matrix between $y$'s exactly.
3. `cov2cor(hatC)` reproduces the sample time-to-time correlation matrix.
4. Conclusion: The unbiased estimate of the time-to-time correlation matrix among the (unobserved) measurements at the knots:

```{r}
cov <- get_omega(fit)
chat <- cov + diag(fit$sigma2, nrow(cov))
r <- cov2cor(chat)
r
```

```{r echo = FALSE}
options(old)
```