# Mixed-Effects Models Using nlmefit and nlmefitsa

This example fits mixed-effects models, plots predictions and residuals, and interprets the results.

load indomethacin

The data in indomethacin.mat records concentrations of the drug indomethacin in the bloodstream of six subjects over eight hours.

Plot the scatter plot of indomethacin in the bloodstream grouped by subject.

gscatter(time,concentration,subject) xlabel('Time (hours)') ylabel('Concentration (mcg/ml)') title('{\bf Indomethacin Elimination}') hold on

Including random effects in a model is effective when data falls into natural groups. In this data, the groups are simply the individuals under study. For details on mixed-effects models, which account for fixed effects and random effects, see Mixed-Effects Models.

Construct the model via an anonymous function.

model = @(phi,t)(phi(1)*exp(-exp(phi(2))*t) + ... phi(3)*exp(-exp(phi(4))*t));

Use the nlinfit function to fit the model to all of the data, ignoring subject-specific effects.

phi0 = [1 2 1 1]; [phi,res] = nlinfit(time,concentration,model,phi0);

Compute the mean squared error.

numObs = length(time); numParams = 4; df = numObs-numParams; mse = (res'*res)/df
mse = 0.0304 

Superimpose the model on the scatter plot of data.

tplot = 0:0.01:8; plot(tplot,model(phi,tplot),'k','LineWidth',2) hold off

Draw the box plot of residuals by subject.

colors = 'rygcbm'; h = boxplot(res,subject,'colors',colors,'symbol','o'); set(h(~isnan(h)),'LineWidth',2) hold on boxplot(res,subject,'colors','k','symbol','ko') grid on xlabel('Subject') ylabel('Residual') hold off

The box plot of residuals by subject shows that the boxes are mostly above or below zero, indicating that the model has failed to account for subject-specific effects.

To account for subject-specific effects, fit the model separately to the data for each subject.

phi0 = [1 2 1 1]; PHI = zeros(4,6); RES = zeros(11,6); for I = 1:6 tI = time(subject == I); cI = concentration(subject == I); [PHI(:,I),RES(:,I)] = nlinfit(tI,cI,model,phi0); end PHI
PHI = 4×6 2.0293 2.8277 5.4683 2.1981 3.5661 3.0023 0.5794 0.8013 1.7498 0.2423 1.0408 1.0882 0.1915 0.4989 1.6757 0.2545 0.2915 0.9685 -1.7878 -1.6354 -0.4122 -1.6026 -1.5069 -0.8731 

Compute the mean squared error.

numParams = 24; df = numObs-numParams; mse = (RES(:)'*RES(:))/df
mse = 0.0057 

Plot the scatter plot of the data and superimpose the model for each subject.

gscatter(time,concentration,subject) xlabel('Time (hours)') ylabel('Concentration (mcg/ml)') title('{\bf Indomethacin Elimination}') hold on for I = 1:6 plot(tplot,model(PHI(:,I),tplot),'Color',colors(I)) end axis([0 8 0 3.5]) hold off

PHI gives estimates of the four model parameters for each of the six subjects. The estimates vary considerably, but taken as a 24-parameter model of the data, the mean-squared error of 0.0057 is a significant reduction from 0.0304 in the original four-parameter model.

Draw the box plot of residuals by subject.

h = boxplot(RES,'colors',colors,'symbol','o'); set(h(~isnan(h)),'LineWidth',2) hold on boxplot(RES,'colors','k','symbol','ko') grid on xlabel('Subject') ylabel('Residual') hold off

Now the box plot shows that the larger model accounts for most of the subject-specific effects. The spread of the residuals (the vertical scale of the box plot) is much smaller than in the previous box plot, and the boxes are now mostly centered on zero.

While the 24-parameter model successfully accounts for variations due to the specific subjects in the study, it does not consider the subjects as representatives of a larger population. The sampling distribution from which the subjects are drawn is likely more interesting than the sample itself. The purpose of mixed-effects models is to account for subject-specific variations more broadly, as random effects varying around population means.

Use the nlmefit function to fit a mixed-effects model to the data. You can also use nlmefitsa in place of nlmefit .

The following anonymous function, nlme_model , adapts the four-parameter model used by nlinfit to the calling syntax of nlmefit by allowing separate parameters for each individual. By default, nlmefit assigns random effects to all the model parameters. Also by default, nlmefit assumes a diagonal covariance matrix (no covariance among the random effects) to avoid overparameterization and related convergence issues.

nlme_model = @(PHI,t)(PHI(:,1).*exp(-exp(PHI(:,2)).*t) + ... PHI(:,3).*exp(-exp(PHI(:,4)).*t)); phi0 = [1 2 1 1]; [phi,PSI,stats] = nlmefit(time,concentration,subject, ... [],nlme_model,phi0)
phi = 4×1 2.8277 0.7729 0.4606 -1.3459 
PSI = 4×4 0.3264 0 0 0 0 0.0250 0 0 0 0 0.0124 0 0 0 0 0.0000 
stats = struct with fields: dfe: 57 logl: 54.5882 mse: 0.0066 rmse: 0.0787 errorparam: 0.0815 aic: -91.1765 bic: -93.0506 covb: [4x4 double] sebeta: [0.2558 0.1066 0.1092 0.2244] ires: [66x1 double] pres: [66x1 double] iwres: [66x1 double] pwres: [66x1 double] cwres: [66x1 double] 

The mean-squared error of 0.0066 is comparable to the 0.0057 of the 24-parameter model without random effects, and significantly better than the 0.0304 of the four-parameter model without random effects.

The estimated covariance matrix PSI shows that the variance of the fourth random effect is essentially zero, suggesting that you can remove it to simplify the model. To do this, use the 'REParamsSelect' name-value pair to specify the indices of the parameters to be modeled with random effects in nlmefit .

[phi,PSI,stats] = nlmefit(time,concentration,subject, ... [],nlme_model,phi0, ... 'REParamsSelect',[1 2 3])
phi = 4×1 2.8277 0.7728 0.4605 -1.3460 
PSI = 3×3 0.3270 0 0 0 0.0250 0 0 0 0.0124 
stats = struct with fields: dfe: 58 logl: 54.5875 mse: 0.0066 rmse: 0.0780 errorparam: 0.0815 aic: -93.1750 bic: -94.8410 covb: [4x4 double] sebeta: [0.2560 0.1066 0.1092 0.2244] ires: [66x1 double] pres: [66x1 double] iwres: [66x1 double] pwres: [66x1 double] cwres: [66x1 double] 

The log-likelihood logl is almost identical to what it was with random effects for all of the parameters, the Akaike information criterion aic is reduced from -91.1765 to -93.1750, and the Bayesian information criterion bic is reduced from -93.0506 to -94.8410. These measures support the decision to drop the fourth random effect.

Refitting the simplified model with a full covariance matrix allows for identification of correlations among the random effects. To do this, use the CovPattern parameter to specify the pattern of nonzero elements in the covariance matrix.

[phi,PSI,stats] = nlmefit(time,concentration,subject, ... [],nlme_model,phi0, ... 'REParamsSelect',[1 2 3], ... 'CovPattern',ones(3))
phi = 4×1 2.8148 0.8293 0.5613 -1.1407 
PSI = 3×3 0.4767 0.1152 0.0499 0.1152 0.0321 0.0032 0.0499 0.0032 0.0236 
stats = struct with fields: dfe: 55 logl: 58.4731 mse: 0.0061 rmse: 0.0782 errorparam: 0.0781 aic: -94.9463 bic: -97.2369 covb: [4x4 double] sebeta: [0.3028 0.1103 0.1179 0.1662] ires: [66x1 double] pres: [66x1 double] iwres: [66x1 double] pwres: [66x1 double] cwres: [66x1 double] 

The estimated covariance matrix PSI shows that the random effects on the first two parameters have a relatively strong correlation, and both have a relatively weak correlation with the last random effect. This structure in the covariance matrix is more apparent if you convert PSI to a correlation matrix using corrcov .

RHO = corrcov(PSI)
RHO = 3×3 1.0000 0.9316 0.4707 0.9316 1.0000 0.1178 0.4707 0.1178 1.0000 
clf; imagesc(RHO) set(gca,'XTick',[1 2 3],'YTick',[1 2 3]) title('{\bf Random Effect Correlation}') h = colorbar; set(get(h,'YLabel'),'String','Correlation');

Incorporate this structure into the model by changing the specification of the covariance pattern to block-diagonal.

P = [1 1 0;1 1 0;0 0 1] % Covariance pattern
P = 3×3 1 1 0 1 1 0 0 0 1 
[phi,PSI,stats,b] = nlmefit(time,concentration,subject, ... [],nlme_model,phi0, ... 'REParamsSelect',[1 2 3], ... 'CovPattern',P)
phi = 4×1 2.7830 0.8981 0.6581 -1.0000 
PSI = 3×3 0.5180 0.1069 0 0.1069 0.0221 0 0 0 0.0454 
stats = struct with fields: dfe: 57 logl: 58.0804 mse: 0.0061 rmse: 0.0768 errorparam: 0.0782 aic: -98.1608 bic: -100.0350 covb: [4x4 double] sebeta: [0.3171 0.1073 0.1384 0.1453] ires: [66x1 double] pres: [66x1 double] iwres: [66x1 double] pwres: [66x1 double] cwres: [66x1 double] 
b = 3×6 -0.8507 -0.1563 1.0427 -0.7559 0.5652 0.1550 -0.1756 -0.0323 0.2152 -0.1560 0.1167 0.0320 -0.2756 0.0519 0.2620 0.1064 -0.2835 0.1389 

The block-diagonal covariance structure reduces aic from -94.9462 to -98.1608 and bic from -97.2368 to -100.0350 without significantly affecting the log-likelihood. These measures support the covariance structure used in the final model. The output b gives predictions of the three random effects for each of the six subjects. These are combined with the estimates of the fixed effects in phi to produce the mixed-effects model.

Plot the mixed-effects model for each of the six subjects. For comparison, the model without random effects is also shown.

PHI = repmat(phi,1,6) + ... % Fixed effects [b(1,:);b(2,:);b(3,:);zeros(1,6)]; % Random effects RES = zeros(11,6); % Residuals colors = 'rygcbm'; for I = 1:6 fitted_model = @(t)(PHI(1,I)*exp(-exp(PHI(2,I))*t) + ... PHI(3,I)*exp(-exp(PHI(4,I))*t)); tI = time(subject == I); cI = concentration(subject == I); RES(:,I) = cI - fitted_model(tI); subplot(2,3,I) scatter(tI,cI,20,colors(I),'filled') hold on plot(tplot,fitted_model(tplot),'Color',colors(I)) plot(tplot,model(phi,tplot),'k') axis([0 8 0 3.5]) xlabel('Time (hours)') ylabel('Concentration (mcg/ml)') legend(num2str(I),'Subject','Fixed') end

If obvious outliers in the data (visible in previous box plots) are ignored, a normal probability plot of the residuals shows reasonable agreement with model assumptions on the errors.

figure normplot(RES(:))