R中缺失和审查数据的多重插补

Question

R中缺失和审查数据的多重插补

我有一个同时包含随机缺失（MAR）和审查数据的数据集。这些变量是相关的，因此我尝试有条件地估算缺失的数据，以便可以估计相关的多元正态分布的分布参数。我想使用Gibbs MCMC，但是很难执行该程序。我的数据框有5个变量（表示为x1：x5），1099个样本，其中包含MAR，检查值和观察值的某种组合。到目前为止，这是我尝试过的：

# packages
library(msm, tmvtnorm, MCMCpack)

# priors 
theta0<-c(rep(0, 5))
Sigma0<-S0<-diag(5)  
nu0<-4 

# initialize parameters
theta<-c(rep(0, 5))
Tau<-diag(5) 

# initialize output matrix
n_samples <- 1000
mu_MCMC <- matrix(0, nrow = n_samples, ncol = 5)
mu_MCMC[1,] <- theta
cov_MCMC <- matrix(0, nrow = n_samples, ncol = 25)
cov_MCMC[1,] <- c(diag(5))

# detection limits
det_lim <- matrix(c(-1.7, 0, 0, 0, 0), nrow = 1, ncol = 5)

# function to detect NaN (i.e., below detection data)
is.nan.data.frame <- function(x)
    do.call(cbind, lapply(x, is.nan))

for(i in 2:n_samples){
    imputedDF <- data.frame()
    for(r in 1:nrow(originalDF)){
        # variables that are MAR or censored
        mis <- r[, is.na(r) & is.nan(r)]    
        # variables that are observed
        obs <- r[, !is.na(r)]

        # subset mu for missing, observed
        mu1 <- mu[, names(r) %in% names(mis)]
        mu2 <- mu[, names(r) %in% names(obs)]

        # calculate sigmas for MVN partitions of mis, obs
        sigma11 <- sigma[names(r) %in% names(mis), names(r) %in% names(mis)]
        sigma22 <- sigma[names(r) %in% names(obs), names(r) %in% names(obs)]
        sigma12 <- sigma[names(r) %in% names(obs), names(r) %in% names(mis)]
        sigma21 <- t(sigma12)

        # create matrix for detection limits based on missing values
        ## if NaN, use detection limit; if NA use Inf
        dl <- c(ifelse("x1" %in% names(is.nan(r)), det_lim[1, "x1"], Inf), 
                ifelse("x2" %in% names(is.nan(r)), det_lim[1, "x2"], Inf), 
                ifelse("x3" %in% names(is.nan(r)), det_lim[1, "x3"], Inf), 
                ifelse("x4" %in% names(is.nan(r)), det_lim[1, "x4"], Inf), 
                ifelse("x5" %in% names(is.nan(r)), det_lim[1, "x5"], Inf))

        # compute mu, sigma to use for conditional MVN
        ## if all values are missing
        if(length(names(obs) == 0) {
            mu_mis <- mu1
            sigma_mis <- sigma11
        ## otherwise
            } else {
                mu_mis <- mu1 + sigma12 %*% solve(sigma22) * (obs - t(mu2))
                sigma_mis <- sigma11 - sigma12 %*% solve(sigma22) %*% sigma21
        }

        # imputation
        ## if all data are observed, missing is empty
        if(length(obs) == 0) {
            mis_impute <- data.frame()
        ## only need to impute a single value
            } else if(length(names(mis)) == 1) {       
                  mis_impute <- rtnorm(1, mean = mu_mis, sd = sigma_mis, lower = -Inf, upper = dl)
        ## have more than one missing value         
                  } else {
                      mis_impute <- rtmvnorm(1, mean = mu_mis, sigma = sigma_mis, lower = rep(-Inf, length = length(names(mis))), upper = dl)
        }

       # merge observed values with simulated 
       ## if all values observed   
       if(length(names(mis)) == 0) {
           sim_result <- obs
           } else {
                 sim_result <- cbind(mis_impute, obs) 
       }

       imputedDF <- rbind(imputedDF, sim_result)
    }

    # update theta
    v <- solve(solve(Sigma0) + nrow(sim_result)*Tau)
    m <- v %*% (solve(Sigma0) %*% theta0 + Tau %*% apply(sim_result,2,sum))
    mu <- as.data.frame(rmvnorm(1,m,v))
    mu_MCMC[i,] <- mu

    # update Sigma
    tmp <- t(sim_result) - mu
    Tau <- rwish(nu0 + nrow(sim_result), solve(S0 + t(tmp) %*% tmp)) 
    sigma <- matrix(c(solve(Tau)), nrow = 5, ncol = 5, byrow = TRUE)
    cov_MCMC[i,] <- c(solve(Tau))
}

归档时间：	9 年，1 月前
查看次数：	706 次
最近记录：	6 年，11 月前

R中缺失和审查数据的多重插补

前言：

归因

生成示例数据集。

资料摘要

数据是完全随机丢失还是不是随机丢失？

控制台输出

数据估算

安慰：

安慰：

MCMC回归

控制台输出

引文：