missing values in two response (joint model)

Duoduo

unread,

Oct 1, 2023, 8:38:38 PM10/1/23

to nimble-users

Hi, everyone. When I tried to use two response (y and logT below) with missing values in both, if the iteration is small like 1000, I can get the result I want. But if the iteration is bigger like 50000, some parameters show NA in the result. I am wondering if there is any solution for this? Thanks.

code=nimbleCode({
for (n in 1:N) {
for (i in 1:I) {
logit(prob[n, i]) <- beta[i] + delta[i] * prod(w[n, i])
y[n, i] ~ dbern(prob[n, i])
logT[n, i] ~ dnorm(zeta[i] + theta[n], den_epsilon[i])
}}
for (n in 1:N) {theta[n] ~ dnorm(0,1)}
for (i in 1:I) {
item_parameter[i, 1:3] ~ dmnorm(item_mu[1:3], item_den[1:3, 1:3])
beta[i] <- item_parameter[i, 1]
delta[i] <- item_parameter[i, 2]
zeta[i] <- item_parameter[i, 3]
den_epsilon[i] ~ dgamma(1, 1)
}
item_mu[1] ~ dnorm(0, 1)
item_mu[2] ~ T(dnorm(0, 1),0,)
item_mu[3] ~ dnorm(0,1)
item_den[1:3, 1:3] ~ dwish(R[1:3, 1:3], 3)
})

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

monitors=c("beta","delta","zeta"," den_epsilon","y","logT")

Daniel Turek

unread,

Oct 2, 2023, 6:05:24 AM10/2/23

to Duoduo, nimble-users

Duoduo, are you able to provide a reproducible script that demonstrates this behavior? That would include the data (or proxy data), and the calls you used to build the model, and run the MCMC.

Thank you,

Daniel

--
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/16617ee3-e387-47af-a4c1-84e1da0a8d6cn%40googlegroups.com.

Duoduo

unread,

Oct 4, 2023, 5:29:08 PM10/4/23

to nimble-users

Hi, Dan. The data I have is attached. I am wondering if I can delete those NA case from response manually(only consider the non-NA) case. But I always have the errors and I am not sure how to fix that. THanks!!

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

library(data.table)

y_NA=fread("missing_response.csv")
setDF(y_NA)
vec1=colSums(!is.na(t(y_NA)))
vec1=c(0,vec1)
vec1=as.vector(vec1)

loc=which(!is.na(y_NA),arr.ind = T)
loc=data.frame(loc)
colnames(loc)=c("person","item")
try=loc %>%
group_by(person) %>%
summarise(score = list(item))

trytry=try[[2]]
vec2=as.vector(unlist(trytry))

y=c(t(y_NA))
y=na.omit(y)

#model we fit
code=nimbleCode({
for (n in 1:1000) {
for (i in 1:vec1[n+1]) {
add=vec1[n]
index=vec2[add+i]
for (k in 1:7) {w[n, index, k] <- pow(attribute[n, k], Q[index, k])}
logit(prob[1+(i-1)+vec1[n]]) <- beta[index] + delta[index] * prod(w[n,index, 1:7])
y[1+(i-1)+vec1[n]] ~ dbern(prob[1+(i-1)+vec1[n]])
}}
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,vec1=vec1,vec2=vec2)
inits <- list(item_den = diag(2))
monitors=c("gamma","lambda","beta","delta")
model<-nimbleModel(code,data = data,inits=inits)
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 = 100, nburnin = 10))

missing_response.csv

Daniel Turek

unread,

Oct 6, 2023, 6:41:29 AM10/6/23

to Duoduo, nimble-users

Duoduo, I gave this a quick go, and it seems to not be able to find the "Q" variable. Is there something I overlooked, or would you mind trying to make it fully reproducible?

Thanks,

Daniel

To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/50bbe2a4-edf7-4d8a-b59f-a596b904ea2cn%40googlegroups.com.

Duoduo

unread,

Oct 6, 2023, 8:29:57 AM10/6/23

to nimble-users

Hi, Dan. Sorry! The Q matrix is below and I am trying to find some way to only consider observed response in model and skip the NA index, Thanks!

Q=data.frame(matrix(ncol = 7,nrow = 30))
Q[1,]=c(1,0,0,0,0,0,0)
Q[2,]=c(0,1,0,0,0,0,0)
Q[3,]=c(0,0,1,0,0,0,0)
Q[4,]=c(0,0,0,1,0,0,0)
Q[5,]=c(0,0,0,0,1,0,0)
Q[6,]=c(0,0,0,0,0,1,0)
Q[7,]=c(0,0,0,0,0,0,1)
Q[8,]=c(1,0,0,0,0,0,0)
Q[9,]=c(0,1,0,0,0,0,0)
Q[10,]=c(0,0,1,0,0,0,0)
Q[11,]=c(0,0,0,1,0,0,0)
Q[12,]=c(0,0,0,0,1,0,0)
Q[13,]=c(0,0,0,0,0,1,0)
Q[14,]=c(0,0,0,0,0,0,1)
Q[15,]=c(1,1,0,0,0,0,0)
Q[16,]=c(0,0,1,1,0,0,0)
Q[17,]=c(0,0,0,1,1,0,0)
Q[18,]=c(0,0,1,0,0,1,1)
Q[19,]=c(0,0,0,1,1,1,0)
Q[20,]=c(0,1,0,0,1,1,0)
Q[21,]=c(0,0,1,1,0,0,1)
Q[22,]=c(1,0,1,0,0,0,1)
Q[23,]=c(0,1,1,1,0,1,0)
Q[24,]=c(1,1,0,0,1,1,0)
Q[25,]=c(1,1,1,1,0,0,0)
Q[26,]=c(0,0,0,1,1,1,1)
Q[27,]=c(1,1,0,0,1,0,1)
Q[28,]=c(1,1,1,0,0,1,1)
Q[29,]=c(1,0,1,0,1,1,1)
Q[30,]=c(1,1,0,1,1,0,1)

Daniel Turek

unread,

Oct 8, 2023, 7:16:59 AM10/8/23

to Duoduo, nimble-users

No problem, Duoduo. Thanks for sending the Q matrix. First, I don't see any NA's in the response (data) value of y:

any(is.na(data$y))
## [1] FALSE

And the same is true of "vec1" and "vec2" also, so I'm not entirely sure what you mean by "only consider observed response in model and skip the NA index".

That said, there's an error or two in your model code. For one, you can't have declarations using an equal sign, like:

add=vec1[n]
index=vec2[add+i]

in your model code. For deterministic model declarations, you should use the R assignment operator, "<-".

But more important, the BUGS language is declarative, not imperative, meaning each declaration (and each iteration inside a for-loop) must declare a unique node. So this code:

for (n in 1:1000) {
for (i in 1:vec1[n+1]) {
add=vec1[n]
index=vec2[add+i]

.....
}}

Results in multiple declarations of "add" and "index", which you can't have. But, we can easily take care of this by defining the equivalent constants, as:

index <- array(NA, c(1000, max(data$vec1)))

for(n in 1:1000) {
for(i in 1:max(data$vec1)) {
index[n,i] <- data$vec2[data$vec1[n] + i]
}
}

(granted, this is not the most efficient way of doing this operation in R, but it's the most clear).

And then everywhere that "index" is used in these nested loops in your model code, replace it with "index[n,i]"

Also, you'll need to move "vec1" and "Q" to the constants list, and add "index" to the constants list. Some updated code is below:

#model we fit
code=nimbleCode({
for (n in 1:1000) {
for (i in 1:vec1[n+1]) {

for (k in 1:7) {

w[n, index[n,i], k] <- pow(attribute[n, k], Q[index[n,i], k])
}
logit_prob[1+(i-1)+vec1[n]] <- beta[index[n,i]] + delta[index[n,i]] * prod(w[n,index[n,i], 1:7])
prob[1+(i-1)+vec1[n]] <- exp(logit_prob[1+(i-1)+vec1[n]])

y[1+(i-1)+vec1[n]] ~ dbern(prob[1+(i-1)+vec1[n]])
}}
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)

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

constants <- list(index = index, Q=Q,vec1=vec1)

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

I think this will get you running, and hopefully you can move forward from here.

Cheers,
Daniel

To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/e262efeb-50f4-4def-859f-c39c3d20cfb9n%40googlegroups.com.

Duoduo

unread,

Oct 10, 2023, 12:31:47 PM10/10/23

to nimble-users

Hi, Dan. Thanks for the detailed response! I am wondering if you could tell me what is the "data" from " any(is.na(data$y)) "? I have the NA of response y in the file I upload, then I try to find the location of NA. When I use y in the model, at the first time, I think maybe I could consider no NA (then I only have jumping index not complete 1000*30 index). I know that if I directly include NA in response, our package can help us set the NA part as a latent variable, I am trying to see if I can only include the observed response in the model and get the posterior estimate of parameter first. THanks!

Daniel Turek

unread,

Oct 13, 2023, 9:38:22 AM10/13/23

to Duoduo, nimble-users

Duoduo, so is it that you want to write your model, so that only the non-NA y values are used for fitting the parameters? (and the NA values are not included in the model at all?)

If so, then you can accomplish that using a little clever indexing, which might be overall simpler than what it looks like you've done in your code with vec1, loc, try, trytry, etc.

I sketched out the main idea, but the details of the arguments to dbern(...) need to be finished:

y_NA=fread("missing_response.csv")
numNonNAs <- apply(y_NA, 1, function(x) sum(!is.na(x)))
nonNAindices <- array(NA, dim(y_NA))
for(i in 1:nrow(y_NA)) nonNAindices[i,1:numNonNAs[i]] <- which(!is.na(y_NA[i,]))

## include numNonNAs and nonNAindices in your constants list:
constants$numNonNAs <- numNonNAs
constants$nonNAindices <- nonNAindices

## include y_NA in your data list:
data$y_NA <- y_NA

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

##for (i in 1:vec1[n+1]) {
## add=vec1[n]
## index=vec2[add+i]
## for (k in 1:7) {w[n, index, k] <- pow(attribute[n, k], Q[index, k])}
## logit(prob[1+(i-1)+vec1[n]]) <- beta[index] + delta[index] * prod(w[n,index, 1:7])
## y[1+(i-1)+vec1[n]] ~ dbern(prob[1+(i-1)+vec1[n]])
##}}
for(i in 1:numNonNAs[n]) {
y_NA[n, nonNAindices[n,i]] ~ dbern(....) ## 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)
})

I hope this helps get you started.

Cheers,

Daniel

To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/226a58a0-3ea3-4ac3-b2e8-ffc3de2f139en%40googlegroups.com.

Duoduo

unread,

Oct 17, 2023, 6:28:52 PM10/17/23

to nimble-users

Hi, Dan. This idea to ignore the missing part index is excellent! It is exactly what I am struggling with. Thanks for the help!

Duoduo

unread,

Oct 23, 2023, 1:10:13 PM10/23/23

to nimble-users

Hi, Dan. I tried the method you mentioned to run the model. 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. Do you have any idea about this? 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,