# How to compute a non linear mix model with links between the parameters and their variance in R?

I'm working on a biological problem where we want to modelize the impact of a treatment and a chemical modification on the measured quantity of a molecule (peptides).

We will note $x_{ijk}$ the measure of this quantity, where $i\in\{1,2\}$ is the treatment, $j\in\{1,2\}$ is the chemical modification and $k\in[|1,n_k|]$ are the replicates. Let $\Pi = \{ (\alpha_j)_j,(b_{ij})_{ij},(\sigma_i)_i,\sigma\}$ be the set of parameters (both random and fixed) we want to estimate.

I want to compute the following mixed non-linear model :

(M1) : $\forall i,j,k : \left\{ \begin{array}{ll} x_{ijk} = \mu + g(\Pi) + \alpha_j + b_{ij} + \theta_{ik} + \epsilon_{ijk} \\ b_{1j} = 0 ,& \theta_{ij} \sim \mathcal{N}(0,\sigma_i)\\ \alpha_1 = 0, & \epsilon_{ijk} \sim \mathcal{N}(0,\sigma) \end{array} \right.$

And $g$ is defined thanks to a biological equation as :

$g : \Pi \rightarrow \left\{ \begin{array}{lcl} e^{0.5\times log^2(2)\times(\sigma_2^2-\sigma_1^2)}\times\frac{1+2^{b_{22}}}{1+2^{b_{21}}} & if & i = 2 \\ 0 & if & i = 1 \end{array} \right.$

To make it more understandable : $\alpha_j$ is the common effect of the chemical modification, $b_{ij}$ is the effect of the chemical modification in each treatment, $\theta_{ik}$ is the random effect (due to the replicates) in each treatment and $\epsilon_{ijk}$ is the residual error. The term $g(\Pi)$ can be seen as the effect of the treatment or a constraint between the parameters.

My question is : How can I compute the parameters and the variance in R (basically i want to estimate $\Pi$)?

I tried to use the packages nlme and saemix but they do not seems to provide solutions for non-linear mixed models when, in the model, there are links between the parameters of the model and their variance. Is there a solution in with a bayesian approach?

NB: I am not interested in a solution with an unspecified $g$ function but as i'm not really sure that this constraint is enough to have a unique solution i'll probably have to change a bit the $g$ function so solutions with tricky substitution will not help me a lot.

Others view of the problem (bonus)

I tried to formulate (M1) as a ill posed linear mixed model (ill posed because there is not a unique solution) plus a non linear constraint between my parameters (to reduce the number of solutions) :

Here we change a little bit $\Pi$ and $g$ to $\Pi'$ and $g'$ defined as :

$\Pi'= \{ (c_i)_i,(\alpha_j)_j,(b_{ij})_{ij},(\sigma_i)_i,\sigma\}$, we add a $c_i$ the effect of the treatment.

$g' : \Pi' \rightarrow c_2 - e^{0.5\times log^2(2)\times(\sigma_2^2-\sigma_1^2)}\times\frac{1+2^{b_{22}}}{1+2^{b_{21}}}$ is g but view as a constraint between the parameters.

Then we can rewrite (M1) as :

(M2) : $\forall i,j,k : \left\{ \begin{array}{ll} x_{ijk} = \mu + c_i + \alpha_j + b_{ij} + \theta_{ik} + \epsilon_{ijk} \\ b_{1j} = 0 ,& \theta_{ij} \sim \mathcal{N}(0,\sigma_i)\\ \alpha_1 = 0, & \epsilon_{ijk} \sim \mathcal{N}(0,\sigma)\\ c_1=0 & g'(\Pi') = 0 \end{array} \right.$

The first three lines of (M2) is the linear mixed model and the last line is the constraint.

This problem seems to be a problem of optimization with constraints but with random terms. Is there some package in R to solve those type of problems. I tried to look at solutions like a lasso model where my penality is $g'$ but I didn't found any R package which allows to do those type of resolution.

Thank a lot, if anyone have an idea or some papers/books I can read it would help me a lot!

I suspect you'll have to fit the model yourself using a more general purpose optimization package. The power of a non-linear mixed model is that the parameters of the non-linear function can be specified as varying across individuals and groups. This, however, is no the case in your model, since $\theta_{ik}$ is not a part of $g$. Assuming that both errors are iid, then instead you can fit the model $$x_{ijk}= g(z_{ijk},\Pi, \sigma_1, \sigma_2)+\epsilon_i$$ where $\epsilon_i\sim\mathcal{N}(0, \sigma^2+\sigma_i^2)$ and $z_{ijk}$ is a vector of indicator variables (I moved the variance parameters outside of $\Pi$ for notational convenience). As a log-likelihood, we had that $$\log P(X|Z,\Pi) = \frac{1}{2}\sum_{i,j,k}\frac{(x_{ijk}-g(z_{ijk},\Pi, \sigma_1, \sigma_2))^2}{\sigma^2+\sigma_i^2}$$. If, however, you want the $\theta_{ij}$ jointly distributed as $\theta_i\sim\mathcal{N}(0,\Sigma_i)$ and the $e_{ijk}$ as $\epsilon\sim\mathcal{N}(0,\Sigma)$, in vector notation, $$\log P(X|Z,\Pi)=\frac{1}{2}(x-g(Z,\Pi,\text{diag}(C)))^TC^{-1}(x-g(Z,\Pi,\text{diag}(C)))$$ where $C = \Sigma + \Sigma_1\oplus\Sigma_2$, we constrain all matrices to have a constant diagonal, and the $x$ vector is partitioned appropriately. This case is a little harder, but still workable, though requires some tricks from convex optimization to maintain positive definiteness in the covariance matrices.