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

    %% ANOVA with Random Effects
% This example shows how to use |anovan| to fit models where a factor's
% levels represent a random selection from a larger (infinite) set of
% possible levels.
%%
% In an ordinary ANOVA model, each grouping variable represents a fixed
% factor. The levels of that factor are a fixed set of values. The goal is
% to determine whether different factor levels lead to different response
% values.
%% Set Up the Model
% Load the sample data.

% Copyright 2015 The MathWorks, Inc.

load mileage
%%
% The |anova2| function works only with balanced data, and it infers the
% values of the grouping variables from the row and column numbers of the
% input matrix. The |anovan| function, on the other hand, requires you to
% explicitly create vectors of grouping variable values. Create these
% vectors in the following way.
%%
% Create an array indicating the factory for each value in mileage. This
% array is 1 for the first column, 2 for the second, and 3 for the third.
factory  = repmat(1:3,6,1);
%%
% Create an array indicating the car model for each mileage value. This
% array is 1 for the first three rows of mileage, and 2 for the remaining
% three rows.
carmod = [ones(3,3); 2*ones(3,3)];
%%
% Turn these matrices into vectors and display them.
mileage = mileage(:);
factory = factory(:);
carmod = carmod(:);
[mileage factory carmod]
%% Fit a Random Effects Model
% Suppose you are studying a few factories but you want information about
% what would happen if you build these same car models in a different
% factory, either one that you already have or another that you might
% construct. To get this information, fit the analysis of variance model,
% specifying a model that includes an interaction term and that the factory
% factor is random.
[pvals,tbl,stats] = anovan(mileage, {factory carmod}, ... 
'model',2, 'random',1,'varnames',{'Factory' 'Car Model'});
%%
% In the fixed effects version of this fit, which you get by omitting the
% inputs |'random',1| in the preceding code, the effect of car model is
% significant, with a _p_ -value of 0.0039. But in this example, which takes
% into account the random variation of the effect of the variable |'Car
% Model'| from one factory to another, the effect is still significant, but
% with a higher _p_ -value of 0.0136.
%% F-Statistics for Models with Random Effects
% The _F_ -statistic in a model having random effects is defined
% differently than in a model having all fixed effects. In the fixed
% effects model, you compute the _F_ -statistic for any term by taking the
% ratio of the mean square for that term with the mean square for error. In
% a random effects model, however, some _F_ -statistics use a different
% mean square in the denominator.
%%
% In the example described in *Setting Up the Model* , the effect of the
% variable |'Factory'| could vary across car models. In this case, the
% interaction mean square takes the place of the error mean square in the
% _F_ -statistic.
%%
% Find the _F_ -statistic.
F = 26.6756 / 0.02
%%
% The degrees of freedom for the statistic are the degrees of freedom for
% the numerator (2) and denominator (2) mean squares.
%%
% Find the _p_ -value.
pval = 1 - fcdf(F,2,2)
%%
% With random effects, the expected value of each mean square depends not
% only on the variance of the error term, but also on the variances
% contributed by the random effects. You can see these dependencies by
% writing the expected values as linear combinations of contributions from
% the various model terms. 
%%
% Find the coefficients of these linear combinations.
stats.ems
%%
%  This returns the |ems| field of the |stats| structure.
%%
% Display text representations of the linear combinations.
stats.txtems
%%
% The expected value for the mean square due to car model (second term)
% includes contributions from a quadratic function of the car model
% effects, plus three times the variance of the interaction term's effect,
% plus the variance of the error term. Notice that if the car model effects
% were all zero, the expression would reduce to the expected mean square
% for the third term (the interaction term). That is why the _F_ -statistic
% for the car model effect uses the interaction mean square in the
% denominator.
%%
% In some cases there is no single term whose expected value matches the
% one required for the denominator of the _F_ -statistic. In that case, the
% denominator is a linear combination of mean squares. The stats structure
% contains fields giving the definitions of the denominators for each _F_
% -statistic. The |txtdenom| field, |stats.txtdenom| , contains a text
% representation, and the |denom| field contains a matrix that defines a
% linear combination of the variances of terms in the model. For balanced
% models like this one, the |denom| matrix, |stats.denom| , contains zeros
% and ones, because the denominator is just a single term's mean square.
%%
% Display the |txtdenom| field.
stats.txtdenom
%%
% Display the |denom| field.
stats.denom
%% Variance Components
% For the model described in *Setting Up the Model* , consider the mileage
% for a particular car of a particular model made at a random factory. The
% variance of that car is the sum of components, or contributions, one from
% each of the random terms.
%%
% Display the names of the random terms.
stats.rtnames
%%
% You do not know the variances, but you can estimate them from the data.
% Recall that the |ems| field of the |stats| structure expresses the expected
% value of each term's mean square as a linear combination of unknown
% variances for random terms, and unknown quadratic forms for fixed terms.
% If you take the expected mean square expressions for the random terms,
% and equate those expected values to the computed mean squares, you get a
% system of equations that you can solve for the unknown variances. These
% solutions are the variance component estimates.
%%
% Display the variance component estimate for each term.
stats.varest
%%
% Under some conditions, the variability attributed to a term is unusually
% low, and that term's variance component estimate is negative. In those
% cases it is common to set the estimate to zero, which you might do, for
% example, to create a bar graph of the components.
%%
% Create a bar graph of the components.
bar(max(0,stats.varest))
gca.xtick = 1:3;
gca.xticklabel = stats.rtnames;
%%
% You can also compute confidence bounds for the variance estimate. The
% |anovan| function does this by computing confidence bounds for the variance
% expected mean squares, and finding lower and upper limits on each
% variance component containing all of these bounds. This procedure leads
% to a set of bounds that is conservative for balanced data. (That is, 95%
% confidence bounds will have a probability of at least 95% of containing
% the true variances if the number of observations for each combination of
% grouping variables is the same.) For unbalanced data, these are
% approximations that are not guaranteed to be conservative.
%%
% Display the variance estimates and the confidence limits for the variance
% estimates of each component.
[{'Term' 'Estimate' 'Lower' 'Upper'};
 stats.rtnames, num2cell([stats.varest stats.varci])]