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}