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.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.