Why is a mixed model a non-linear statistical model? I was wondering if anyone could elaborate more on this statement which I came across whilst reading a book on non-linear mixed effects model. I know that you can have linear, generalised linear and non-linear  mixed effect models. 
The book says 

mixed models are non-linear statistical models, due mainly to the presence of variance parameters.

How do I explain this to someone who understands very little statistics?
For instance, how can a linear mixed effect model be non-linear?
 A: And, here's a "street" version of the above:
a) what is a linear model? It's one which can be expressed in the form of sums and scalar products of the inputs (y = ax + b at the simplest).  If "a" is a function of some other factor ("z" not "x") perhaps a random function, A(z), then we have y = A(z)x + b which is no longer a linear model.
b) what is a "mixed-model?" One which has linear parts and non-linear parts:  y = ax + A(Z)x + b
A: Hopefully the amount of notation suppressed and corners cut in what follows still leaves something intelligible:
On what 'mixed' means. Imagine somewhere at the heart of the model we have a line looking something like this.
$$\eta = \beta_0 +  \beta_1 x_1 + ... +  \beta_p x_p$$
Where the $x_k$ are our covariates. Focus on the coefficients, the $\beta_k$. We could either: think of these as fixed numbers (fixed effects) or as random numbers (random effects). 
Why might we want to think of a $\beta_k$ as random? For example, imagine a situation where each person gets their own intercept term $\beta_0$. This might represent them being naturally healthier, faster, smarter or whatever. We could then model what effect a drug / treatment had on top of their individual natural level.  We might want to think of the effect of the drug or treatment though as being fixed and bumping $\eta$ up or down by the same amount for every individual. 
Hence the term mixed effects where some are random (e.g., person's natural level) and some are fixed (e.g., effect of treatment).
On linear and non-linear.
The 'linear' in generalised linear models refers to the equation for $\eta$ above.  
If $\mu$ is the mean of the variable we are modelling then we can do non-linear things with $\eta$ like:
$$\mu = \frac{1}{1+e^{-\eta}}$$
Which is involved in logistic regression. Or 
$$\mu = e^{\eta}$$
As in models for count data.
This is where having random effects in $\eta$ can cause some real bumps. Suppose we have a prior distribution for $\beta_0\sim\mathcal N(0,g)$. Then for example with logistic data:
$$E[Y]=E[E[Y|\beta_0]]=\int \frac{1}{1+e^{-\eta}} (2\pi)^{-1/2}e^{-\frac{1}{2}\beta_0^2} d \beta_0$$
Which is pretty hairy, despite the fact that $\beta_0$ was added as a linear term to $\eta$. Think about the distribution of $\beta_0$ which is symmetric, but taking $\exp()$ of it is not.
However if we model $\mu=\eta$ things work out nice and the random terms disappear from the mean upon marginalisation. They don't, however, disappear from the variance.
Suppose for some observation $i$, with a random intercept (which is independent of the error term $\varepsilon_i \sim\mathcal N(0,\sigma^2)$):
$$Y_i|\beta_0 = \beta_0 + \beta_1 x_1 + \varepsilon_i$$
Then 
$$E[Y_i]=E[E[Y_i|\beta_0]]=E[\beta_0 + \beta_1 x_1]=\beta_1 x_1$$
But 
$${\rm Var}(Y_i)={\rm Var}(\beta_0 + \beta_1 x_1 + \varepsilon_i) = {\rm Var}(\beta_0) + {\rm Var}( \varepsilon_i) = g + \sigma^2$$
So the 'variance parameters' still show up.
