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

    %% Two Grouping Variables  

% Copyright 2015 The MathWorks, Inc.


%% 
% Load the sample data. 
load carbig  

%% 
% Fit a linear mixed-effects model for miles per gallon (MPG), with fixed
% effects for acceleration and weight, a potentially correlated random effect
% for intercept and acceleration grouped by model year, and an independent
% random effect for weight, grouped by the origin of the car. This model
% corresponds to

%%
% $\begin{array}{l}
% MP{G_{imk}} = {\beta _0} + {\beta _1}Ac{c_i} + {\beta _2}Weigh{t_i} + {b_{10m}} + {b_{11}}_mAc{c_i} + {b_{21k}}Weigh{t_i} + {\varepsilon _{imk}},\quad \\
% \quad \quad m = 1,2,...,13,\quad k = 1,2,...,8,
% \end{array}$ 


%%
% where _m_ represents the levels for the variable |Model_Year|, and _k_
% represents the levels for the variable |Origin|. $MPG_{imk}$
% is the miles per gallon for the _i_ th observation, _m_ th model year, and
% _k_ th origin that correspond to the _i_ th observation. The random-effects
% terms and the observation error have the following prior distributions:

%%
% ${b_{1m}} = \left( \begin{array}{l}
% {b_{10m}}\\
% {b_{11m}}
% \end{array} \right) \sim N\left( {0,\left( {\begin{array}{*{20}{c}}
% {\sigma _{10}^2}&{\sigma _{10,11}^{}}\\
% {\sigma _{10,11}^{}}&{\sigma _{11}^2}
% \end{array}} \right)} \right)$

%%
% $b_{2k}\,\sim N\left( {0,\sigma _2^2} \right)$

%%
% ${\varepsilon _{imk}} \sim N\left( {0,{\sigma ^2}} \right)$

%%
% Here, the random-effects term $b_{1m}$ represents the
% first random effect at level _m_ of the first grouping variable. The random-effects
% term $b_{10m}$ corresponds to the first random effects
% term (1), for the intercept (0), at the _m_ th level (_m_) of the first
% grouping variable. Likewise $b_{11m}$ is the level
% _m_ for the first predictor (1) in the first random-effects term (1). 
%
% Similarly, $b_{2k}$ stands for the second random effects-term
% at level _k_ of the second grouping variable. 
%
% $\theta^{2}_{10}$ is the variance of the
% random-effects term for the intercept, $\theta^{2}_{11}$
% is the variance of the random effects term for the predictor acceleration,
% and $\theta_{10,11}$ is the covariance of the random-effects
% terms for the intercept and the predictor acceleration. $\theta^{2}_{2}$
% is the variance of the second random-effects term, and $\theta^{2}$
% is the residual variance. 
%
% First, prepare the design matrices for fitting the linear mixed-effects
% model. 
X = [ones(406,1) Acceleration Weight];
Z = {[ones(406,1) Acceleration],[Weight]};
Model_Year = nominal(Model_Year);
Origin = nominal(Origin);
G = {Model_Year,Origin};  

%% 
% Fit the model using the design matrices. 
lme = fitlmematrix(X,MPG,Z,G,'FixedEffectPredictors',....
{'Intercept','Acceleration','Weight'},'RandomEffectPredictors',...
{{'Intercept','Acceleration'},{'Weight'}},'RandomEffectGroups',{'Model_Year','Origin'});  

%% 
% Compute the estimates of covariance parameters for the random effects. 
[psi,mse,stats] = covarianceParameters(lme) 

%%
% The residual variance |mse| is 9.0755. |psi| is a 2-by-1 cell array, and
% |stats| is a 3-by-1 cell array. To see the contents, you must index into
% these cell arrays.  

%% 
% First, index into the first cell of |psi|. 
psi{1} 

%%
% The first cell of |psi| contains the covariance parameters for the correlated
% random effects for intercept $\theta^{2}_{10}$
% as 8.5160, and for acceleration $\theta^{2}_{11}$
% as 0.1087. The estimate for the covariance of the random-effects terms
% for the intercept and acceleration $\theta_{10,11}$ is -0.8387.  

%% 
% Now, index into the second cell of |psi|. 
psi{2} 

%%
% The second cell of |psi| contains the estimate for the variance of the
% random-effects term for weight $\theta^{2}_{2}$.

%% 
% Index into the first cell of |stats|. 
stats{1} 

%%
% This table shows the standard deviation estimates for the random-effects
% terms for intercept and acceleration. Note that the standard deviations
% estimates are the square roots of the diagonal elements in the first cell
% of |psi|. Specifically, 2.9182 = sqrt(8.5160) and 0.32968 = sqrt(0.1087).
% The correlation is a function of the covariance of intercept and acceleration,
% and the standard deviations of intercept and acceleration. The covariance
% of intercept and acceleration is the off-diagonal value in the first cell
% of psi, -0.8387. So, the correlation is -0.8387/(0.32968*2.92182) = -0.87. 

%%
% The grouping variable for intercept and acceleration is |Model_Year|.  

%% 
% Index into the second cell of |stats|. 
stats{2} 

%%
% The second cell of |stats| has the standard deviation estimate and the
% 95% confidence limits for the standard deviation of the random-effects
% term for |Weight|. The grouping variable is |Origin|.  

%% 
% Index into the third cell of |stats|. 
stats{3} 

%%
% The third cell of |stats| contains the estimate for residual standard
% deviation and the 95% confidence limits. The estimate for residual standard
% deviation is the square root of |mse|, sqrt(9.0755) = 3.0126.  

%% 
% Construct 99% confidence intervals for the covariance parameters. 
[~,~,stats] = covarianceParameters(lme,'Alpha',0.01);
stats{1}  

%%  
stats{2}  

%%  
stats{3}