将混合模型的结果保存在数据集中
Save the results of a mixed model in a dataset
我正在拟合下面的混合模型:
. mixed y trt || clst:trt, nocons reml dfmethod(sat)
Performing EM optimization:
Performing gradient-based optimization:
Iteration 0: log restricted-likelihood = -1295.3123
Iteration 1: log restricted-likelihood = -1295.3098
Iteration 2: log restricted-likelihood = -1295.3098
Computing standard errors:
Computing degrees of freedom:
Mixed-effects REML regression Number of obs = 919
Group variable: clst Number of groups = 49
Obs per group:
min = 1
avg = 18.8
max = 30
DF method: Satterthwaite DF: min = 888.00
avg = 900.91
max = 913.83
F(1, 913.83) = 0.40
Log restricted-likelihood = -1295.3098 Prob > F = 0.5251
------------------------------------------------------------------------------
y | Coef. Std. Err. t P>|t| [95% Conf. Interval]
-------------+----------------------------------------------------------------
trt | .1455914 .2290005 0.64 0.525 -.3038366 .5950193
_cons | .3951269 .2241477 1.76 0.078 -.0447941 .835048
------------------------------------------------------------------------------
------------------------------------------------------------------------------
Random-effects Parameters | Estimate Std. Err. [95% Conf. Interval]
-----------------------------+------------------------------------------------
clst: Identity |
var(trt) | .0341507 .0173905 .0125877 .092652
-----------------------------+------------------------------------------------
var(Residual) | .9546016 .0453034 .8698131 1.047655
------------------------------------------------------------------------------
LR test vs. linear model: chibar2(01) = 9.46 Prob >= chibar2 = 0.0010
. return list
scalars:
r(level) = 95
matrices:
r(table) : 9 x 4
接下来,我计算ICC
如下:
. nlcom (icc_est: (exp(_b[lns1_1_1:_cons])^2)/((exp(_b[lns1_1_1:_cons])^2)+(exp(_b[lnsig_e:_cons])^2)))
icc_est: (exp(_b[lns1_1_1:_cons])^2)/((exp(_b[lns1_1_1:_cons])^2)+(exp(_b[lnsig_e:_cons])^2))
------------------------------------------------------------------------------
y | Coef. Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
icc_est | .0345392 .0171907 2.01 0.045 .0008461 .0682323
------------------------------------------------------------------------------
如何将结果保存到数据集中?
我想保留显示的所有三个表:固定效应、随机效应和 ICC 结果。
考虑以下使用 Stata 的 pig
玩具数据集的可重现示例:
webuse pig, clear
mixed weight week || id:week, nocons reml dfmethod(sat)
nlcom (icc_est: (exp(_b[lns1_1_1:_cons])^2)/((exp(_b[lns1_1_1:_cons])^2)+(exp(_b[lnsig_e:_cons])^2))), post
------------------------------------------------------------------------------
weight | Coef. Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
icc_est | .1380299 .0265754 5.19 0.000 .0859431 .1901167
------------------------------------------------------------------------------
以下对我有用:
generate double coef = _b[icc_est]
generate double se = _se[icc_est]
generate p = string(2 * (normal(-(_b[icc_est] / _se[icc_est]))), "%9.3f")
generate double upper = _b[icc_est] + _se[icc_est] * invnormal(0.025)
generate double lower = _b[icc_est] + _se[icc_est] * invnormal(0.975)
list coef se p upper lower in 1
+-------------------------------------------------------+
| coef se p upper lower |
|-------------------------------------------------------|
1. | .13802987 .02657538 0.000 .08594308 .19011667 |
+-------------------------------------------------------+
save mydata.dta
主模型的结果过程类似。
作为后续,轻松获取随机截距方差和SE以及残差方差和SE将需要多一行代码。但是正如之前的回复所指出的,主模型的结果是通过与ICC结果相同的方式获得的。请参阅下面的代码。
mixed y trt || clst:trt, nocons reml dfmethod(sat)
gen double fixedcoef = _b[trt]
gen double fixedse = _se[trt]
_diparm lns1_1_1, f(exp(@)^2) d(2*exp(@)^2)
gen double randomcoef = r(est)
gen double randomse = r(se)
_diparm lnsig_e, f(exp(@)^2) d(2*exp(@)^2)
gen double residcoef = r(est)
gen double residse = r(se)
我正在拟合下面的混合模型:
. mixed y trt || clst:trt, nocons reml dfmethod(sat)
Performing EM optimization:
Performing gradient-based optimization:
Iteration 0: log restricted-likelihood = -1295.3123
Iteration 1: log restricted-likelihood = -1295.3098
Iteration 2: log restricted-likelihood = -1295.3098
Computing standard errors:
Computing degrees of freedom:
Mixed-effects REML regression Number of obs = 919
Group variable: clst Number of groups = 49
Obs per group:
min = 1
avg = 18.8
max = 30
DF method: Satterthwaite DF: min = 888.00
avg = 900.91
max = 913.83
F(1, 913.83) = 0.40
Log restricted-likelihood = -1295.3098 Prob > F = 0.5251
------------------------------------------------------------------------------
y | Coef. Std. Err. t P>|t| [95% Conf. Interval]
-------------+----------------------------------------------------------------
trt | .1455914 .2290005 0.64 0.525 -.3038366 .5950193
_cons | .3951269 .2241477 1.76 0.078 -.0447941 .835048
------------------------------------------------------------------------------
------------------------------------------------------------------------------
Random-effects Parameters | Estimate Std. Err. [95% Conf. Interval]
-----------------------------+------------------------------------------------
clst: Identity |
var(trt) | .0341507 .0173905 .0125877 .092652
-----------------------------+------------------------------------------------
var(Residual) | .9546016 .0453034 .8698131 1.047655
------------------------------------------------------------------------------
LR test vs. linear model: chibar2(01) = 9.46 Prob >= chibar2 = 0.0010
. return list
scalars:
r(level) = 95
matrices:
r(table) : 9 x 4
接下来,我计算ICC
如下:
. nlcom (icc_est: (exp(_b[lns1_1_1:_cons])^2)/((exp(_b[lns1_1_1:_cons])^2)+(exp(_b[lnsig_e:_cons])^2)))
icc_est: (exp(_b[lns1_1_1:_cons])^2)/((exp(_b[lns1_1_1:_cons])^2)+(exp(_b[lnsig_e:_cons])^2))
------------------------------------------------------------------------------
y | Coef. Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
icc_est | .0345392 .0171907 2.01 0.045 .0008461 .0682323
------------------------------------------------------------------------------
如何将结果保存到数据集中?
我想保留显示的所有三个表:固定效应、随机效应和 ICC 结果。
考虑以下使用 Stata 的 pig
玩具数据集的可重现示例:
webuse pig, clear
mixed weight week || id:week, nocons reml dfmethod(sat)
nlcom (icc_est: (exp(_b[lns1_1_1:_cons])^2)/((exp(_b[lns1_1_1:_cons])^2)+(exp(_b[lnsig_e:_cons])^2))), post
------------------------------------------------------------------------------
weight | Coef. Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
icc_est | .1380299 .0265754 5.19 0.000 .0859431 .1901167
------------------------------------------------------------------------------
以下对我有用:
generate double coef = _b[icc_est]
generate double se = _se[icc_est]
generate p = string(2 * (normal(-(_b[icc_est] / _se[icc_est]))), "%9.3f")
generate double upper = _b[icc_est] + _se[icc_est] * invnormal(0.025)
generate double lower = _b[icc_est] + _se[icc_est] * invnormal(0.975)
list coef se p upper lower in 1
+-------------------------------------------------------+
| coef se p upper lower |
|-------------------------------------------------------|
1. | .13802987 .02657538 0.000 .08594308 .19011667 |
+-------------------------------------------------------+
save mydata.dta
主模型的结果过程类似。
作为后续,轻松获取随机截距方差和SE以及残差方差和SE将需要多一行代码。但是正如之前的回复所指出的,主模型的结果是通过与ICC结果相同的方式获得的。请参阅下面的代码。
mixed y trt || clst:trt, nocons reml dfmethod(sat)
gen double fixedcoef = _b[trt]
gen double fixedse = _se[trt]
_diparm lns1_1_1, f(exp(@)^2) d(2*exp(@)^2)
gen double randomcoef = r(est)
gen double randomse = r(se)
_diparm lnsig_e, f(exp(@)^2) d(2*exp(@)^2)
gen double residcoef = r(est)
gen double residse = r(se)