Load the sample data.

load flu

The `flu` dataset array has a `Date` variable,
and 10 variables containing estimated influenza rates (in 9 different
regions, estimated from Google^{®} searches, plus a nationwide estimate
from the Center for Disease Control and Prevention, CDC).

To fit a linear-mixed effects model, your data must be
in a properly formatted dataset array. To fit a linear mixed-effects
model with the influenza rates as the responses and region as the
predictor variable, combine the nine columns corresponding to the
regions into a tall array. The new dataset array, `flu2`,
must have the response variable, `FluRate`, the nominal
variable, `Region`, that shows which region each
estimate is from, and the grouping variable `Date`.

flu2 = stack(flu,2:10,'NewDataVarName','FluRate',...
'IndVarName','Region');
flu2.Date = nominal(flu2.Date);

Fit a linear mixed-effects model with fixed effects for
region and a random intercept that varies by `Date`.

Region is a categorical variable. You can specify the contrasts
for categorical variables using the `DummyVarCoding` name-value
pair argument when fitting the model. When you do not specify the
contrasts, `fitlme` uses the `'reference'` contrast
by default. Because the model has an intercept, `fitlme` takes
the first region, `NE`, as the reference and creates
eight dummy variables representing the other eight regions. For example, *I*[*MidAtl*]
is the dummy variable representing the region `MidAtl`.
For details, see Dummy Indicator Variables.

The corresponding model is

where y_{im} is
the observation *i* for level *m* of
grouping variable `Date`, *β*_{j}, *j* =
0, 1, ..., 8, are the fixed-effects coefficients, with *β*_{0} being
the coefficient for region `NE`. *b*_{0m} is
the random effect for level *m* of the grouping variable `Date`,
and *ε*_{im} is
the observation error for observation *i*. The random
effect has the prior distribution, *b*_{0m} ~
N(0,σ_{b}^{2})
and the error term has the distribution, *ε*_{im} ~
N(0,σ^{2}).

lme = fitlme(flu2,'FluRate ~ 1 + Region + (1|Date)')

Linear mixed-effects model fit by ML
Model information:
Number of observations 468
Fixed effects coefficients 9
Random effects coefficients 52
Covariance parameters 2
Formula:
FluRate ~ 1 + Region + (1|Date)
Model fit statistics:
AIC BIC LogLikelihood Deviance
318.71 364.35 -148.36 296.71
Fixed effects coefficients (95% CIs):
Name Estimate SE tStat DF pValue Lower Upper
'(Intercept)' 1.2233 0.096678 12.654 459 1.085e-31 1.0334 1.4133
'Region_MidAtl' 0.010192 0.052221 0.19518 459 0.84534 -0.092429 0.11281
'Region_ENCentral' 0.051923 0.052221 0.9943 459 0.3206 -0.050698 0.15454
'Region_WNCentral' 0.23687 0.052221 4.5359 459 7.3324e-06 0.13424 0.33949
'Region_SAtl' 0.075481 0.052221 1.4454 459 0.14902 -0.02714 0.1781
'Region_ESCentral' 0.33917 0.052221 6.495 459 2.1623e-10 0.23655 0.44179
'Region_WSCentral' 0.069 0.052221 1.3213 459 0.18705 -0.033621 0.17162
'Region_Mtn' 0.046673 0.052221 0.89377 459 0.37191 -0.055948 0.14929
'Region_Pac' -0.16013 0.052221 -3.0665 459 0.0022936 -0.26276 -0.057514
Random effects covariance parameters (95% CIs):
Group: Date (52 Levels)
Name1 Name2 Type Estimate Lower Upper
'(Intercept)' '(Intercept)' 'std' 0.6443 0.5297 0.78368
Group: Error
Name Estimate Lower Upper
'Res Std' 0.26627 0.24878 0.285

The *p*-values 7.3324e-06 and 2.1623e-10 respectively
show that the fixed effects of the flu rates in regions `WNCentral` and `ESCentral` are
significantly different relative to the flu rates in region `NE`.

The confidence limits for the standard deviation of the random-effects
term, σ^{2}_{b},
do not include 0 (0.5297, 0.78368), which indicates that the random-effects
term is significant. You can also test the significance of the random-effects
terms using the `compare` method.

The conditional fitted response from the model at a given
observation includes contributions from fixed and random effects.
For example, the estimated best linear unbiased predictor (BLUP) of
the flu rate for region `WNCentral` in week 10/9/2005
is

This is the fitted conditional response, since
it includes contributions to the estimate from both the fixed and
random effects. You can compute this value as follows.

beta = fixedEffects(lme);
[~,~,STATS] = randomEffects(lme); % Compute the random-effects statistics (STATS)
STATS.Level = nominal(STATS.Level);
y_hat = beta(1) + beta(4) + STATS.Estimate(STATS.Level=='10/9/2005')

y_hat =
1.2884

In the previous calculation, `beta(1)` corresponds
to the estimate for *β*_{0} and `beta(4)` corresponds
to the estimate for *β*_{3} You
can simply display the fitted value using the `fitted` method.

F = fitted(lme);
F(flu2.Date == '10/9/2005' & flu2.Region == 'WNCentral')

ans =
1.2884

The estimated marginal response for region `WNCentral` in
week 10/9/2005 is

Compute the fitted marginal response.

F = fitted(lme,'Conditional',false);
F(flu2.Date == '10/9/2005' & flu2.Region == 'WNCentral')

ans =
1.4602