Load the sample data.

```
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

`Specifying Mixed-Effects Models`

page discusses a useful model for this type of data.

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

Super impose 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 = 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 overparametrization 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 = 2.8277 0.7729 0.4606 -1.3459 PSI = 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 = 2.8277 0.7728 0.4605 -1.3460 PSI = 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 = 2.8148 0.8293 0.5613 -1.1407 PSI = 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.9462 bic: -97.2368 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) 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');

RHO = 1.0000 0.9316 0.4707 0.9316 1.0000 0.1178 0.4707 0.1178 1.0000

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 [phi,PSI,stats,b] = nlmefit(time,concentration,subject, ... [],nlme_model,phi0, ... 'REParamsSelect',[1 2 3], ... 'CovPattern',P)

P = 1 1 0 1 1 0 0 0 1 phi = 2.7830 0.8981 0.6581 -1.0000 PSI = 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 = -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.

clf; normplot(RES(:))