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

    %% LME Model for Split-Plot Experiment  

%% 
% Load the sample data. 
load(fullfile(matlabroot,'examples','stats','fertilizer.mat'))

%%
% The dataset array includes data from a split-plot experiment, where soil
% is divided into three blocks based on the soil type: sandy, silty, and
% loamy. Each block is divided into five plots, where five types of tomato
% plants (cherry, heirloom, grape, vine, and plum) are randomly assigned
% to these plots. The tomato plants in the plots are then divided into subplots,
% where each subplot is treated by one of four fertilizers. This is simulated
% data.  

%% 
% Store the data in a dataset array called |ds|, and define |Tomato|, |Soil|,
% and |Fertilizer| as categorical variables. 
ds = fertilizer;
ds.Tomato = nominal(ds.Tomato);
ds.Soil = nominal(ds.Soil);
ds.Fertilizer = nominal(ds.Fertilizer);  

%% 
% Fit a linear mixed-effects model, where |Fertilizer| and |Tomato| are
% the fixed-effects variables, and the mean yield varies by the block (soil
% type) and the plots within blocks (tomato types within soil types) independently. 
%
% This model corresponds to
%
% $$\begin{array}{l}
% {y_{imjk}} = {\beta _0} + \sum\limits_{m = 2}^4 {{\beta _{1m}}I{{\left[ F \right]}_{im}}}  + \sum\limits_{j = 2}^5 {{\beta _{2j}}I{{\left[ T \right]}_{ij}}}  + \sum\limits_{j = 2}^5 {\sum\limits_{m = 2}^4 {{\beta _{3mj}}I{{\left[ F \right]}_{im}}I{{\left[ T \right]}_{ij}}} } \\
% \quad \quad \quad  + {b_{0k}}{S_k} + {b_{0jk}}{(S*T)_{jk}} + {\varepsilon _{imjk}},
% \end{array}$$
%
% where $i$ = 1, 2, ..., 60, index $m$ corresponds to the fertilizer types,
% $j$ corresponds to the tomato types, and $k$ = 1, 2, 3 corresponds to
% the blocks (soil). $S_{k}$ represents the $k$ th soil
% type, and $(S*T)_{jk}$ represents the $j$ th tomato
% type nested in the $k$ th soil type. $I[F]_{im}$ is
% the dummy variable representing level $m$ of the fertilizer. Similarly,
% $I[T]_{ij}$ is the dummy variable representing level
% $j$ of the tomato type. 
%
% The random effects and observation error have these prior distributions:
% $b_{0k}$~N(0, $\sigma^{2}_{S}$ ),
% $b_{0jk}$~N(0, $\sigma^{2}_{S*T}$ ),
% and $\epsilon_{imjk}$ ~ N(0, $\sigma^{2}$ ). 
lme = fitlme(ds,'Yield ~ Fertilizer * Tomato + (1|Soil) + (1|Soil:Tomato)') 

%%
% The $p$-values corresponding to the last 12 rows in the fixed-effects
% coefficients display (0.82701 to 0.47881) indicate that interaction coefficients
% between the tomato and fertilizer types are not significant. To test for
% the overall interaction between tomato and fertilizer, use the |<docid:stats_ug.btvrk5h
% anova>| method after refitting the model using |'effects'| contrasts. 

%%
% The confidence interval for the standard deviations of the random-effects
% terms ( $\sigma^{2}_{S}$ ), where the intercept
% is grouped by soil, is very large. This term does not appear significant.  

%% 
% Refit the model after removing the interaction term |Tomato:Fertilizer|
% and the random-effects term |(1 | Soil)|. 
lme = fitlme(ds,'Yield ~ Fertilizer + Tomato + (1|Soil:Tomato)') 

%%
% You can compare the two models using the |<docid:stats_ug.btvrk6m compare>|
% method with the simulated likelihood ratio test since both a fixed-effect
% and a random-effect term are tested.