## An application of non-linear mixed-effects models using MATLAB

During the last year of my Ph.D., I did some research collaborations with scientists from the USDA-ARS (Agricultural Research Service) at Fort Collins, CO. In addition to enjoying FoCo’s ample gastronomic options, I helped them analyze nitrogen use in corn production using data from long term experiments.

So, how did we do this? The answer is very simple: by doing the most classic AgEcon thing someone can do… estimating corn production functions! This collaboration openned up the doors for me to learn more about Non-linear models and combine them with the so-called Mixed-effects models a.k.a Random-effects models / Random coefficient models / Random parameter models / Hierarchical models / Multilevel models / Nested data models… phew, that is a lot of names but they basically mean the same thing: models that allow parameter values to vary across the population, according to some pre-specified distribution.

This type of model specification (non-linear + random-effects) makes the model estimation much more complex, and as you might suspect, it is not commonly taught in the regular econometrics courses one takes during grad school (but I might be wrong about this…). During my research and learning curve, I quickly realized most of the people were using SAS to solve non-linear mixed effects (nlme) models. This represented an additional challenge for me, as I have never been a big fan of SAS and I didn’t want to learn it – my big excuse: “I was trying to finish my dissertation and I didn’t have the time for that”.

The objective of this post is by no means to criticize SAS, neither to discuss the theory behind non-linear mixed-effects models, rather I would like to use these lines to show you how you can use MATLAB to estimate this type of models. As mentioned before, I didn’t want to use SAS, so I had to develope all the code for estimating the nlme in MATLAB and I will provide it here in this post. Why MATLAB? Simply because the previous semester I took a bayesian econometrics class and it was “very fresh” in my head. In addition, back then I didn’t have the latest Stata version that now allows you to estimate nlme (haven’t tried it yet). Of course our good friend R is also capable of estimating nlme, and I did use it a couple times, but in R the nlme models (i) converged much slower when compared to MATLAB, and (ii) required the re-scaling of the data. These are some clear disadvantages of R when compared to MATLAB when estimating nlme, especially if you have to do some bootstrapping to estimate your standard errors.

The production function we were trying to estimate was the Linear Response Stochastic Plateau (LRSP), originally developed in this paper published in 2008 in AJAE:

Tembo et al. (2008) proposed a yield response function to a single input represented by:

This function describes a Linear plus Plateau production function with year-based random effects of this form:

LRSP specifications are highly flexible and can incorporate additional random effects (see Tumussime et al., 2011). However, the use of the most flexible specifications is limited by data availability and the convergence of the nlme model.

To provide a practical guidance to estimate this nlme model in MATLAB I generated artifical data to simulate this exercise. The data I generated simulate corn yields over 10 years, using nine nitrogen levels, and each level with five repetitions. The data look like this:

When grouped by year, it looks like this (I know it is a little hard to visualize, but you can get the idea that each year have different effects on the yields):

I loaded the data and generate the previous figure using this code:

First I estimated the linear plus plateau model without random effects, that is, a plain non-linear model. Here, I provide very vague initial values for the estimation. In addition I create a figure to see the model results super-imposed on the data. I will later use these results as initial values for the more complex nlme model.

To make a quick assessment of the residuals, I draw the box-plot of residuals by year.

The box plot of residuals by year shows that the boxes are mostly above or below zero, indicating that the model has failed to account for year-specific effects. This motivates the use of year-based random effects and nlme models.

Now I estimate the linear plus plateau model and include random effects based on years. This is the non-linear mixed effects model. I use the |nlmefit| function to fit a non-linear mixed-effects model to the data. You can also use |nlmefitsa| in place of |nlmefit|. I use nlmefit because its output and reports are easier to follow.

The following anonymous function, |NLMEmodel| , adapts the three-parameter model used by |fitnlm| to the calling syntax of |nlmefit| by allowing separate parameters for each YEAR.

By default, |nlmefit| assigns random effects to all the model parameters and assumes a diagonal covariance matrix of the random effects (no covariance among the random effects). This with the objective to avoid overparametrization and related convergence issues. Here is how you use the |nlmefit| function to fit the LRSP model to all of the data, including year-specific effects. For this nlme model I used the initial values previously calculated in NL1 (betaNL1).

Next, I use REParamsSelect to specify Random effects for only the first parameter (intercept) and third parameter (yield plateau). This specification is close to that presented in Tembo et al., 2008.

I also estimated a model with a full covariance matrix. This allows for identification of correlations among the random effects. To do this, use the |CovPattern| option to specify the pattern of nonzero elements in the covariance matrix. This specification is close to that presented in Tumusiime et al., 2011.

The default method to approximate the likelihood of the model is LME. LME uses the likelihood for the linear mixed-effects model at the current conditional estimates of beta and the random effects for the m groups (years). NLMEFIT also provides the option to use FOCE “ First order Laplacian approximation” at the conditional estimates of the random effects, which is equivalent to Gaussian quadrature.

I plotted the residuals to check if the new model actually accounted for year-specific effects:

The box plot of residuals by year shows that the boxes are centered at zero, indicating that the model has accounted for year-specific effects.

To make the analysis of random effects easier, the structure in the covariance matrix is more apparent if you convert |PSI|to a correlation matrix using |corrcov|, i.e.:

RHO4 = corrcov(PSI4)

You can also plot a histogram of the individual errors for a visual assessment of the pdf:

In addition, you can create a normal probability plot of the individual errors using:

And why not calculate the Skewness and Kurtosis?