Missing Data: multivariate normal, how to make it work?

[Markian Pahuta]

Hi,

I have written a hierarchical model with a MVN likelihood. I am able to get the model to run when no data is missing (use file full.dat). But when I run the model with missing data (missing.dat), I get the error:

Error : Error in function stan::prob::multi_normal_log(N4stan5agrad3varE): Random variable[0]  is 1.#QNAN:0, but must not be nan!

error occurred during calling the sampler; sampling not done

 The RSTAN model is: 

library(rstan)

mcode <- '

data{

int <lower =1> Nrep[2];

int <lower =1> Nsubjects[2];

int <lower=1> nRowMan;

int <lower=1> man_ind[nRowMan,3];

real <lower=0> manlongx[nRowMan];

int <lower=1> nRowGrav;

int <lower=1> grav_ind[nRowGrav, 3];

real <lower=0> gravlongx[nRowGrav];

vector[2] basemu;

}

transformed data{

vector [2]  mu_ln2;

vector [2]  sigma_ln2;

real eta;

mu_ln2[1] <- 0.0;

mu_ln2[2] <- 0.0;

sigma_ln2[1] <- 0.5;

sigma_ln2[2] <- 0.5;

eta <- 1.0;

}

parameters{

vector <lower=0>[2] muoverall;

vector <lower=0>[2] mustudy[2];

vector <lower=0>[2] musubject[2, 30];

cov_matrix[2] A;

cov_matrix[2] B;

cov_matrix[2] C [2];

cov_matrix[2] D [2];

}

transformed parameters{

vector[2] xobserved[2,30,3];

for (a in 1:nRowMan) xobserved[man_ind[a,1],man_ind[a,2],man_ind[a,3],1] <- manlongx[a];

for (b in 1:nRowGrav) xobserved[grav_ind[b,1],grav_ind[b,2],grav_ind[b,3],2] <- gravlongx[b];

}

model {

A ~ lkj_cov(mu_ln2, sigma_ln2, eta);

B ~ lkj_cov(mu_ln2, sigma_ln2, eta);

muoverall~multi_normal(basemu, A);

for (study in 1:2){

mustudy[study] ~ multi_normal (muoverall, B);

C[study] ~ lkj_cov(mu_ln2, sigma_ln2, eta);

D[study] ~ lkj_cov(mu_ln2, sigma_ln2, eta);

for (subject in 1:Nsubjects[study]){

musubject[study, subject] ~ multi_normal(mustudy[study], C[study]);

xobserved[study, subject,1] ~ multi_normal(musubject[study, subject], D[study]);

}

}

}

generated quantities{

real bias;

bias<-muoverall[2]-muoverall[1];

}

'




load('full.dat')

#load('missing.dat')




datalist = list(

Nsubjects=c(25,26),

nRowMan =nrow(man_ind),

man_ind =man_ind,

manlongx =manlongx,

nRowGrav =nrow(grav_ind),

grav_ind =grav_ind,

gravlongx =gravlongx,

basemu=c(4,4),

Nrep=c(1,3))




modfit <- stan(model_code = mcode, data = datalist,iter = 1000, chains = 4)

I was hoping someone could please help me figure this out?

Thank you

Mark

[Bob Carpenter]

Unlike BUGS, Stan does not have a way to represent R's NA (undefined/missing)
values. So you can't provide them as part of input in RStan and therefore
can't directly translate BUGS-style missing-data models. There's a chapter
of the manual on coding missing-data models in Stan. It requires preprocessing
the data into a different format and also fiddling with the model.

- Bob

On 8/10/13 4:38 PM, Markian Pahuta wrote:
> Hi,
>
> I have written a hierarchical model with a MVN likelihood. I am able to get the model to run when no data is missing
> (use file full.dat). But when I run the model with missing data (missing.dat), I get the error:
>

> /Error : Error in function stan::prob::multi_normal_log(N4stan5agrad3varE): Random variable[0] is 1.#QNAN:0,
> but must not be nan!/
>
> /error occurred during calling the sampler; sampling not done/

> --
> You received this message because you are subscribed to the Google Groups "stan users mailing list" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.
> For more options, visit https://groups.google.com/groups/opt_out.
>
>

[Markian Pahuta]

Hi Bob,

Thank you for your reply. I am aware of this difference, and used the approach in section 6.4 of the manual. My input data format has no NAs.

Does the error mean that I haven't implemented the transformed parameters correctly?

[Bob Carpenter]

Sorry -- I should've looked at the model more closely.

The error you're seeing arises when a local or transformed
variable is not properly initialized or when data is initialized
as NaN. Parameters are always initialized, so you don't need to
worry about them. So the problem has to be xobserved, because it's
the only non-parameter implicated:

> xobserved[study, subject,1] ~ multi_normal(musubject[study, subject], D[study]);

It must be xobserved because everything else is a parameter.
This definition is suspicious because it might not initialize
all of xobserved depending on the values of grav_ind[b].

> for (b in 1:nRowGrav) xobserved[grav_ind[b,1],grav_ind[b,2],grav_ind[b,3],2] <- gravlongx[b];

You can try printing it before you use it to debug where this
is happening.

- Bob

> stan-users+...@googlegroups.com <javascript:>.
> > For more options, visit https://groups.google.com/groups/opt_out <https://groups.google.com/groups/opt_out>.

[Markian Pahuta]

Hi Bob,

Thank you for suggesting the print function - I hadn't thought of that!

So my code is correctly transforming parameters/data, and the indexes are all working correctly in the likelihood loop. In my data vector for the MVN, occasionally one of the observations will be missing:

[1.#QNAN:0,3.5:0]

and sometimes both will be missing

[1.#QNAN:0,1.#QNAN:0]

Do I have to do further processing to get the MVN to work? It seems that the values of my missing observations "#QNAN:0" reflect the error I am getting when I run the model?

Thanks

Mark

[Bob Carpenter]

On 8/12/13 10:00 PM, Markian Pahuta wrote:
> Hi Bob,
>

> Thank you for suggesting the print function - I hadn't thought of that!
>
> So my code is correctly transforming parameters/data, and the indexes are all working correctly in the likelihood loop.
> In my data vector for the MVN, occasionally one of the observations will be missing:
> [1.#QNAN:0,3.5:0]
> and sometimes both will be missing
> [1.#QNAN:0,1.#QNAN:0]
>
> Do I have to do further processing to get the MVN to work? It seems that the values of my missing observations "#QNAN:0"
> reflect the error I am getting when I run the model?

Yes, you can't put NaN into multivariate normal because
it leads to a NaN result in the density.

- Bob

> <https://groups.google.com/groups/opt_out> <https://groups.google.com/groups/opt_out

[Markian Pahuta]

Hi,

 

I thought STAN allowed missing values in the MVN? How do I code these values then?

[Bob Carpenter]

I don't know how to make it work for multi-variate normal, but
others might know.

You can read the chapter on missing data in the manual, but I
don't have any advice other than that for your particular case.

- Bob

[Markian Pahuta]

Hi,

I figured out my problem: missing observations must be declared as a parameter variable, and then assigned to the correct indices of transformed parameter variable in place of observed data. Everything works well now.

[Kuk-Hyun Ahn]

Dear Markian Pahuta

I have same issue with you.

If possible, can you provide a sample code for this problem?

I do now know how to deal with missing value in multivariate normal distributions


2013년 8월 17일 토요일 오후 3시 13분 48초 UTC-4, Markian Pahuta 님의 말:

[Bob Carpenter]

Is the issue that you have something like

y ~ multi_normal(...);

where y has missing data? It's straightforward, but involves
some index fiddling --- there's a chapter in the manual
that explains how to deal with missing data in Stan. The basic idea
is that you have to declare the missing data as raw parameters,
code up the y components that exist as data, then put the two together
into the y used in the above statement.

- Bob

> You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.

> To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.

> To post to this group, send email to stan-...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.

[Kuk-Hyun Ahn]

Yes, only Y has missinh values.
Could you suggest me any example for this case?
I have a trouble to understand it

[Bob Carpenter]

It's the very first example in chapter 8 of the manual,
"Missing Data and Partially Known Parameters". You can
use exactly the same approach, only you'll have a matrix
or array of vectors rather than a single array.

- Bob