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.