Modeling formulation, fitting of model and reduction of model

During this section multiple models are considered. The models included fixed effect, general and generalized random effect models. These are listed in the table below:

The best model would end up being chosen based on the BIC and AIC likelihood criterion where the lowest numbers would decide which model is the best. This is then further compared to the goodness of fit. Note that all the less successful models along with their diagnostics can be seen in the modeling appendix.

First the general model reduction and diagnostic workflow is presented. The method of forward selection was utilized. This is started by creating a simple main effect model and gradually building up with higher order terms. The reduction of the model is done one step at a time utilizing the partial partition ANOVA (type II ANOVA test) for the main effects, and afterwards utilizing a likelihood ratio test along with AIC selection for the higher order terms. These was chosen since neither of the methods relies on the design being completely balanced (1). The methods are shown in detail in the likelihood likelihood appendix.

Fixed gaussian effect model formulation, reduction and diagnostic

The fixed gaussian model can be formulated as previously stated in the DOE section, in the following general matrix nation: \begin{equation} \textbf{Y}=X\boldsymbol{\beta}+\boldsymbol{\varepsilon},\\ \boldsymbol{\varepsilon} \sim N(0,\sigma^{2}\textbf{I}) \end{equation} For the response variable $\textbf{Y}$, it is assumed that the $\textbf{Y} \sim N(\mu,\sigma^{2})$.
The design matrix of the experiments is denoted by $\textbf{X}$ and $\beta$ is a vector of regression coefficients spanning between $\beta_{0}...\beta_{n}$. The $\varepsilon$ is the random residual component which is assumed to be i.i.d and follow a normal distribution with an equally weighted $\sigma^{2}$ around a mean of 0.

The full main effect models of the replicated and unreplicated fixed gaussian linear models can be seen below:

Pilot model

\begin{equation} Y_{i}=\mu+\alpha(bake.time_{i})+\beta(bake.temp_{i})\\+\gamma(mix_{i})+\delta(incu.time_{i})+\zeta(sub_{i})\\+eta(W_{i})+ \varepsilon_{i} \ ,where \ \varepsilon_{i} \sim N(0,\sigma^{2}) \end{equation} The model is reduced via the type 2 ANOVA as shown below:

The pilot model is now fully reduced and has no need of higher order terms since this only had the purpose of screening for significant levels of main effects (2). The reason as to why the substrate is kept in the model is due to the fact that this was the main factor of interest (see more in post hoc section). The replicated model of this type can be seen in the appendix.

Gamma model with gaussian second stage model

The next model is more complicated than the previous models (lnk to the model appendix) as no explicit likelihood solution exists for it (2). In general the models of this type can be formulated as: \begin{equation} \boldsymbol{y} = X \boldsymbol{\beta} + Z \boldsymbol{u} + \boldsymbol{\epsilon} \end{equation} where X is the design matrix of the fixed effects, and the Z is the design matrix of the random effects. Where \boldsymbol{y} follows a Gamma distribution and the \boldsymbol{u} follows a gaussian distribution.
The coefficients are stored in the $\boldsymbol{\beta}$ for the fixed effect and the \boldsymbol{u} contains the random effect coefficient. The explicit likelihood formulation and Laplace approximation solution to this type of the model is seen in the likelihood appendix.

Using the same link function (log link) as presented in the gamma model (link to the model appendix), the main effect model of this type can be written as: \begin{equation} \eta_{i}=log(\mu_{i})\\ where\ \mu_{i}=\beta x_{i}+\Psi_{i} \end{equation} This translates into: \begin{equation} log(Y_{i})=\mu+\alpha(bake.time_{i})+\beta(bake.temp_{i})\\+\gamma(mix_{i})+\delta(incu.time_{i})+\zeta(sub_{i})\\+eta(W_{i})+ A(id_{i})+ \varepsilon_{i} \\ where \ \varepsilon_{i} \sim G(\alpha,\beta) \ and\\ \ id_{i} \sim N(0,\Psi^{2}) \end{equation} The response variable of the model is:

• $Y_{i}$ refers to the load measured in pascal for the i’th fungal cube

The fixed factor effects of the model are:

• $\mu$ is the overall mean
• $\alpha$ is the term describing the different baking times at the level of i=1..2
• $\beta$ is the term describing the different baking temperatures at the level of i=1..2
• $\gamma$ is the term describing the different mixing ratios at the level of i=1..3
• $\delta$ is the term describing the different incubation times at the level of i=1..2
• $\zeta$ is the term describing the different substrates at the level of i=1..2

Each term is to be understood as deviations from the overall mean

The numerical effect of the model is

• $\eta$ which is the weight of each cube $W_{i}$ corresponding to the $Y_{i}$ measurement

The random effects of the model are:

• $A(id_{i})$ which is the individual cubes tested on each tray.

From this the following reduced model is found via LRT testing, resulting in the following:

It seems the random cube to cube variation does contribute to the model as the LRT is significant and the AIC goes up if the term is removed. Thus the following model is identified: \begin{equation} Y_{i}=\mu+\beta(bake.temp_{i})+\delta(incu.time_{i})+\zeta(sub_{i})+ A(id_{i})+ \varepsilon_{i} \\ where \ \varepsilon_{i} \sim G(\alpha,\beta) \ and \ id_{i} \sim N(0,\Psi^{2}) \end{equation}

Inserting the 2 factor interactions resulted in a final model where no interaction terms were significant both for random and fixed terms. So, the final model of the gamma mixed model is the one presented above.

The residual analysis of the model is carried out the same way as the gaussian mixed model.

From the model it does not seem as if no autocorrelation is present in residuals or random effects, thus the assumption of independence is fulfilled. The normality assumption seems to hold for the random effect, but could be better for the residuals. The residual vs fitted shows that the model reasonably models out most variation in the data, but the variance is not constant as it is seen for the higher predicted values and the higher variance.

Model comparison/selection

As mentioned in the beginning, to choose the best models several parameters comes into play. The direct comparison via the AIC, BIC and goodness of fit the table is shown here:

From the table is it seen that the gamma mixed effect model is the better choice due to the goodness of fit of the data and the lower BIC and AIC. With the model choice in place, it is now time to do the post hoc analysis of the gamma mixed model.