Select the method of your choice.

Mixed effect models contain both fixed and random effects. The general form of the mixed effect model is:

**y **=** Xβ **+** Z**_{1}**μ**_{1}** **+** Z**_{2}**μ**_{2} + ... + **Z**_{c}**μ**_{c} + **ε**

Term | Description |
---|---|

y | the n x 1 vector of response values |

X | the n x p design matrix for the fixed effects, p ≤ n |

Z_{i} | the n x mdesign matrix for the _{i} i^{th} random effect in the model |

β | a p x 1 vector of unknown parameters |

μ_{i} | an m x 1 vector of independent variables from N(0, _{i}σ^{2})_{i} |

ε | an n x 1 vector of independent variables from N(0, σ^{2})_{i} |

c | the number of random effects in the model |

Stability studies fits two models with a random batch factor. The largest model contains time, the random batch factor, and the random interaction between time and batch.

**y **=** Xβ **+** Z**_{1}**μ**_{1}** **+** Z**_{2}**μ**_{2} + **ε**

The smaller model contains time and the random batch factor.

**y **=** Xβ **+** Z**_{1}**μ**_{1}** **+ **ε**

The general variance-covariance matrix of the response vector, **y**, is:

V(**σ**^{2}) = V(σ^{2}, σ^{2}_{1}, ... , σ^{2}_{c}) = σ^{2}I_{n} + σ^{2}_{1}**Z**_{1}**Z'**_{1} + ... + σ^{2}_{c}**Z**_{c}**Z'**_{c}

where

**σ**^{2} = (σ^{2}, σ^{2}_{1}, ... , σ^{2}_{c})'

σ^{2}, σ^{2}_{1}, ... , σ^{2}_{c} are called variance components.

By factoring from the variance, you can find a representation of H(*θ*), which is in the computation of the log-likelihood of mixed models.

V(**σ**^{2}) = σ^{2}H(*θ*) = σ^{2}**[I**_{n} + *θ*_{1}**Z**_{1}**Z'**_{1} + ... + *θ*_{c}**Z**_{c}**Z'**_{c}]

When batch is a random factor, the unknown parameter estimates come from minimizing twice the negative of the restricted log-likelihood function. The minimization is equivalent to maximizing the restricted log-likelihood function. The function to minimize is:

Term | Description |
---|---|

n | the number of observations |

p | the number of parameters in β, 2 for stability studies |

σ^{2} | the error variance component |

X | the design matrix ––for the fixed terms, the constant and time |

H(θ) | I_{n} + θ_{1}Z_{1}Z'_{1} + ... + θ_{c}Z_{c}Z'_{c} |

I_{n} | the identity matrix with n rows and columns |

θ_{i} | the ratio of the variance of the i^{th} random term over the error variance |

Z_{i} | the n x m matrix of known codings for the _{i}i^{th} random effect in the model |

m_{i} | the number of levels for the i^{th} random effect |

c | the number of random effects in the model |

|H(θ)| | the determinant of H(θ) |

X' | the transpose of X |

H^{-1}(θ) | the inverse of H(θ) |

Box-Cox transformation selects lambda values, as shown below, which minimize the residual sum of squares. The resulting transformation is *Y* ^{λ} when λ ≠ 0 and ln(*Y*) when λ = 0. When λ < 0, Minitab also multiplies the transformed response by −1 to maintain the order from the untransformed response.

Minitab searches for an optimal value between −2 and 2. Values that fall outside of this interval might not result in a better fit.

Here are some common transformations where *Y*′ is the transform of the data *Y*:

Lambda (λ) value | Transformation |
---|---|

λ = 2 | Y′ = Y ^{2} |

λ = 0.5 | Y′ = |

λ = 0 | Y′ = ln(Y ) |

λ = −0.5 | |

λ = −1 | Y′ = −1 / Y |

The model selection determines whether the shelf life depends on batch and whether the effect of time depends on the batch. Minitab considers the following three models in sequence:

- Time + Batch + Batch*Time (unequal slopes and intercepts for batches)
- Time + Batch (equal slopes and unequal intercepts for batches)
- Time (equal slopes and intercepts for batches)

If the Batch*Time interaction is significant, the analysis fits the first model. If the interaction is not significant but the Batch term is significant in the second model, the analysis fits the second model. Otherwise, the analysis fits the third model.

The test for whether to pool batches is slightly different from the test to include batch, although both depend on the chi-square distribution. The formulas for the test statistics and p-values are as follow.

difference = −2L_{2} − (−2L_{1})

p = 0.5 * Prob(χ^{2}_{1} > difference) + 0.5 * Prob(χ^{2}_{2} > difference)

difference = −2L_{3} − (−2L_{2})

p = 0.5 * Prob(χ^{2}_{1} > difference)

Term | Description |
---|---|

L_{a} | the log-likelihood for model a |

p | the p-value for the test |

Prob(χ^{2}_{1}> difference) | the probability that a random variable from a chi-square distribution with 1 degree of freedom is greater than the difference |

Prob(χ^{2}_{2}> difference) | the probability that a random variable from a chi-square distribution with 2 degrees of freedom is greater than the difference |

- Searle, S.R. Casella, G. and McCuloch, C.E. (1992). Variance Components
- West, B.T., Welch, K.B. and Galecki, A.T. (2007). Linear Mixed Models: A Practical Guide Using Statistical Software.
- Chow, S. (2007). Statistical Design and Analysis of Stability Studies.