Missing value

[Duoduo]

Hi, everyone. I tried to use only non-NA y values (the response has several missing values, but we want to skip those index). I changed the index in Jags code.  It can get the MCMC result with no error. However, the estimated parameter is very different than the true parameter I set. In my opinion, the method that directly see the missing response as latent variable is similar to the method that skip the index of missing response(we did right now, only consider observed response in the MCMC). I give the code and data below and I am not sure if I made any mistake here.  Thanks!

library(data.table)

library(scales)
library(MCMCvis)
library(nimble)
library(MASS)
library(writexl)
library(readxl)
library(grid)
library(dplyr)

set.seed(1994)

theta=c()

for (n in  1:1000) {

  theta[n]=rnorm(1,0,1)
}
alpha=rep(1.5,7)
lamda=c(-1.5,-1,-0.5,0,0.5,1,1.5)
beta_zeta <- as.data.frame(mvrnorm(n=30,
                                   mu=c(-2.197, 4.394),
                                   Sigma=matrix(c(1, -0.8, -0.8, 1), ncol=2)))

parameter_DINA_alpha_lamda<-cbind(alpha,lamda)
parameter_DINA_beta_zeta<-beta_zeta
colnames(parameter_DINA_beta_zeta) <- c('beta','zeta')
parameter_DINA_theta<-as.data.frame(theta)
parameter_DINA_alpha_lamda<-as.data.frame(parameter_DINA_alpha_lamda)

lamdaD<-lamda
thetaD<-theta
betaD<-parameter_DINA_beta_zeta$beta
zetaD<-parameter_DINA_beta_zeta$zeta
gammaD<-alpha

y=fread("D:/try_new/response_ex.csv")
setDF(y)
#good method
#DINA model
numNonNAs <- apply(y, 1, function(x) sum(!is.na(x)))
nonNAindices <- array(NA, dim(y))
for(i in 1:nrow(y))   nonNAindices[i,1:numNonNAs[i]] <- which(!is.na(y[i,]))

#model we fit
code=nimbleCode({
  for (n in 1:1000) {

    for(i in 1:numNonNAs[n]) {

      for (k in 1:7) {w[n, nonNAindices[n,i], k] <- pow(attribute[n, k], Q[nonNAindices[n,i], k])}
      logit(prob[n, nonNAindices[n,i]]) <- beta[nonNAindices[n,i]] + delta[nonNAindices[n,i]] * prod(w[n, nonNAindices[n,i], 1:7])
      y[n, nonNAindices[n,i]] ~ dbern(prob[n, nonNAindices[n,i]])   ## you'll need to complete this

     
    }}
  for (n in 1:1000) {
    for (k in 1:7) {
      logit(att_prob[n, k]) <- gamma[k] * theta[n] - lambda[k]
      attribute[n, k] ~ dbern(att_prob[n, k])
    }}
  for (n in 1:1000) {theta[n]~ dnorm(0,1)}
  for (i in 1:30) {
    item_parameter[i, 1:2] ~ dmnorm(item_mu[1:2], item_den[1:2, 1:2])
    beta[i] <- item_parameter[i, 1]
    delta[i] <- item_parameter[i, 2]
  }
  for (k in 1:7) {
    lambda[k] ~ dnorm(0, 1)
    gamma[k] ~ T(dnorm(0, 1),0,)
  }
  item_mu[1] ~ dnorm(0, 1)
  item_mu[2] ~ T(dnorm(0, 1),0,)
  R[1, 1] <- 1
  R[2, 2] <- 1
  R[1, 2] <- 0
  R[2, 1] <- 0
  item_den[1:2, 1:2] ~ dwish(R[1:2, 1:2], 2)
})

data <- list(y=y,Q=Q)

inits <- list(item_den = diag(2))

Consts <- list(numNonNAs = numNonNAs,
               nonNAindices=nonNAindices)

monitors=c("gamma","lambda","beta","delta")

model<-nimbleModel(code,data = data,inits=inits,constants = Consts)

cModel_DINA <-compileNimble(model)
conf_DINA <- configureMCMC(model, monitors = monitors)
model_DINAMCMC <- buildMCMC(conf_DINA)
cModel_DINAMCMC <- compileNimble(model_DINAMCMC, project = model)
system.time(samples_DINA<-runMCMC(cModel_DINAMCMC,

                                  nchains=2,niter = 80000, nburnin = 10000))

[Chris Paciorek]

Without diving deep into the statistical formulation of the model, my only thought here is that there might be identifiability issues with your model.

It may be that various parameter combinations could about equally well fit the data. 

-Chris

--
You received this message because you are subscribed to the Google Groups "nimble-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to nimble-users...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/e0afdd55-6d6f-459b-a5ee-60204b5057a5n%40googlegroups.com.

[Duoduo]

Hi,  Chris. Thanks for the response! I am wondering if you could explain more about the comment you made: It may be that various parameter combinations could about equally well fit the data?  If I include the NA in response then to run the model, nimble will directly see the missing part as a latent variable (this can run). But now I prefer to only consider the non missing part of response in the model, I need to skip the missing index, for example, if the Y12 is missing value, then Y12 will not be included in the model. But the result is not what I expect. So I am not clear if I made any mistake of nimble code here. Thanks!

[Chris Paciorek]

Let's back up so I can better understand what you are seeing.

I think you are saying that when you fit the model where you skip the missing values you get a result that is not similar to when you treat the missing values as latent variables. Is that correct?

Also, you mention JAGS code. Are you getting the surprising result in JAGS or NIMBLE or both?

I don't see anything that looks wrong in terms of how you set up the model that considers only the non-missing values. 

Regarding non-identifiability, I just meant that it may be that different combinations of parameters can reasonably well fit the data. If that's the case, it may be that the MCMC for the full data model is drawing samples in part of the parameter space and the MCMC for the observed data model is drawing samples in another part of the space, and the mixing is not good enough for the two MCMCs to fully move around the space.

One strategy you could try is to simulate a dataset where you know the true parameter values and then do MCMC under both of the approaches and see whether one or both do a good job of recovering the true parameter values.

-chris

To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/9f311290-9512-4bb9-870c-d53d4a25abc3n%40googlegroups.com.