www.gusucode.com > stats 源码程序 matlab案例代码 > stats/RandomInterceptModelExample.m

    %% Fit Random Intercept LME Model  

%% 
% 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(R)
% searches, plus a nationwide estimate from the Centers 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, combine the nine columns corresponding to the
% regions into an array. The new dataset array, |flu2|, must have the
% new response variable |FluRate|, the nominal variable |Region| that shows
% which region each estimate is from, the nationwide estimate |WtdILI|,
% and the grouping variable |Date|. 
flu2 = stack(flu,2:10,'NewDataVarName','FluRate',...
    'IndVarName','Region');
flu2.Date = nominal(flu2.Date);  

%% 
% Display the first six rows of |flu2|. 
flu2(1:6,:)  

%% 
% Fit a linear mixed-effects model with a fixed-effects term for the nationwide
% estimate, |WtdILI|, and a random intercept that varies by |Date|. The
% model corresponds to  
%
% $${y_{im}} = {\beta _0} + {\beta _1}WtdIL{I_{im}} + {b_{0m}} +
% {\varepsilon _{im}},\quad i = 1,2,...,468,\quad m = 1,2,...,52,$$
%
% where $y_{im}$ is the observation $i$ for level $m$
% of grouping variable |Date|. $b_{0m}$ is the random
% effect for level $m$ of the grouping variable |Date| and $\epsilon_{im}$
% is the observation error for observation $i$. The random effect has the
% prior distribution, $b$ ~ N(0, $\sigma^{2}_{b}$ )
% and the error term has the distribution, $\epsilon$ ~ N(0, $\sigma^{2}$ ). 
lme = fitlme(flu2,'FluRate ~ 1 + WtdILI + (1|Date)') 

%%
% Estimated covariance parameters are displayed in the section titled "Random
% effects covariance parameters". The estimated value of $\sigma_{b}$
% is 0.17146 and its 95% confidence interval is [0.13227, 0.22226]. Since
% this interval does not include 0, the random-effects term is significant.
% You can formally test the significance of any random-effects term using
% a likelihood ratio test via the |<docid:stats_ug.btvrk6m compare>| method.  
%% 
% The estimated response at an observation is the sum of the fixed effects
% and the random-effect value at the grouping variable level corresponding
% to that observation. For example, the estimated flu rate for observation
% 28 is 
% 
% $$\begin{array}{l}
% {{\hat y}_{28}} = {{\hat \beta }_0} + {{\hat \beta }_1}WtdIL{I_{28}} + {{\hat b}_{10/30/2005}}\\
% \quad \quad = 0.1639 + 0.7236*(1.343) + 0.3318\\
% \quad \quad = 1.46749,
% \end{array}$$  
%
% where ${\hat b}$  is the estimated
% best linear unbiased predictor (BLUP) of the random effects for the intercept.
% 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(2)*flu2.WtdILI(28) + STATS.Estimate(STATS.Level=='10/30/2005')  

%% 
% You can display the fitted value using the |fitted| method. 
F = fitted(lme);
F(28)