Re: [Gretl-users] Multiple Imputation analysis
On Fri, 5 Mar 2010, Nathan Paxton wrote:
I do work with multiple imputation datasets, and I generate them
in R, using the Amelia package. But the sorts of tests I want to
run aren't easily done using any of the more canned packages, so
I end up having to write long code files to implement Rubin's
(1986) combination rules for MI parameter estimates.
I don't need to use gretl to do the actual MI, but I wonder if
there's any facility in gretl for doing analysis on separate
MI-generated datasets and then combining the results according
to the Rubin rules. That is, can gretl be told easily to run a
procedure m times for the m existing datasets and then combine
the results to get 1 parameter estimate?
One approach would be to store each dataset as a pair of matrices,
one for the dependent variable and one for the independent
variables: say y1.mat, X1.mat; y2.mat, X2.mat; and so on up to
ym.mat and Xm.mat. (See the help for mwrite and mread.)
Then you could script the analysis along these lines:
<script>
matrix y X beta
matrix Beta = {}
matrix W = {}
matrix V = {}
# loop, appending a row of OLS coeffs and a
# row of variances to the matrices Beta and W
# respectively, at each step
loop i=1..m
y = mread("y$i.mat")
X = mread("X$i.mat")
beta = mols(y, X, null, &V)
Beta |= beta'
W |= diag(V)'
endloop
# means of coefficients
bbar = meanc(Beta)
print bbar
# means of within-imputation variances
Wbar = meanc(W)
print Wbar
# deviations of coefficients from means across imputations
Bdev = cdemean(Beta)
# Between-imputation variance
B = sumc(Bdev .* Bdev) / (m - 1)
# overall variance
T = Wbar + (1 + 1/m) * B
# overall standard errors
se = sqrt(T)
se
</script>
An alternative to loading the datasets as matrices would be to use
the "append" command within the loop to grab the data as regular
series. If the names of the variables are the same in all the
datasets, and they each have an initial "obs" column that
identifies the observations (e.g. "A", "B", "C",...) then
"append"
will in fact overwrite the old values with the new.
Trivial example of suitable data for this trick:
1.txt:
obs v1 v2
A 1 6
B 2 7
C 3 8
D 4 9
E 5 10
2.txt:
A 6 7
B 7 8
C 8 9
D 9 10
E 10 11
Allin Cottrell