Reproducing a structural equation model and conditional density in RStan

[Robert Ness]

Hello,

I am new to Stan and learning the ropes.

I have two variables X and Y

It turns out that there is a structural relationship where X determines Y: E(Y|X) = M* wX / (1 + wX).  M is a constant, w is the parameter that I want to do inference on.

This is not evident in the plot because X itself is random in the data and values concentrate around the modes.  But you can see the relationship a bit if you condition on X.  The following graph was generated by conditional kernel density estimation;

So if I fit the curve to the mode of Y for each value of X, I get nice estimates for w.  The conditional kernel density approach I used is nonparametric.  I want a parametric approach, namely the Gamma distribution: Y ~ Gamma(shape, shape / E(Y|X)).  I want to model Y directly, not the modes of Y|X spit out by some blackbox function. 

I tried the following Stan approach 

library(rstan)

options(mc.cores = parallel::detectCores())

stan_code <- "

data {

  int<lower=0> N;

  real<lower=0> M;

  vector[N] x;

  vector[N] y;

}

parameters {

  real <lower=0.001,upper=10> w;

  real <lower=0.001,upper=10> shape;

}

model {

  for(i in 1:N)

    y[i] ~ gamma(shape, shape / (M * (  w * x[i]  / (1 + w * x[i]))));

  }

"

bayes_model <- stan(model_code = stan_code,

    data = list(x = .data$X,

                y = .data$Y,

                N = nrow(.data),

                M = M))  

 The code runs but it is not working out.  I want the approach to take into account that X is a random variable (I assume Gamma(a, b)), and that less probable observations of X are just as important in inference on w as observations of X around the modes.  

More specifically, the result is giving me a posterior distribution for w with a mean not to far off from the output of glm().  I need to reflect the fact that X is a random variable as well, and somehow add "weight" to observations where the marginal density of X is low, just like conditional kernel density function f(y|x) "weights" the joint kernel density f(x, y) with the inverse of the marginal kernel density of x,  1/f(x).  

Any ideas?

[Bob Carpenter]

First, we tend to think like Bayesians and treat X as a random variable
all the time. The question I think you have is whether it's measured
with noise or not. You can check out the manual chapter on measurement
error models.

You might want to check out our manual chapter on priors, too. You
won't need that lower=0.001 business in Stan and we don't recommend
hard upper-bound constraints on positive-constrained parameters as they
don't do what people expect them to do. I need to get a case study on
this written because I'm writing it too frequently in the mailing list
(second time in the last ten messages or so), because people copy
the practices of BUGS.

I don't quite understand what you're asking w.r.t. kernel density
estimation, but maybe somebody else will.

- Bob

> On Aug 13, 2016, at 3:15 AM, Robert Ness <rober...@gmail.com> wrote:
>
> Hello,
>
> I am new to Stan and learning the ropes.
>
> I have two variables X and Y
>
>
>
>

> It turns out that there is a structural relationship where X determines Y: E(Y|X) = M* wX / (1 + wX). M is a constant, w is the parameter that I want to do inference on.
>
> This is not evident in the plot because X itself is random in the data and values concentrate around the modes. But you can see the relationship a bit if you condition on X. The following graph was generated by conditional kernel density estimation;
>
>
>
>
>
>

> So if I fit the curve to the mode of Y for each value of X, I get nice estimates for w. The conditional kernel density approach I used is nonparametric. I want a parametric approach, namely the Gamma distribution: Y ~ Gamma(shape, shape / E(Y|X)). I want to model Y directly, not the modes of Y|X spit out by some blackbox function.
>
> I tried the following Stan approach
>
>
> library(rstan)
> options(mc.cores = parallel::detectCores())
> stan_code <- "
> data {
> int<lower=0> N;
> real<lower=0> M;
> vector[N] x;
> vector[N] y;
> }
> parameters {
> real <lower=0.001,upper=10> w;
> real <lower=0.001,upper=10> shape;
> }
> model {
> for(i in 1:N)
> y[i] ~ gamma(shape, shape / (M * ( w * x[i] / (1 + w * x[i]))));
> }
> "
> bayes_model <- stan(model_code = stan_code,
> data = list(x = .data$X,
> y = .data$Y,
> N = nrow(.data),
> M = M))
>
> The code runs but it is not working out. I want the approach to take into account that X is a random variable (I assume Gamma(a, b)), and that less probable observations of X are just as important in inference on w as observations of X around the modes.
>
> More specifically, the result is giving me a posterior distribution for w with a mean not to far off from the output of glm(). I need to reflect the fact that X is a random variable as well, and somehow add "weight" to observations where the marginal density of X is low, just like conditional kernel density function f(y|x) "weights" the joint kernel density f(x, y) with the inverse of the marginal kernel density of x, 1/f(x).
>
> Any ideas?
>

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