#simulate a partially observed dataset
set.seed(1234)
n <- 100
x <- rnorm(n)
y <- x+rnorm(n)
y[1:50] <- NA
temp <- data.frame(x,y)
#impute using normUniImp, without posterior draws
imps <- normUniImp(temp, y~x, M=10, pd=FALSE)
#define a function which performs our desired analysis on a dataset, returning
#the parameter estimates
yonx <- function(inputData) {
fitmod <- lm(y~x, data=inputData)
list(est=c(fitmod$coef,sigma(fitmod)^2))
}
#define a function which when passed a dataset and parameter
#vector, calculates the likelihood score vector
myScore <- function(inputData, parm) {
beta0 <- parm[1]
beta1 <- parm[2]
sigmasq <- parm[3]
res <- inputData$y - beta0 - beta1*inputData$x
cbind(res/sigmasq, (res*inputData$x)/sigmasq, res^2/(2*sigmasq^2)-1/(2*sigmasq))
}
#call scoreBased to perform variance estimation
scoreBased(imps, analysisFun=yonx, scoreFun=myScore)
Run the code above in your browser using DataLab