www.gusucode.com > stats 源码程序 matlab案例代码 > stats/RandomInterceptModelExample.m
%% Fit Random Intercept LME Model %% % Load the sample data. load flu %% % The |flu| dataset array has a |Date| variable, and 10 variables containing % estimated influenza rates (in 9 different regions, estimated from Google(R) % searches, plus a nationwide estimate from the Centers for Disease Control % and Prevention, CDC). %% % To fit a linear-mixed effects model, your data must be in a properly formatted % dataset array. To fit a linear mixed-effects model with the influenza % rates as the responses, combine the nine columns corresponding to the % regions into an array. The new dataset array, |flu2|, must have the % new response variable |FluRate|, the nominal variable |Region| that shows % which region each estimate is from, the nationwide estimate |WtdILI|, % and the grouping variable |Date|. flu2 = stack(flu,2:10,'NewDataVarName','FluRate',... 'IndVarName','Region'); flu2.Date = nominal(flu2.Date); %% % Display the first six rows of |flu2|. flu2(1:6,:) %% % Fit a linear mixed-effects model with a fixed-effects term for the nationwide % estimate, |WtdILI|, and a random intercept that varies by |Date|. The % model corresponds to % % $${y_{im}} = {\beta _0} + {\beta _1}WtdIL{I_{im}} + {b_{0m}} + % {\varepsilon _{im}},\quad i = 1,2,...,468,\quad m = 1,2,...,52,$$ % % where $y_{im}$ is the observation $i$ for level $m$ % of grouping variable |Date|. $b_{0m}$ is the random % effect for level $m$ of the grouping variable |Date| and $\epsilon_{im}$ % is the observation error for observation $i$. The random effect has the % prior distribution, $b$ ~ N(0, $\sigma^{2}_{b}$ ) % and the error term has the distribution, $\epsilon$ ~ N(0, $\sigma^{2}$ ). lme = fitlme(flu2,'FluRate ~ 1 + WtdILI + (1|Date)') %% % Estimated covariance parameters are displayed in the section titled "Random % effects covariance parameters". The estimated value of $\sigma_{b}$ % is 0.17146 and its 95% confidence interval is [0.13227, 0.22226]. Since % this interval does not include 0, the random-effects term is significant. % You can formally test the significance of any random-effects term using % a likelihood ratio test via the |<docid:stats_ug.btvrk6m compare>| method. %% % The estimated response at an observation is the sum of the fixed effects % and the random-effect value at the grouping variable level corresponding % to that observation. For example, the estimated flu rate for observation % 28 is % % $$\begin{array}{l} % {{\hat y}_{28}} = {{\hat \beta }_0} + {{\hat \beta }_1}WtdIL{I_{28}} + {{\hat b}_{10/30/2005}}\\ % \quad \quad = 0.1639 + 0.7236*(1.343) + 0.3318\\ % \quad \quad = 1.46749, % \end{array}$$ % % where ${\hat b}$ is the estimated % best linear unbiased predictor (BLUP) of the random effects for the intercept. % You can compute this value as follows. beta = fixedEffects(lme); [~,~,STATS] = randomEffects(lme); % Compute the random-effects statistics (STATS) STATS.Level = nominal(STATS.Level); y_hat = beta(1) + beta(2)*flu2.WtdILI(28) + STATS.Estimate(STATS.Level=='10/30/2005') %% % You can display the fitted value using the |fitted| method. F = fitted(lme); F(28)