Hi Folks,
I'm using the 'norm' package (based on Shafer's NORM)
on some  data. In outline, (X,Y) are bivariate normal,
var(X)=0.29, var(Y)=24.4, cov(X,Y)=-0.277,
there are some 900 cases, and some 170 values of Y
have been set "missing" (NA).

The puzzle is that, repeating the multiple imputation
starting from the same random seed, I get different
answers from the repeats depending if I do an odd number
of imputations, but the same answer on the repeats
if I do en even number (which includes the second
repeat of an odd number).

It may possibly have something to do with how I've
written the code for the loop, but if so then I'm
not seeing it!

CODE:

## Set up the situation:
Data<-read.csv("MyData.csv")
X<-Data$X; Y<-Data$Y
##(If you want to try it, set your own data here)
Raw<-cbind(X,Y)
library(norm)

## Initialise stuff
s<-prelim.norm(Raw)
t0<-em.norm(s)

##########################
## Set the Random Seed
rngseed(31425)

## Do the first imputation:
t     <-  da.norm(s,t0,steps=20)
Imp   <- imp.norm(s,t, Raw)
X.Imp <- Imp[,1]; Y.Imp<-Imp[,2]

## Now do the rest, and accumulate lists of the results
## Est.Imp = list of estimated coeffs
## SE.Imp  = list of SEs of estimated coeffs:
Est.Imp <- list(summary(lm(Y.Imp~X.Imp))$coef[,1])
SE.Imp  <- list(summary(lm(Y.Imp~X.Imp))$coef[,2])
N=4
for(i in (2:N)){
  t<-da.norm(s,t,steps=20)
  Imp<-imp.norm(s,t,Raw)
  X.Imp<-Imp[,1]; Y.Imp<-Imp[,2]
  Est.Imp<-c(Est.Imp,list(summary(lm(Y.Imp~X.Imp))$coef[,1]))
  SE.Imp <-c( SE.Imp,list(summary(lm(Y.Imp~X.Imp))$coef[,2]))
}

## Finally, combine the imputations:
mi.inference(Est.Imp,SE.Imp)


I've illustrated N=4 (even) above, for 4 imputations.

Now, I run the code repeatedly from "## Set the Random Seed"
down to "mi.inference(Est.Imp,SE.Imp)"

With N=4, I always get the same result.

If I set N=5, I alternately get different results:
The second run is different from the first, but the
third is the same as the first, and the fourth is the
same as the second, ...

In general, for even N, it is as for N=4, and for odd N
it is as for N=5.

Any ideas???

Thanks,
Ted.

--------------------------------------------------------------------
E-Mail: (Ted Harding) <[EMAIL PROTECTED]>
Fax-to-email: +44 (0)870 094 0861
Date: 21-Sep-07                                       Time: 14:13:27
------------------------------ XFMail ------------------------------

______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Reply via email to