MATLAB Examples

## EAXAMPLE: Split-Plot Data

`Data analysis by Linear Mixed Model with MATLAB function HPMIXED`
```% Author: Viktor Witkovsky (witkovsky@savba.sk) % Ver.: 06-Jan-2014 13:52:11 ```

## Data Description

Consider the following data from Stroup (1989a), which arise from a balanced split-plot design with the whole plots arranged in a randomized complete-block design. The variable A is the whole-plot factor, and the variable B is the subplot factor. A traditional analysis of these data involves the construction of the whole-plot error (A*Block) to test A and the pooled residual error (B*Block and A*B*Block) to test B and A*B.

Documentation reproduced from SAS/STAT® 9.22 User's Guide, Example 56.1 Split-Plot Design.

```clear load dsSplitPlotData ```

## Create the model design matrices

```formula = 'y ~ A + B + A:B + (1 | Block) + (1 | Block:A)'; model = hpmixedmodel(SplitPlotData,formula); model.Description = 'SplitPlotData: SAS Example'; ```

## Fit the linear mixed model by HPMIXED with limitted output

```opts.verbose = false; lmefit1 = hpmixed(model,opts); disp(lmefit1) ```
``` Title: 'Linear mixed-effects model fit by REML' DataInfo: 'SplitPlotData: SAS Example' ModelInfo: [1x1 struct] fixedEffects: [1x1 struct] randomEffects: [1x1 struct] varianceComponents: [1x1 struct] Details: [1x1 struct] loops: 47 tictoc: 0.3526 ```

## Fit the linear mixed model by HPMIXED with complete output

```opts.verbose = true; lmefit = hpmixed(model,opts); disp(lmefit) ```
``` MODEL INFORMATION Formula: 'y ~ 1 + A*B + (1 | Block) + (1 | Block:A)' FitMethod: 'REML' I_REML_used: 'Fisher Information Matrix / Exact' ddfMethod: 'Satterthwaite' observationsNum: 24 fixedEffectsLength: 6 randomEffectsLength: 16 randomEffectsNum: 2 randomEffectsDim: [4 12] varianceComponentsNum: 3 MMEmatRows: 22 MMEmatNonzeros: 109 MMEmatSparsity: 0.2252 loglik_REML: -59.8809 neg2loglik_REML: 119.7618 AIC: 125.7618 BIC: 123.9207 covergenceCrit: 2.0060e-11 covergenceTol: 2.3182e-11 FIXED EFFECTS Estimate SE tStat DF Quantile (Intercept) 37 4.6674 7.9273 4.6782 2.6247 A_2 1 3.5173 0.28431 8.8066 2.2698 A_3 -11 3.5173 -3.1274 8.8066 2.2698 B_2 -8.25 2.1635 -3.8133 9 2.2622 A_2:B_2 0.5 3.0596 0.16342 9 2.2622 A_3:B_2 7.75 3.0596 2.533 9 2.2622 pValue Lower Upper (Intercept) 0.00069637 24.75 49.25 A_2 0.78275 -6.9834 8.9834 A_3 0.012498 -18.983 -3.0166 B_2 0.0041319 -13.144 -3.3559 A_2:B_2 0.8738 -6.4213 7.4213 A_3:B_2 0.032076 0.82872 14.671 RANDOM EFFECTS Estimate SE tStat DF Quantile REcoef _1 10.763 4.4865 2.399 3.8212 2.8284 REcoef _2 -0.52686 4.4865 -0.11743 3.8212 2.8284 REcoef _3 -5.645 4.4865 -1.2582 3.8212 2.8284 REcoef _4 -4.5912 4.4865 -1.0233 3.8212 2.8284 REcoef _5 3.7276 3.0331 1.229 5.2134 2.5393 REcoef _6 -1.4476 3.0331 -0.47726 5.2134 2.5393 REcoef _7 0.37331 3.0331 0.12308 5.2134 2.5393 REcoef _8 -3.7171 3.0331 -1.2255 5.2134 2.5393 REcoef _9 -1.2253 3.0331 -0.40397 5.2134 2.5393 REcoef _10 4.8125 3.0331 1.5866 5.2134 2.5393 REcoef _11 0.59034 3.0331 0.19463 5.2134 2.5393 REcoef _12 0.39867 3.0331 0.13144 5.2134 2.5393 REcoef _13 -2.3806 3.0331 -0.78487 5.2134 2.5393 REcoef _14 -0.6009 3.0331 -0.19811 5.2134 2.5393 REcoef _15 2.2742 3.0331 0.7498 5.2134 2.5393 REcoef _16 -2.8052 3.0331 -0.92484 5.2134 2.5393 pValue Lower Upper REcoef _1 0.077447 -1.9268 23.453 REcoef _2 0.91243 -13.217 12.163 REcoef _3 0.27971 -18.335 7.0449 REcoef _4 0.3665 -17.281 8.0987 REcoef _5 0.27165 -3.9742 11.43 REcoef _6 0.65251 -9.1495 6.2543 REcoef _7 0.90665 -7.3286 8.0752 REcoef _8 0.27285 -11.419 3.9848 REcoef _9 0.70228 -8.9272 6.4766 REcoef _10 0.17107 -2.8894 12.514 REcoef _11 0.85304 -7.1115 8.2922 REcoef _12 0.90035 -7.3032 8.1005 REcoef _13 0.46668 -10.083 5.3213 REcoef _14 0.85045 -8.3028 7.101 REcoef _15 0.48582 -5.4277 9.9761 REcoef _16 0.39584 -10.507 4.8967 VARIANCE COMPONENTS Estimate SE tStat DF Quantile pValue sigma2_1 62.396 56.538 1.1036 Inf 1.96 0.26977 sigma2_2 15.382 11.791 1.3045 Inf 1.96 0.19206 sigma2_3 9.3611 4.4129 2.1213 Inf 1.96 0.033895 Lower Upper sigma2_1 -48.417 173.21 sigma2_2 -7.7287 38.493 sigma2_3 0.71204 18.01 FURTHER Details Matrices: [1x1 struct] Options: [1x1 struct] Model: [1x1 struct] Title: 'Linear mixed-effects model fit by REML' DataInfo: 'SplitPlotData: SAS Example' ModelInfo: [1x1 struct] fixedEffects: [1x1 struct] randomEffects: [1x1 struct] varianceComponents: [1x1 struct] Details: [1x1 struct] loops: 47 tictoc: 0.0201 ```

## EXAMPLE 1: (Statistics for FIXED and RANDOM effects and FITTED values)

```STAT_FE = getStats('fixed',lmefit); disp(STAT_FE) STAT_RE = getStats('random',lmefit); disp(STAT_RE) STAT_FIT = getStats('fitted',lmefit); disp(STAT_FIT) figure plot(model.y,STAT_FIT.TABLE.Estimate,'o') grid xlabel('y observed') ylabel('y fitted') title('SplitPlot Data Fitted by HPMIXED') ```
``` Estimate SE tStat DF Quantile (Intercept) 37 4.6674 7.9273 4.6782 2.6247 A_2 1 3.5173 0.28431 8.8066 2.2698 A_3 -11 3.5173 -3.1274 8.8066 2.2698 B_2 -8.25 2.1635 -3.8133 9 2.2622 A_2:B_2 0.5 3.0596 0.16342 9 2.2622 A_3:B_2 7.75 3.0596 2.533 9 2.2622 pValue Lower Upper (Intercept) 0.00069637 24.75 49.25 A_2 0.78275 -6.9834 8.9834 A_3 0.012498 -18.983 -3.0166 B_2 0.0041319 -13.144 -3.3559 A_2:B_2 0.8738 -6.4213 7.4213 A_3:B_2 0.032076 0.82872 14.671 verbose: 1 alpha: 0.0500 Description: 'Estimated FIXED effects' colnames: {6x1 cell} effectID: [] ddfMethod: 'Satterthwaite' TABLE: [6x8 dataset] Details: [1x1 struct] options: [1x1 struct] Estimate SE tStat DF Quantile REcoef _1 10.763 4.4865 2.399 3.8212 2.8284 REcoef _2 -0.52686 4.4865 -0.11743 3.8212 2.8284 REcoef _3 -5.645 4.4865 -1.2582 3.8212 2.8284 REcoef _4 -4.5912 4.4865 -1.0233 3.8212 2.8284 REcoef _5 3.7276 3.0331 1.229 5.2134 2.5393 REcoef _6 -1.4476 3.0331 -0.47726 5.2134 2.5393 REcoef _7 0.37331 3.0331 0.12308 5.2134 2.5393 REcoef _8 -3.7171 3.0331 -1.2255 5.2134 2.5393 REcoef _9 -1.2253 3.0331 -0.40397 5.2134 2.5393 REcoef _10 4.8125 3.0331 1.5866 5.2134 2.5393 REcoef _11 0.59034 3.0331 0.19463 5.2134 2.5393 REcoef _12 0.39867 3.0331 0.13144 5.2134 2.5393 REcoef _13 -2.3806 3.0331 -0.78487 5.2134 2.5393 REcoef _14 -0.6009 3.0331 -0.19811 5.2134 2.5393 REcoef _15 2.2742 3.0331 0.7498 5.2134 2.5393 REcoef _16 -2.8052 3.0331 -0.92484 5.2134 2.5393 pValue Lower Upper REcoef _1 0.077447 -1.9268 23.453 REcoef _2 0.91243 -13.217 12.163 REcoef _3 0.27971 -18.335 7.0449 REcoef _4 0.3665 -17.281 8.0987 REcoef _5 0.27165 -3.9742 11.43 REcoef _6 0.65251 -9.1495 6.2543 REcoef _7 0.90665 -7.3286 8.0752 REcoef _8 0.27285 -11.419 3.9848 REcoef _9 0.70228 -8.9272 6.4766 REcoef _10 0.17107 -2.8894 12.514 REcoef _11 0.85304 -7.1115 8.2922 REcoef _12 0.90035 -7.3032 8.1005 REcoef _13 0.46668 -10.083 5.3213 REcoef _14 0.85045 -8.3028 7.101 REcoef _15 0.48582 -5.4277 9.9761 REcoef _16 0.39584 -10.507 4.8967 verbose: 1 alpha: 0.0500 Description: 'Estimated RANDOM effects' colnames: 'REcoef _' effectID: [] ddfMethod: 'Satterthwaite' TABLE: [16x8 dataset] Details: [1x1 struct] options: [1x1 struct] Estimate SE tStat DF Quantile pValue yhat_1 51.491 2.2975 22.412 11.098 2.1986 1.3602e-10 yhat_2 47.315 2.2975 20.595 11.098 2.1986 3.4067e-10 yhat_3 37.136 2.2975 16.164 11.098 2.1986 4.6208e-09 yhat_4 32.756 2.2975 14.258 11.098 2.1986 1.7534e-08 yhat_5 36.248 2.2975 15.777 11.098 2.1986 5.9833e-09 yhat_6 30.286 2.2975 13.182 11.098 2.1986 4.0024e-08 yhat_7 31.945 2.2975 13.905 11.098 2.1986 2.2843e-08 yhat_8 32.754 2.2975 14.257 11.098 2.1986 1.7548e-08 yhat_9 17.974 2.2975 7.8236 11.098 2.1986 7.6466e-06 yhat_10 31.808 2.2975 13.845 11.098 2.1986 2.3905e-08 yhat_11 35.683 2.2975 15.532 11.098 2.1986 7.0731e-09 yhat_12 18.604 2.2975 8.0975 11.098 2.1986 5.502e-06 yhat_13 43.241 2.2975 18.821 11.098 2.1986 9.0195e-10 yhat_14 39.565 2.2975 17.221 11.098 2.1986 2.3452e-09 yhat_15 36.636 2.2975 15.947 11.098 2.1986 5.3402e-09 yhat_16 24.506 2.2975 10.667 11.098 2.1986 3.5765e-07 yhat_17 28.498 2.2975 12.404 11.098 2.1986 7.5567e-08 yhat_18 29.786 2.2975 12.965 11.098 2.1986 4.7647e-08 yhat_19 23.695 2.2975 10.314 11.098 2.1986 5.0334e-07 yhat_20 25.004 2.2975 10.883 11.098 2.1986 2.9134e-07 yhat_21 17.474 2.2975 7.606 11.098 2.1986 9.9922e-06 yhat_22 23.558 2.2975 10.254 11.098 2.1986 5.3389e-07 yhat_23 27.933 2.2975 12.158 11.098 2.1986 9.3049e-08 yhat_24 18.104 2.2975 7.8799 11.098 2.1986 7.1421e-06 Lower Upper yhat_1 46.439 56.542 yhat_2 42.264 52.367 yhat_3 32.085 42.188 yhat_4 27.705 37.807 yhat_5 31.197 41.299 yhat_6 25.234 35.337 yhat_7 26.894 36.997 yhat_8 27.702 37.805 yhat_9 12.923 23.026 yhat_10 26.757 36.859 yhat_11 30.632 40.734 yhat_12 13.552 23.655 yhat_13 38.189 48.292 yhat_14 34.514 44.617 yhat_15 31.585 41.688 yhat_16 19.455 29.557 yhat_17 23.447 33.549 yhat_18 24.734 34.837 yhat_19 18.644 28.747 yhat_20 19.952 30.055 yhat_21 12.423 22.526 yhat_22 18.507 28.609 yhat_23 22.882 32.984 yhat_24 13.052 23.155 verbose: 1 alpha: 0.0500 Description: 'FITTED yhat' colnames: 'yhat_' effectID: [] ddfMethod: 'Satterthwaite' TABLE: [24x8 dataset] Details: [1x1 struct] options: [1x1 struct] ```

## EXAMPLE 2: (t-statistics for the LSMEANS CONTRASTS for interaction A:B)

```options.STAT.inference = 'contrasts'; options.STAT.inferenceSpace = 'broad'; options.STAT.ddfMethod = 'Satterthwaite'; [Lambda,options] = getLambda({'A' 'B' },model,SplitPlotData,options); STAT_AB_CONTRASTS = getStats(Lambda,lmefit,options); disp(STAT_AB_CONTRASTS) ```
``` Estimate SE tStat DF Quantile A:B_1 1 - A:B_1 2 8.25 2.1635 3.8133 9 2.2622 A:B_1 1 - A:B_2 1 -1 3.5173 -0.28431 8.8066 2.2698 A:B_1 1 - A:B_2 2 6.75 3.5173 1.9191 8.8066 2.2698 A:B_1 1 - A:B_3 1 11 3.5173 3.1274 8.8066 2.2698 A:B_1 1 - A:B_3 2 11.5 3.5173 3.2695 8.8066 2.2698 pValue Lower Upper A:B_1 1 - A:B_1 2 0.0041319 3.3559 13.144 A:B_1 1 - A:B_2 1 0.78275 -8.9834 6.9834 A:B_1 1 - A:B_2 2 0.087899 -1.2334 14.733 A:B_1 1 - A:B_3 1 0.012498 3.0166 18.983 A:B_1 1 - A:B_3 2 0.0099755 3.5166 19.483 inference: 'contrasts' inferenceSpace: 'broad' ddfMethod: 'Satterthwaite' includedZcols: [] includedXcols: [] VarDescription: 'A:B' Description: [1x69 char] colnames: {5x1 cell} grpname: 'A:B' verbose: 1 alpha: 0.0500 effectID: [] TABLE: [5x8 dataset] Details: [1x1 struct] options: [1x1 struct] ```

## EXAMPLE 3: (FTEST for the LSMEANS CONTRASTS for interaction A:B)

```options.STAT.inference = 'contrasts'; options.STAT.inferenceSpace = 'broad'; options.STAT.ddfMethod = 'FaiCornelius'; [Lambda,options] = getLambda({'A' 'B' },model,SplitPlotData,options); STAT_AB_FTEST = getStats(Lambda,lmefit,options); disp(STAT_AB_FTEST) ```
``` Fstat NDF DDF Quantile pValue A:B 7.1133 5 8.3553 3.6072 0.0072091 inference: 'contrasts' inferenceSpace: 'broad' ddfMethod: 'FaiCornelius' includedZcols: [] includedXcols: [] VarDescription: 'A:B' Description: [1x69 char] colnames: {5x1 cell} grpname: 'A:B' verbose: 1 alpha: 0.0500 effectID: [] TABLE: [1x5 dataset] Details: [1x1 struct] options: [1x1 struct] ```

## EXAMPLE 4: (LSMEANS PAIRWISE COMPARISONS for the effect A)

```options.STAT.inference = 'pairwise'; options.STAT.inferenceSpace = 'narrow'; options.STAT.ddfMethod = 'Satterthwaite'; [Lambda,options] = getLambda({'A'},model,SplitPlotData,options); STAT_A_COMPARISONS = getStats(Lambda,lmefit,options); disp(STAT_A_COMPARISONS) ```
``` Estimate SE tStat DF Quantile pValue A_1 - A_2 -1.25 1.5298 -0.8171 9 2.2622 0.43497 A_1 - A_3 7.125 1.5298 4.6575 9 2.2622 0.0011895 A_2 - A_3 8.375 1.5298 5.4746 9 2.2622 0.00039291 Lower Upper A_1 - A_2 -4.7106 2.2106 A_1 - A_3 3.6644 10.586 A_2 - A_3 4.9144 11.836 inference: 'pairwise' inferenceSpace: 'narrow' ddfMethod: 'Satterthwaite' includedZcols: [] includedXcols: [] VarDescription: 'A' Description: [1x71 char] colnames: {3x1 cell} grpname: 'A' verbose: 1 alpha: 0.0500 effectID: [] TABLE: [3x8 dataset] Details: [1x1 struct] options: [1x1 struct] ```

## EXAMPLE 5: (Statistics for the LSMEANS for the effect A)

```options.STAT.inference = 'lsmeans'; options.STAT.inferenceSpace = 'broad'; [Lambda,options] = getLambda({'A'},model,SplitPlotData,options); STAT_A_LSMEANS = getStats(Lambda,lmefit,options); disp(STAT_A_LSMEANS) ```
``` Estimate SE tStat DF Quantile pValue A_1 32.875 4.5403 7.2407 4.1955 2.7262 0.0016042 A_2 34.125 4.5403 7.516 4.1955 2.7262 0.0013847 A_3 25.75 4.5403 5.6714 4.1955 2.7262 0.0041243 Lower Upper A_1 20.497 45.253 A_2 21.747 46.503 A_3 13.372 38.128 inference: 'lsmeans' inferenceSpace: 'broad' ddfMethod: 'Satterthwaite' includedZcols: [] includedXcols: [] VarDescription: 'A' Description: [1x57 char] colnames: {3x1 cell} grpname: 'A' verbose: 1 alpha: 0.0500 effectID: [] TABLE: [3x8 dataset] Details: [1x1 struct] options: [1x1 struct] ```

## EXAMPLE_6: (BROAD INFERENCE SPACE for LSMEANS of the factor Block)

```options.STAT.inference = 'lsmeans'; options.STAT.inferenceSpace = 'broad'; [Lambda,options] = getLambda({'Block'},model,SplitPlotData,options); options.STAT.alpha = 0.01; STAT_Block_LSMEANS_BROAD = getStats(Lambda,lmefit,options); disp(STAT_Block_LSMEANS_BROAD) ```
``` Estimate SE tStat DF Quantile pValue Block_1 30.917 4.1558 7.4394 3 5.8409 0.005027 Block_2 30.917 4.1558 7.4394 3 5.8409 0.005027 Block_3 30.917 4.1558 7.4394 3 5.8409 0.005027 Block_4 30.917 4.1558 7.4394 3 5.8409 0.005027 Lower Upper Block_1 6.6429 55.19 Block_2 6.6429 55.19 Block_3 6.6429 55.19 Block_4 6.6429 55.19 inference: 'lsmeans' inferenceSpace: 'broad' ddfMethod: 'Satterthwaite' includedZcols: [] includedXcols: [] VarDescription: 'Block' Description: [1x61 char] colnames: {4x1 cell} grpname: 'Block' alpha: 0.0100 verbose: 1 effectID: [] TABLE: [4x8 dataset] Details: [1x1 struct] options: [1x1 struct] ```

## EXAMPLE 7: (NARROW INFERENCE SPACE for LSMEANS of the factor Block)

```options.STAT.inference = 'lsmeans'; options.STAT.inferenceSpace = 'narrow'; options.STAT.alpha = 0.01; [Lambda,options] = getLambda({'Block'},model,SplitPlotData,options); STAT_Block_LSMEANS_NARROW = getStats(Lambda,lmefit,options); disp(STAT_Block_LSMEANS_NARROW) ```
``` Estimate SE tStat DF Quantile pValue Block_1 42.564 1.2385 34.369 9.3098 3.2227 3.984e-11 Block_2 30.347 1.2385 24.504 9.3098 3.2227 9.0026e-10 Block_3 24.808 1.2385 20.031 9.3098 3.2227 5.6917e-09 Block_4 25.948 1.2385 20.952 9.3098 3.2227 3.7769e-09 Lower Upper Block_1 38.573 46.555 Block_2 26.355 34.338 Block_3 20.817 28.799 Block_4 21.957 29.939 inference: 'lsmeans' inferenceSpace: 'narrow' ddfMethod: 'Satterthwaite' includedZcols: [] includedXcols: [] VarDescription: 'Block' Description: [1x62 char] colnames: {4x1 cell} grpname: 'Block' alpha: 0.0100 verbose: 1 effectID: [] TABLE: [4x8 dataset] Details: [1x1 struct] options: [1x1 struct] ```

## EXAMPLE 8: (ITERMEDIATE INFERENCE SPACE for LSMEANS of the factor Block)

```options.STAT.inference = 'lsmeans'; options.STAT.inferenceSpace = 'intermediate'; options.STAT.includedXcols = (1:size(model.X,2)); options.STAT.includedZcols = model.dimRE(1) + (1:model.dimRE(2)); options.STAT.alpha = 0.01; [Lambda,options] = getLambda({'Block'},model,SplitPlotData,options); STAT_Block_LSMEANS_ITERMEDIATE = getStats(Lambda,lmefit,options); disp(STAT_Block_LSMEANS_ITERMEDIATE) ```
``` Estimate SE tStat DF Quantile pValue Block_1 31.801 4.4214 7.1925 3.678 4.8914 0.0027118 Block_2 30.873 4.4214 6.9827 3.678 4.8914 0.0030053 Block_3 30.453 4.4214 6.8875 3.678 4.8914 0.0031515 Block_4 30.539 4.4214 6.9071 3.678 4.8914 0.0031206 Lower Upper Block_1 10.174 53.428 Block_2 9.2466 52.5 Block_3 8.826 52.08 Block_4 8.9126 52.166 inference: 'lsmeans' inferenceSpace: 'intermediate' ddfMethod: 'Satterthwaite' includedZcols: [5 6 7 8 9 10 11 12 13 14 15 16] includedXcols: [1 2 3 4 5 6] VarDescription: 'Block' Description: [1x91 char] colnames: {4x1 cell} grpname: 'Block' alpha: 0.0100 verbose: 1 effectID: [] TABLE: [4x8 dataset] Details: [1x1 struct] options: [1x1 struct] ```

## Comparison with MATLAB function FITLME

```tic; lme = fitlme(SplitPlotData,formula,'FitMethod','REML'); toc disp(lme) ```
```Elapsed time is 1.764975 seconds. Linear mixed-effects model fit by REML Model information: Number of observations 24 Fixed effects coefficients 6 Random effects coefficients 16 Covariance parameters 3 Formula: y ~ 1 + A*B + (1 | Block) + (1 | Block:A) Model fit statistics: AIC BIC LogLikelihood Deviance 137.76 145.78 -59.881 119.76 Fixed effects coefficients (95% CIs): Name Estimate SE tStat DF pValue '(Intercept)' 37 4.6674 7.9273 18 2.7903e-07 'A_2' 1 3.5173 0.28431 18 0.77942 'A_3' -11 3.5173 -3.1274 18 0.0058198 'B_2' -8.25 2.1635 -3.8133 18 0.001273 'A_2:B_2' 0.5 3.0596 0.16342 18 0.87201 'A_3:B_2' 7.75 3.0596 2.533 18 0.020826 Lower Upper 27.194 46.806 -6.3896 8.3896 -18.39 -3.6104 -12.795 -3.7047 -5.928 6.928 1.322 14.178 Random effects covariance parameters (95% CIs): Group: Block (4 Levels) Name1 Name2 Type Estimate Lower '(Intercept)' '(Intercept)' 'std' 7.8991 3.2504 Upper 19.197 Group: Block:A (12 Levels) Name1 Name2 Type Estimate Lower '(Intercept)' '(Intercept)' 'std' 3.922 1.8503 Upper 8.313 Group: Error Name Estimate Lower Upper 'Res Std' 3.0596 1.9277 4.8562 ```