partially missing data in bivariate normal distribution

[Jaslene Lin]

Dear Stan users:

I am new to using Rstan and I would like to implement a previous working model coded in R2WinBUGS to Rstan.

My data is bivariate and normally distributed. yi=(yi1,yi2) for i=1...n and I would like to fit a multivariate normal model. However there are some missing data on yi2 for some observations but their yi1 are all present. I read chapter 8 and searched for relevant topics within the group but I do not seem to find a solution to import my data set into Rstan as my some of my observations are partially missing. Are there any examples or documents that I could read up to solve my very trivial problem?

Thanks in advance,

Jaslene

[Daniel Lee]

Have you looked at Ch. 8 in the v2.9.0 version of the manual titled "Missing Data & Partially Known Parameters"?

Daniel

--
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.

[Jaslene Lin]

Hi Daniel:

Thank you for your quick reply. Yes I have. In that example, if I understand correctly, y is only univariate, so we separately modelling observed and missing cases.

Mine is bivariate case with a snap shot of data looking like this,

> cr
          [,1]         [,2]
 [1,] 3.149149 -0.140450679
 [2,] 4.431096  0.324657435
 [3,] 3.534658           NA
 [4,] 4.322317  2.515527025
 [5,] 4.579680  0.135935029
 [6,] 3.829151  1.442398431
 [7,] 3.340404 -0.551143435
 [8,] 4.107840  2.255137882
 [9,] 2.125653 -0.007589424
[10,] 3.469233 -0.068156790
[11,] 3.906750  0.734196341
[12,] 3.854476  0.804081973
[13,] 3.541880           NA
[14,] 3.903313  2.529434036
[15,] 3.829010  0.904271246
[16,] 3.885068           NA
[17,] 2.962250  0.358900085
[18,] 4.340386  1.461094418
[19,] 2.462504 -0.588568428
[20,] 6.913772  3.848423718
[21,] 2.954179  0.231237494
[22,] 5.089764           NA

what you have suggested is that I can model all complete bivariate observations and then fit a model to observations where I partially observe?

Thanks in advance,

Jaslene

[Daniel Lee]

Yes, you will want to separate out observed from missing and pass those in as separate data structures.

[Jaslene Lin]

Thanks, Daniel.

I will give it a try now.

Cheers,

Jaslene

[Jaslene Lin]

I tried separating my dataset but it gives me the following error message

> cr_miss
         [,1] [,2]
[1,] 3.534658   NA
[2,] 3.541880   NA
[3,] 3.885068   NA
[4,] 5.089764   NA
[5,] 3.051150   NA
> hierarchical_dat1<-list(N=45,M=5,D=2,cr_obs=cr_obs,cr_miss=cr_miss)
> 
> fit1<-stan(file='coding missing data.stan',data=hierarchical_dat1,iter=1000,chains=4)
Error in FUN(X[[i]], ...) : Stan does not support NA (in cr_miss) in data
failed to preprocess the data; sampling not done

 

data {
 
  int<lower=0> N;// num observations with complete midpoints and log(range)
 
  int<lower=0> M;// num observations with missing values on log(range)
 
  int<lower=0> D;// dimension of observations
 
  vector[D] cr_obs[N];//observations with complete midpoints and log(range)
 
  vector[D] cr_miss[M];//observations with missing values on log(range)
 
 
 
 
 }

parameters {
  vector[D] mu; //prior location for (midppoints and log(range))
 
  cholesky_factor_corr[D] L_Omega;// cholesky factor for correlation matrix Omega
 
  vector<lower=0>[D] sigma;//prior scale

}

transformed parameters{
 
  matrix[D,D] L_Sigma;
   
  L_Sigma<-diag_pre_multiply(sigma, L_Omega);// cholesky factor for variance and covariance matrix
   
  //matrix[D,D] Sigma;// variance and covariance matirx
 
 // Sigma<-diag_pre_multiply(sigma, L_Omega)*diag_pre_multiply(sigma,L_Omega)';
 
  
}

 
model {
   
    mu[1]~uniform(0,100);
   
    mu[2]~normal(0,100);
   
    sigma~cauchy(0,2.5);
   
    L_Omega~lkj_corr_cholesky(2);
   
   
    cr_obs~multi_normal_cholesky(mu,L_Sigma);
   
    cr_miss~multi_normal_cholesky(mu,L_Sigma);
                       
    }

 

In this case, I think I need to recode cr_miss as what the manual suggests in chapter 13 sparse and ragged data structures?

Thanks in advance,

Jaslene

[Daniel Lee]

You're providing the first column of cr_miss, so you should make the data you pass in a vector.

In parameters, you're estimating the second column of cr_miss, so that should be a parameter vector.

In transformed parameters, you'll want to declare cr_miss as a matrix with the first column being data and the second being parameters.

There might be a more clever way to do it, but it's late and I'm not so familiar with what you're trying to do. (I didn't actually try to read your Stan program.)

Daniel

--

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

[Jaslene Lin]

Thank you so much for your help(especially at this time, it's afternoon in Sydney) , I am understanding it gradually now.

Cheers,

Jaslene

[Michael Betancourt]

Jaslene,

If everything is bivariate Gaussian then you can do everything analytically.

First, split your data into non-missing (which will be two-dimensional)

and missing the second component (which will be one-dimensional).

The non-missing data will be distributed according to the bivariate

Gaussian that you quoted

(y_1, y_2) ~ multi_normal( (mu_1, mu_2), Sigma ),

where Sigma is the 2x2 covariance matrix.

In order to construct the distribution of the missing data we have to

marginalize y_{2} out of the bivariate Gaussian.  Fortunately, this

is straightforward,

y_1 ~ normal(mu_1, sigma_1),

where sigma_1 is the square root of Sigma_{11}.

Your Stan model would then look _something_ like

data {

  int N_nonmiss;

  vector[2] y_nonmiss[N_nonmiss];

  int N_miss;

  real y_miss[N_miss];

}

parameters {

  vector[2] mu;

  cov_matrix[2] Sigma;

}

model {

  // Prior on mu

  // Prior on Sigma, see the manual for examples

  y_nonmiss ~ multi_normal(mu, Sigma);

  y_miss ~ normal(mu[1], sqrt(Sigma[1, 1]);

[Bob Carpenter]

Cool. I'm going to add that example to the manual:

https://github.com/stan-dev/stan/issues/1238#issuecomment-224378060

I'd rename the variable Michael called y_miss because
it's actually the observed y[n,1], not the missing
y[n,2].

The only drawback is that it doesn't impute the missing
data for you. If you need to do that, there's probably
a way to analytically derive the conditional

p(y[n,2] | y[n,1], mu, Sigma)

Then you could add that to the generated quantities block
for maximum efficiency (both in per iteration time and
N_eff per iteration). Or declare the y[n,2] as parameters if
you need to use them elsewhere.

- Bob

[Michael Betancourt]

https://en.wikipedia.org/wiki/Multivariate_normal_distribution#Conditional_distributions

[Michael Betancourt]

This is a nice example that emphasizes that missing
data is a _modeling_ problem, not something that
can be automated.

On Jun 7, 2016, at 7:52 PM, Bob Carpenter <ca...@alias-i.com> wrote:

[Bob Carpenter]

Indeed.

The default BUGS/JAGS approach is probably not the model of
missingness most people are after. I had a nice discussion of
this with Merlise Clyde of Duke in the last plus one MCMCski --- she
was telling me about the lengths she went to coding up models for missing data
in BUGS/JAGS (but I don't have a reference) conditioned on other data.
When you do that and define a proper joint distribution, the code
looks much more like Stan code than the direct approach of just plugging
NA into a full-data model.

- Bob

P.S. Thanks for the reference to the conditional y[n,2] | y[n,1] form.

[Jaslene Lin]

Hi All:

Thank you so much for your quick replies. I have had a busy morning (OZ time) and have not had a chance to review all your comments. I will read them thoroughly and give them a go now.

Yesterday, I tried what Daniel has suggested but I think there still has some errors in my construction and did not work.

Kindest regards,

Jaslene

[Jaslene Lin]

Hi All:

I tried what Daniel had suggested and I separated my data into complete N observations  and M observations with partially observed y[n,2]

However, I seemed to run into the trouble of reconstructing the matrix of cr_miss

Here is my code

data {
 
  int<lower=0> N;// num observations with complete midpoints and log(range)
 
  int<lower=0> M;// num observations with missing values on log(range)
 
  int<lower=0> D;// dimension of observations
 
  vector[D] cr_obs[N];//observations with complete midpoints and log(range)
 

  veotor[1] cr_miss_midpoints[M];//observations with complete midpoints but missing log(range)
 
 
 
 
 
 }

parameters {
  vector[D] mu; //prior location for (midppoints and log(range))
 
  cholesky_factor_corr[D] L_Omega;// cholesky factor for correlation matrix Omega
 
  vector<lower=0>[D] sigma;//prior scale
 

  vector[1] cr_miss_r[M];// observations with missing log(range);

}

transformed parameters{
 
  matrix[D,D] L_Sigma;
 

  matrix [M,D] cr_miss; matrix for containing observed data (midpoints) and unobserved log(range) now treated as parameters to be estimated
  
 
 
  
   for(m in 1:M)
  cr_miss[m,D-1]<-cr_miss_midpoints;
 
 
 
  for(m in 1:M)
  cr_miss[m,D]<-cr_miss_r;

 
   
  L_Sigma<-diag_pre_multiply(sigma, L_Omega);// cholesky factor for variance and covariance matrix
   
  //matrix[D,D] Sigma;// variance and covariance matirx
 
 // Sigma<-diag_pre_multiply(sigma, L_Omega)*diag_pre_multiply(sigma,L_Omega)';
 

 
 
} 

I highlighted the part I would like some advice on/

I realised that STAN is row-major. How can I assign columns of a matrix to a column vector?

Thanks in advance,

Jaslene

 

 

[Jaslene Lin]

I have tried and applied conditional normal distribution for partially observed data and here is my code

data {
 
  int<lower=0> N;// num observations with complete midpoints and log(range)
 
  int<lower=0> M;// num observations with missing values on log(range)
 
  int<lower=0> D;// dimension of observations
 

  vector[D] cr_obs[N];//two-dimensional observations with complete midpoints and log(range)
 
  real cr_miss[M];//oen-dimensional observations with complete midpoints but missing log(range)

 
 
  }
parameters {
  vector[D] mu; //prior location for (midppoints and log(range))
 
  cholesky_factor_corr[D] L_Omega;// cholesky factor for correlation matrix Omega
 

  //matrix[D,D] Sigma;//variance and covariance matrix

 
  vector<lower=0>[D] sigma;//prior scale
 
 }

transformed parameters{
  matrix[D,D] Sigma;//variance and covariance matrix
 
  matrix[D,D] L_Sigma;
 
  Sigma<-diag_pre_multiply(sigma,L_Omega)*diag_pre_multiply(sigma,L_Omega)';

 
  L_Sigma<-diag_pre_multiply(sigma, L_Omega);// cholesky factor for variance and covariance matrix
}

model {
   
    mu[1]~uniform(0,100);
   
    mu[2]~normal(0,100);
   
    sigma~cauchy(0,2.5);
   
    L_Omega~lkj_corr_cholesky(2);
   
   
    cr_obs~multi_normal_cholesky(mu,L_Sigma);
   

    cr_miss~normal(mu[1],sqrt(Sigma[1,1]));
                       
    }

 

However, here are some messages I received when running the model I am not sure whether these errors are major and implying there is something wrong with the construction of my model or simply just warning messages.

[1] "The following numerical problems occured the indicated number of times after warmup on chain 4"
                                                                                count
Exception thrown at line 49: normal_log: Scale parameter is 0, but must be > 0!     1
[1] "When a numerical problem occurs, the Metropolis proposal gets rejected."
[1] "However, by design Metropolis proposals sometimes get rejected  even when there are no numerical problems."
[1] "Thus, if the number in the 'count' column is small,  do not ask about this message on stan-users."

Could you please offer me some advice in terms of dealing with error messages? or they implying something fatally wrong about my model?

Kindest regards,

Jaslene

 

[Bob Carpenter]

That's just a warning that there was a numerical issue that
was trapped. As the message says, if there aren't very many
of these, you don't have to worry abot it.

To get rid of them, you can try increasing adapt_delta with

control=list(adapt_delta = 0.99)

and see if that helps. These issues arise when numerical
stability problems creep in under the hood with the Hamiltonian
simulation coupled to the unconstraining transforms plus Jacobians
that we apply automatically.

- Bob

[Jaslene Lin]

Thanks Bob. I still would like to know how can I assign a column vector to a column of a matrix. I have been reading the manual but does not seem to find the solution.

I know Rstan is row-major index, so if I specify matrix [M,N] Sigma;

then for (m in 1:M)

Sigma[m]<-assigning some values to the row of the matrix Sigma;

However, can I select column vector of Sigma matrix so that I can assign some values to it?

Thanks for your help.

Kindest regards,

Jaslene

[Bob Carpenter]

matrix[M, N] a;
vector[M] b;

// assign b to k-th column of a
a[, k] <- b;

[Jaslene Lin]

Thanks Bob and all others for helping me. my model finally works. I will continue reading the user manual.

Kindest regards,

Jaslene