09-nonlinear_generalized_linear_mixed.Rmd

# Nonlinear and Generalized Linear Mixed Models

-   NLMMs extend the nonlinear model to include both fixed effects and random effects
-   GLMMs extend the generalized linear model to include both fixed effects and random effects.

A nonlinear mixed model has the form of

$$
Y_{ij} = f(\mathbf{x_{ij} , \theta, \alpha_i}) + \epsilon_{ij}
$$

for the j-th response from cluster (or subject) i ($i = 1,...,n$), where

-   $j = 1,...,n_i$
-   $\mathbf{\theta}$ are the fixed effects
-   $\mathbf{\alpha}_i$ are the random effects for cluster i
-   $\mathbf{x}_{ij}$ are the regressors or design variables
-   $f(.)$ is nonlinear mean response function

A GLMM can be written as:

we assume

$$
y_i |\alpha_i \sim \text{indep } f(y_i | \alpha)
$$

and $f(y_i | \mathbf{\alpha})$ is an exponential family distribution,

$$
f(y_i | \alpha) = \exp [\frac{y_i \theta_i - b(\theta_i)}{a(\phi)} - c(y_i, \phi)]
$$

The conditional mean of $y_i$ is related to $\theta_i$

$$
\mu_i = \frac{\partial b(\theta_i)}{\partial \theta_i}
$$

The transformation of this mean will give us the desired linear model to model both the fixed and random effects.

$$
\begin{aligned}
E(y_i |\alpha) &= \mu_i \\
g(\mu_i) &= \mathbf{x_i' \beta + z'_i \alpha}
\end{aligned}
$$

where $g()$ is a known link function and $\mu_i$ is the conditional mean. We can see similarity to [GLM](#generalized-linear-models)

We also have to specify the random effects distribution

$$
\alpha \sim f(\alpha)
$$

which is similar to the specification for mixed models.

Moreover, law of large number applies to fixed effects so that you know it is a normal distribution. But here, you can specify $\alpha$ subjectively.

Hence, we can show NLMM is a special case of the GLMM

$$
\begin{aligned}
\mathbf{Y}_i &= \mathbf{f}(\mathbf{x}_i, \mathbf{\theta, \alpha}_i) + \mathbf{\epsilon}_i \\
\mathbf{Y}_i &= \mathbf{g}^{-1} (\mathbf{x}_i' \beta + \mathbf{z}_i' \mathbf{\alpha}_i) + \mathbf{\epsilon}_i
\end{aligned}
$$

where the inverse link function corresponds to a nonlinear transformation of the fixed and random effects.

**Note**:

-   we can't derive the analytical formulation of the marginal distribution because nonlinear combination of normal variables is not normally distributed, even in the case of additive error ($e_i$) and random effects ($\alpha_i$) are both normal.

**Consequences of having random effects**

The marginal mean of $y_i$ is

$$
E(y_i) = E_\alpha(E(y_i | \alpha)) = E_\alpha (\mu_i) = E(g^{-1}(\mathbf{x_i' \beta + z_i' \alpha}))
$$

Because $g^{-1}()$ is nonlinear, this is the most simplified version we can go for.

In special cases such as log link ($g(\mu) = \log \mu$ or $g^{-1}() = \exp()$) then

$$
E(y_i) = E(\exp(\mathbf{x_i' \beta + z_i' \alpha})) = \exp(\mathbf{x'_i \beta})E(\exp(\mathbf{z}_i'\alpha))
$$

which is the moment generating function of $\alpha$ evaluated at $\mathbf{z}_i$

**Marginal variance** of $y_i$

$$
\begin{aligned}
var(y_i) &= var_\alpha (E(y_i | \alpha)) + E_\alpha (var(y_i | \alpha)) \\
&= var(\mu_i) + E(a(\phi) V(\mu_i)) \\
&= var(g^{-1} (\mathbf{x'_i \beta + z'_i \alpha})) + E(a(\phi)V(g^{-1} (\mathbf{x'_i \beta + z'_i \alpha})))
\end{aligned}
$$

Without specific assumption about $g()$ and/or the conditional distribution of $\mathbf{y}$, this is the most simplified version.

**Marginal covariance of** $\mathbf{y}$

In a linear mixed model, random effects introduce a dependence among observations which share any random effect in common

$$
\begin{aligned}
cov(y_i, y_j) &= cov_{\alpha}(E(y_i | \mathbf{\alpha}),E(y_j | \mathbf{\alpha})) + E_{\alpha}(cov(y_i, y_j | \mathbf{\alpha})) \\
&= cov(\mu_i, \mu_j) + E(0) \\
&= cov(g^{-1}(\mathbf{x}_i' \beta + \mathbf{z}_i' \mathbf{\alpha}), g^{-1}(\mathbf{x}'_j \beta + \mathbf{z}_j' \mathbf{\alpha}))
\end{aligned}
$$

-   Important: conditioning to induce the covariability

Example:

Repeated measurements on the subjects.

Let $y_{ij}$ be the j-th count taken on the $i$-th subject.

then, the model is $y_{ij} | \mathbf{\alpha} \sim \text{indep } Pois(\mu_{ij})$. Here

$$
\log(\mu_{ij}) = \mathbf{x}_{ij}' \beta + \alpha_i 
$$

where $\alpha_i \sim iid N(0,\sigma^2_{\alpha})$

which is a log-link with a random patient effect.

## Estimation

In linear mixed models, the marginal likelihood for $\mathbf{y}$ is the integration of the random effects from the hierarchical formulation

$$
f(\mathbf{y}) = \int f(\mathbf{y}| \alpha) f(\alpha) d \alpha
$$

For linear mixed models, we assumed that the 2 component distributions were Gaussian with linear relationships, which implied the marginal distribution was also linear and Gaussian and allows us to solve this integral analytically.

On the other hand, GLMMs, the distribution for $f(\mathbf{y} | \alpha)$ is not Gaussian in general, and for NLMMs, the functional form between the mean response and the random (and fixed) effects is nonlinear. In both cases, we can't perform the integral analytically, which means we have to solve it

-   [numerically](#estimation-by-numerical-integration) and/or

-   [linearize the inverse link function](#estimation-by-linearization).

### Estimation by Numerical Integration {#estimation-by-numerical-integration}

The marginal likelihood is

$$
L(\beta; \mathbf{y}) = \int f(\mathbf{y} | \alpha) f(\alpha) d \alpha
$$

Estimation fo the fixed effects requires $\frac{\partial l}{\partial \beta}$, where $l$ is the log-likelihood

One way to obtain the marginal inference is to numerically integrate out the random effects through

-   numerical quadrature

-   Laplace approximation

-   Monte Carlo methods

When the dimension of $\mathbf{\alpha}$ is relatively low, this is easy. But when the dimension of $\alpha$ is high, additional approximation is required.

### Estimation by Linearization {#estimation-by-linearization}

Idea: Linearized version of the response (known as working response, or pseudo-response) called $\tilde{y}_i$ and then the conditional mean is

$$
E(\tilde{y}_i | \alpha) = \mathbf{x}_i' \beta + \mathbf{z}_i' \alpha
$$

and also estimate $var(\tilde{y}_i | \alpha)$. then, apply [Linear Mixed Models] estimation as usual.

The difference is only in how the linearization is done (i.e., how to expand $f(\mathbf{x, \theta, \alpha})$ or the inverse link function

#### Penalized quasi-likelihood

(PQL)

This is the more popular method

$$
\tilde{y}_i^{(k)} = \hat{\eta}_i^{(k-1)} + ( y_i - \hat{\mu}_i^{(k-1)})\frac{d \eta}{d \mu}| \hat{\eta}_i^{(k-1)}
$$

where

-   $\eta_i = g(\mu_i)$ is the linear predictor

-   $k$ = iteration of the optimization algorithm

The algorithm updates $\tilde{y}_i$ after each linear mixed model fit using $E(\tilde{y}_i | \alpha)$ and $var(\tilde{y}_i | \alpha)$

Comments:

-   Easy to implement

-   Inference is only asymptotically correct due to the linearizaton

-   Biased estimates are likely for binomial response with small groups and worst for Bernoulli response. Similarly for Poisson models with small counts. [@faraway2016extending]

-   Hypothesis testing and confidence intervals also have problems.

#### Generalized Estimating Equations

(GEE)

Let a marginal generalized linear model for the mean of y as a function of the predictors, which means we linearize the mean response function and assume a dependent error structure

Example\
Binary data:

$$
logit (E(\mathbf{y})) = \mathbf{X} \beta
$$

If we assume a "working covariance matrix", $\mathbf{V}$ the the elements of $\mathbf{y}$, then the maximum likelihood equations for estimating $\beta$ is

$$
\mathbf{X'V^{-1}y} = \mathbf{X'V^{-1}} E(\mathbf{y})
$$

If $\mathbf{V}$ is correct, then unbiased estimating equations

We typically define $\mathbf{V} = \mathbf{I}$. Solutions to unbiased estimating equation give consistent estimators.

In practice, we assume a covariance structure, and then do a logistic regression, and calculate its large sample variance

Let $y_{ij} , j = 1,..,n_i, i = 1,..,K$ be the j-th measurement on the $i$-th subject.

$$
\mathbf{y}_i = 
\left(
\begin{array}
{c}
y_{i1} \\
. \\
y_{in_i}
\end{array}
\right)
$$

with mean

$$
\mathbf{\mu}_i =
\left(
\begin{array}
{c}
\mu_{i1} \\
. \\
\mu_{in_i}
\end{array}
\right)
$$

and

$$
\mathbf{x}_{ij} = 
\left(
\begin{array}
{c}
X_{ij1} \\
. \\
X_{ijp}
\end{array}
\right)
$$

Let $\mathbf{V}_i = cov(\mathbf{y}_i)$, then based on[@liang1986longitudinal] GEE estimates for $\beta$ can be obtained from solving the equation:

$$
S(\beta) = \sum_{i=1}^K \frac{\partial \mathbf{\mu}_i'}{\partial \beta} \mathbf{V}^{-1}(\mathbf{y}_i - \mathbf{\mu}_i) = 0
$$

Let $\mathbf{R}_i (\mathbf{c})$ be an $n_i \times n_i$ "working" correlation matrix specified up to some parameters $\mathbf{c}$. Then, $\mathbf{V}_i = a(\phi) \mathbf{B}_i^{1/2}\mathbf{R}(\mathbf{c}) \mathbf{B}_i^{1/2}$, where $\mathbf{B}_i$ is an $n_i \times n_i$ diagonal matrix with $V(\mu_{ij})$ on the j-th diagonal

If $\mathbf{R}(\mathbf{c})$ is the true correlation matrix of $\mathbf{y}_i$, then $\mathbf{V}_i$ is the true covariance matrix

The working correlation matrix must be estimated iteratively by a fitting algorithm:

1.  Compute the initial estimate of $\beta$ (using GLM under the independence assumption)

2.  Compute the working correlation matrix $\mathbf{R}$ based upon studentized residuals

3.  Compute the estimate covariance $\hat{\mathbf{V}}_i$

4.  Update $\beta$ according to

    $$
    \beta_{r+1} = \beta_r + (\sum_{i=1}^K \frac{\partial \mathbf{\mu}'_i}{\partial \beta} \hat{\mathbf{V}}_i^{-1} \frac{\partial \mathbf{\mu}_i}{\partial \beta})
    $$

5.  Iterate until the algorithm converges

Note: Inference based on likelihoods is not appropriate because this is not a likelihood estimator

### Estimation by Bayesian Hierarchical Models

Bayesian Estimation

$$
f(\mathbf{\alpha}, \mathbf{\beta} | \mathbf{y}) \propto f(\mathbf{y} | \mathbf{\alpha}, \mathbf{\beta}) f(\mathbf{\alpha})f(\mathbf{\beta})
$$

Numerical techniques (e.g., MCMC) can be used to find posterior distribution. This method is best in terms of not having to make simplifying approximation and fully accounting for uncertainty in estimation and prediction, but it could be complex, time-consuming, and computationally intensive.

Implementation Issues:

-   No valid joint distribution can be constructed from the given conditional model and random parameters

-   The mean/ variance relationship and the random effects lead to constraints on the marginal covariance model

-   Difficult to fit computationally

2 types of estimation approaches:

1.  Approximate the objective function (marginal likelihood) through integral approximation

    1.  Laplace methods

    2.  Quadrature methods

    3.  Monte Carlo integration

2.  Approximate the model (based on Taylor series linearization)

Packages in R

-   GLMM: `MASS:glmmPQL` `lme4::glmer` `glmmTMB`

-   NLMM: `nlme::nlme`; `lme4::nlmer` `brms::brm`

-   Bayesian: `MCMCglmm` ; `brms:brm`

Example: Non-Gaussian Repeated measurements

-   When the data are Gaussian, then [Linear Mixed Models]

-   When the data are non-Gaussian, then [Nonlinear and Generalized Linear Mixed Models]

## Application

### Binomial (CBPP Data)

```{r}
data(cbpp,package = "lme4")
head(cbpp)
```

PQL

Pro:

-   Linearizes the response to have a pseudo-response as the mean response (like LMM)

-   computationally efficient

Cons:

-   biased for binary, Poisson data with small counts

-   random effects have to be interpreted on the link scale

-   can't interpret AIC/BIC value

```{r}
library(MASS)
pql_cbpp <-
    glmmPQL(
        cbind(incidence, size - incidence) ~ period,
        random  = ~ 1 | herd,
        data    = cbpp,
        family  = binomial(link = "logit"),
        verbose = F
    )
summary(pql_cbpp)
```

```{r}
exp(0.556)
```

is how the herd specific outcome odds varies.

We can interpret the fixed effect coefficients just like in GLM. Because we use logit link function here, we can say that the log odds of the probability of having a case in period 2 is -1.016 less than period 1 (baseline).

```{r}
summary(pql_cbpp)$tTable
```

Numerical Integration

Pro:

-   more accurate

Con:

-   computationally expensive

-   won't work for complex models.

```{r}
library(lme4)
numint_cbpp <-
    glmer(
        cbind(incidence, size - incidence) ~ 
            period + (1 | herd),
        data = cbpp,
        family = binomial(link = "logit")
    )
summary(numint_cbpp)
```

For small data set, the difference between two approaches are minimal

```{r}
library(rbenchmark)
benchmark(
    "MASS" = {
        pql_cbpp <-
            glmmPQL(
                cbind(incidence, size - incidence) ~ period,
                random = ~ 1 | herd,
                data = cbpp,
                family = binomial(link = "logit"),
                verbose = F
            )
    },
    "lme4" = {
        glmer(
            cbind(incidence, size - incidence) ~ period + (1 | herd),
            data = cbpp,
            family = binomial(link = "logit")
        )
    },
    replications = 50,
    columns = c("test", "replications", "elapsed", "relative"),
    order = "relative"
)

```

In numerical integration, we can set `nAGQ > 1` to switch the method of likelihood evaluation, which might increase accuracy

```{r}
library(lme4)
numint_cbpp_GH <-
    glmer(
        cbind(incidence, size - incidence) ~ period + (1 | herd),
        data = cbpp,
        family = binomial(link = "logit"),
        nAGQ = 20
    )
summary(numint_cbpp_GH)$coefficients[, 1] - 
    summary(numint_cbpp)$coefficients[, 1]

```

Bayesian approach to GLMMs

-   assume the fixed effects parameters have distribution

-   can handle models with intractable result under traditional methods

-   computationally expensive

```{r}
library(MCMCglmm)
Bayes_cbpp <-
    MCMCglmm(
        cbind(incidence, size - incidence) ~ period,
        random  = ~ herd,
        data    = cbpp,
        family  = "multinomial2",
        verbose = FALSE
    )
summary(Bayes_cbpp)
```

-   `MCMCglmm` fits a residual variance component (useful with dispersion issues)

```{r}
# explains less variability
apply(Bayes_cbpp$VCV,2,sd)
```

```{r}
summary(Bayes_cbpp)$solutions
```

interpret Bayesian "credible intervals" similarly to confidence intervals

Make sure you make post-hoc diagnoses

```{r}
library(lattice)
xyplot(as.mcmc(Bayes_cbpp$Sol), layout = c(2, 2))
```

There is no trend, well-mixed

```{r}
xyplot(as.mcmc(Bayes_cbpp$VCV),layout=c(2,1))
```

For the herd variable, a lot of them are 0, which suggests problem. To fix the instability in the herd effect sampling, we can either

-   modify the prior distribution on the herd variation

-   increases the number of iteration

```{r}
library(MCMCglmm)
Bayes_cbpp2 <-
    MCMCglmm(
        cbind(incidence, size - incidence) ~ period,
        random = ~ herd,
        data   = cbpp,
        family = "multinomial2",
        nitt   = 20000,
        burnin = 10000,
        prior  = list(G = list(list(
            V  = 1, nu = .1
        ))),
        verbose = FALSE
    )
xyplot(as.mcmc(Bayes_cbpp2$VCV), layout = c(2, 1))

```

To change the shape of priors, in `MCMCglmm` use:

-   `V` controls for the location of the distribution (default = 1)

-   `nu` controls for the concentration around V (default = 0)

### Count (Owl Data)

```{r}
library(glmmTMB)
library(dplyr)
data(Owls, package = "glmmTMB")
Owls <- Owls %>% 
    rename(Ncalls = SiblingNegotiation)
```

In a typical Poisson model, $\lambda$ (Poisson mean), is model as $\log(\lambda) = \mathbf{x'\beta}$ But if the response is the rate (e.g., counts per BroodSize), we could model it as $\log(\lambda / b) = \mathbf{x'\beta}$ , equivalently $\log(\lambda) = \log(b) + \mathbf{x'\beta}$ where $b$ is BroodSize. Hence, we "offset" the mean by the log of this variable.

```{r}
owls_glmer <-
    glmer(
        Ncalls ~ offset(log(BroodSize)) 
        + FoodTreatment * SexParent +
            (1 | Nest),
        family = poisson,
        data = Owls
    )
summary(owls_glmer)

```

-   nest explains a relatively large proportion of the variability (its standard deviation is larger than some coefficients)

-   the model fit isn't great (deviance of 5202 on 594 df)

```{r}
# Negative binomial model
owls_glmerNB <-
    glmer.nb(Ncalls ~ offset(log(BroodSize)) 
             + FoodTreatment * SexParent
             + (1 | Nest), data = Owls)

c(Deviance = round(summary(owls_glmerNB)$AICtab["deviance"], 3),
  df = summary(owls_glmerNB)$AICtab["df.resid"])

```

There is an improvement using negative binomial considering over-dispersion

```{r}
hist(Owls$Ncalls,breaks=30)
```

To account for too many 0s in these data, we can use zero-inflated Poisson (ZIP) model.

-   `glmmTMB` can handle ZIP GLMMs since it adds automatic differentiation to existing estimation strategies.

```{r}
library(glmmTMB)
owls_glmm <-
    glmmTMB(
        Ncalls ~ FoodTreatment * SexParent + offset(log(BroodSize)) +
            (1 | Nest),
        ziformula =  ~ 0,
        family = nbinom2(link = "log"),
        data = Owls
    )
owls_glmm_zi <-
    glmmTMB(
        Ncalls ~ FoodTreatment * SexParent + offset(log(BroodSize)) +
            (1 | Nest),
        ziformula =  ~ 1,
        family = nbinom2(link = "log"),
        data = Owls
    )
# Scale Arrival time to use as a covariate for zero-inflation parameter
Owls$ArrivalTime <- scale(Owls$ArrivalTime)
owls_glmm_zi_cov <- glmmTMB(
    Ncalls ~ FoodTreatment * SexParent +
        offset(log(BroodSize)) +
        (1 | Nest),
    ziformula =  ~ ArrivalTime,
    family = nbinom2(link = "log"),
    data = Owls
)
as.matrix(anova(owls_glmm, owls_glmm_zi))
as.matrix(anova(owls_glmm_zi, owls_glmm_zi_cov))
summary(owls_glmm_zi_cov)
```

We can see ZIP GLMM with an arrival time covariate on the zero is best.

-   arrival time has a positive effect on observing a nonzero number of calls

-   interactions are non significant, the food treatment is significant (fewer calls after eating)

-   nest variability is large in magnitude (without this, the parameter estimates change)

### Binomial

```{r}
library(agridat)
library(ggplot2)
library(lme4)
library(spaMM)
data(gotway.hessianfly)
dat <- gotway.hessianfly
dat$prop <- dat$y / dat$n

ggplot(dat, aes(x = lat, y = long, fill = prop)) +
    geom_tile() +
    scale_fill_gradient(low = 'white', high = 'black') +
    geom_text(aes(label = gen, color = block)) +
    ggtitle('Gotway Hessian Fly')
```

-   Fixed effects ($\beta$) = genotype

-   Random effects ($\alpha$) = block

```{r}
flymodel <-
    glmer(
        cbind(y, n - y) ~ gen + (1 | block),
        data   = dat,
        family = binomial,
        nAGQ   = 5
    )
summary(flymodel)
```

Equivalently, we can use `MCMCglmm` , for a Bayesian approach

```{r}
library(coda)
Bayes_flymodel <- MCMCglmm(
    cbind(y, n - y) ~ gen ,
    random  = ~ block,
    data    = dat,
    family  = "multinomial2",
    verbose = FALSE
)
plot(Bayes_flymodel$Sol[, 1], main = dimnames(Bayes_flymodel$Sol)[[2]][1])

```

```{r}
autocorr.plot(Bayes_flymodel$Sol[, 1], main = dimnames(Bayes_flymodel$Sol)[[2]][1])
```

### Example from [@Schabenberger_2001] section 8.4.1

```{r}
dat2 <- read.table("images/YellowPoplarData_r.txt")
names(dat2) <- c('tn', 'k', 'dbh', 'totht',
                 'dob', 'ht', 'maxd', 'cumv')
dat2$t <- dat2$dob / dat2$dbh
dat2$r <- 1 - dat2$dob / dat2$totht
```

The cumulative volume relates to the complementary diameter (subplots were created based on total tree height)

```{r}
library(ggplot2)
library(dplyr)
dat2 <- dat2 %>% group_by(tn) %>% mutate(
    z = case_when(
        totht < 74 & totht >= 0 ~ 'a: 0-74ft',
        totht < 88 & totht >= 74 ~ 'b: 74-88',
        totht < 95 & totht >= 88 ~ 'c: 88-95',
        totht < 99 & totht >= 95 ~ 'd: 95-99',
        totht < 104 & totht >= 99 ~ 'e: 99-104',
        totht < 109 & totht >= 104 ~ 'f: 104-109',
        totht < 115 & totht >= 109 ~ 'g: 109-115',
        totht < 120 & totht >= 115 ~ 'h: 115-120',
        totht < 140 & totht >= 120 ~ 'i: 120-150',
    )
)
ggplot(dat2, aes(x = r, y = cumv)) + 
    geom_point(size = 0.5) + 
    facet_wrap(vars(z))
```

The proposed non-linear model:

$$
V_{id_j} = (\beta_0 + (\beta_1 + b_{1i})\frac{D^2_i H_i}{1000})(\exp[-(\beta_2 + b_{2i})t_{ij} \exp(\beta_3 t_{ij})]) + e_{ij}
$$

where

-   $b_{1i}, b_{2i}$ are random effects

-   $e_{ij}$ are random errors

```{r}
library(nlme)
tmp <-
    nlme(
        cumv ~ (b0 + (b1 + u1) *
                    (dbh * dbh * totht / 1000)) *
            (exp(-(b2 + u2) * (t / 1000) * exp(b3 * t))), 
        data = dat2,
        fixed = b0 + b1 + b2 + b3 ~ 1,
        # 1 on the right hand side of the formula indicates 
        # a single fixed effects for the corresponding parameters
        random = list(pdDiag(u1 + u2 ~ 1)),
        #uncorrelated random effects
        groups = ~ tn,
        #group on trees so each tree w/ have u1 and u2
        start = list(fixed = c(
            b0 = 0.25,
            b1 = 2.3,
            b2 = 2.87,
            b3 = 6.7
        ))
    )
summary(tmp)
nlme::intervals(tmp)
```

-   Little different from the book because of different implementation of nonlinear mixed models.

```{r}
library(cowplot)
nlmmfn <- function(fixed,rand,dbh,totht,t){
  b0 <- fixed[1]
  b1 <- fixed[2]
  b2 <- fixed[3]
  b3 <- fixed[4]
  u1 <- rand[1]
  u2 <- rand[2]
  #just made so we can predict w/o random effects
  return((b0+(b1+u1)*(dbh*dbh*totht/1000))*(exp(-(b2+u2)*(t/1000)*exp(b3*t))))
}



#Tree 1
pred1 <- data.frame(seq(1, 24, length.out = 100))
names(pred1) <- 'dob'
pred1$tn <- 1
pred1$dbh <- unique(dat2[dat2$tn == 1, ]$dbh)
pred1$t <- pred1$dob / pred1$dbh
pred1$totht <- unique(dat2[dat2$tn == 1, ]$totht)
pred1$r <- 1 - pred1$dob / pred1$totht


pred1$test <- predict(tmp, pred1)
pred1$testno <-
    nlmmfn(
        fixed = tmp$coefficients$fixed,
        rand = c(0, 0),
        pred1$dbh,
        pred1$totht,
        pred1$t
    )

p1 <-
    ggplot(pred1) + 
    geom_line(aes(x = r, y = test, color = 'with random')) +
    geom_line(aes(x = r, y = testno, color = 'No random')) + 
    labs(colour = "") + 
    geom_point(data = dat2[dat2$tn == 1, ], aes(x = r, y = cumv)) + 
    ggtitle('Tree 1') + theme(legend.position = "none")


#Tree 151
pred151        <- data.frame(seq(1, 21, length.out = 100))
names(pred151) <- 'dob'
pred151$tn     <- 151
pred151$dbh    <- unique(dat2[dat2$tn == 151, ]$dbh)
pred151$t      <- pred151$dob / pred151$dbh
pred151$totht  <- unique(dat2[dat2$tn == 151, ]$totht)
pred151$r      <- 1 - pred151$dob / pred151$totht


pred151$test <- predict(tmp, pred151)
pred151$testno <-
    nlmmfn(
        fixed = tmp$coefficients$fixed,
        rand = c(0, 0),
        pred151$dbh,
        pred151$totht,
        pred151$t
    )

p2 <-
    ggplot(pred151) + 
    geom_line(aes(x = r, y = test, color = 'with random')) +
    geom_line(aes(x = r, y = testno, color = 'No random')) + 
    labs(colour = "") + 
    geom_point(data = dat2[dat2$tn == 151,], aes(x = r, y = cumv)) + 
    ggtitle('Tree 151') + 
    theme(legend.position = "none")


#Tree 279
pred279        <- data.frame(seq(1, 9, length.out = 100))
names(pred279) <- 'dob'
pred279$tn     <- 279
pred279$dbh    <- unique(dat2[dat2$tn == 279, ]$dbh)
pred279$t      <- pred279$dob / pred279$dbh
pred279$totht  <- unique(dat2[dat2$tn == 279, ]$totht)
pred279$r      <- 1 - pred279$dob / pred279$totht


pred279$test <- predict(tmp, pred279)
pred279$testno <-
    nlmmfn(
        fixed = tmp$coefficients$fixed,
        rand = c(0, 0),
        pred279$dbh,
        pred279$totht,
        pred279$t
    )

p3 <-
    ggplot(pred279) + 
    geom_line(aes(x = r, y = test, color = 'with random')) +
    geom_line(aes(x = r, y = testno, color = 'No random')) + 
    labs(colour = "") + 
    geom_point(data = dat2[dat2$tn == 279, ], aes(x = r, y = cumv)) + 
    ggtitle('Tree 279') + 
    theme(legend.position = "none")

plot_grid(p1, p2, p3)
```

red line = predicted observations based on the common fixed effects

teal line = tree-specific predictions with random effects

## Summary

::: {style="text-align:center"}
![](images/umbrella_of_models.PNG){width="450" height="300"}
:::