I am using mi estimate: mixed, with reml. In Stata 16.1
Following estimation of the model, I am needing to estimate both the random effect and it's standard error manually. mi doesn't generate these with the predict command.
However, the 'manual' approach I have successfully used previously to estimate the standard error of the random effect (after a mle model), doesn't generate the same result with reml.
See the example below. This is not with mi data, but identifies the problem.
Code:
webuse pig, clear mixed weight week || id:, reml gen IN_MODEL_REML = e(sample) predict double (RE_REML) if IN_MODEL_REML == 1, reffects reses(RE_REML_SE) predict double PRED_XB_REML, xb gen double RESID_PRED_REML = weight - PRED_XB_REML sort id by id: egen double MEAN_RESID_REML = mean(RESID_PRED_REML) by id: egen CLUSTER_OBS_REML = total(IN_MODEL_REML) , missing quietly mixed weight week || id:, reml predictnl double RE_REML_MAN = (exp([lns1_1_1]_cons)^2/(exp([lns1_1_1]_cons)^2 + exp([lnsig_e]_cons)^2/CLUSTER_OBS_REML))*MEAN_RESID_REML // Predict RE and between imputation RE variance. Formula from Rabe-Hesketh et al (2012), p. 112 predictnl double RE_VAR_REML_MAN = (1-(exp([lns1_1_1]_cons)^2/(exp([lns1_1_1]_cons)^2 + exp([lnsig_e]_cons)^2/(CLUSTER_OBS_REML)))) *exp([lns1_1_1]_cons)^2 //estimate within RE variance . Formula from Rabe-Hesketh et al (2012), p. 113. gen double RE_SE_REML_MAN = RE_VAR_REML_MAN^0.5 by id: gen OB_NUM = _n sum RE_REML RE_REML_MAN if OB_NUM == 1 sum RE_REML_SE RE_SE_REML_MAN if OB_NUM == 1
Standard errors for BLUPs are calculated based on the iterative technique of Bates and Pinheiro (1998, sec. 3.3) for estimating the BLUPs themselves. If estimation is done by REML, these standard errors account for uncertainty in the estimate of , while for ML the standard errors treat as known. As such, standard errors of REML-based BLUPs will usually be larger.
If not I may need to hit the text books to better understand how REML works. Using REML was not my choice...I am needing to replicate other work undertaken in R.
Regards,
Andrew
0 Response to Manually estimating standard error of the random effect, after the mi estimate: mixed command, with reml
Post a Comment