Modelling latent correlated variables in repeated-measures (aka. "within-subjects") design

[Mike Lawrence]

Consider data-generating code as follows:

transformed data{

  vector[2] coef_means ; #intercept & effect means

  vector[2] coef_sds ; #intercept & effect sds

  corr_matrix[2] cor_mat ; #correlation btwn intercept & effect

  matrix[2,2] contrasts ; #contrasts to generate observed data

  #filling correlation matrix

  cor_mat[1,1] = 1 ;

  cor_mat[1,2] = .6 ;

  cor_mat[2,1] = .6 ;

  cor_mat[2,2] = 1 ;

  #fill means & sds

  coef_means[1] = 0 ; #intercept mean

  coef_means[2] = 0 ; #effect mean

  coef_sds[1] = 1 ; #intercept sd

  coef_sds[2] = 1 ; #effect sd

  #fill contrasts

  contrasts[1,1] = 1 ; #intercept

  contrasts[2,1] = 1 ; #intercept

  contrasts[1,2] = -1 ; #effect * -1

  contrasts[2,2] = 1 ; #effect * 1

}

parameters{

  vector[2] coefs ; #observed intercept & effect

}

model{

  #observed intercept & effect as multivariate normal

  coefs ~ multi_normal(

      coef_means

    , quad_form_diag( cor_mat , coef_sds )

  ) ;

}

generated quantities{

  matrix[2,3] dat ;

  dat[,1:2] = contrasts ; #save the contrasts used

  dat[,3] = contrasts * coefs ; #compute observed vals

}

saved as "make_data.stan" and used to generate data via R:

library(rstan)

nSubj = 100
#generate nSubj "subjects" of observed data

from_stan = rstan::stan(

  file = 'make_data.stan'

  , iter = nSubj*2

  , seed = 1

  , chains = 1

)

dat = rstan::extract(from_stan)$dat

dat = matrix(dat,ncol=3) #unwrap array

X = cbind(dat[,1],dat[,2]) #contrasts

Y = dat[,3] #observed values

subjID = rep(1:nSubj,times=2) #subject labels

Now, given this data (contrasts X and values Y for each of subject) I'd like to recover the original data-generating parameters (coef_means, coef_sds, cor), but I'm unsure how to structure the model. 

I know I can compute estimates of the coefs for each subject outside stan via lm:

est_coefs = matrix(NA,ncol=2,nrow=nSubj)

for(i in 1:nSubj){

  est_coefs[i,] = coef(lm(Y[subjID==i]~0+X[subjID==i,]))

}

data_for_stan = list(

  nSubj = nSubj

  , est_coefs = est_coefs

)

then run stan model like: 

data {

  int<lower=1> nSubj ;

  matrix[nSubj,2] est_coefs ;

}

parameters {

  vector[2] coef_means ;

  vector<lower=0>[2] coef_sds ;

  corr_matrix[2] cor_mat ;

}

model {

  coef_means ~ normal(0,1) ;

  coef_sds ~ weibull(2,1) ;

  cor_mat ~ lkj_corr(2) ;

  for(s in 1:nSubj){

    est_coefs[s,] ~ multi_normal( coef_means , quad_form_diag(cor_mat, coef_sds) ) ;

  }

}

But I'm wondering if there's a way to pass the X & Y & subjID data objects in directly rather than computing the coefficient estimates by hand in R first, not only because I'm lazy, but also because I'm eventually going to have more complicated coefficient designs and missing data at the observation level. 

Thoughts?

[Mike Lawrence]

Managed to work it out!

Given the data generated as previously described, one can compute the coefficients for each subject inside stan via:

data {

  int<lower=1> nSubj ;

  int<lower=1> nY ;

  int<lower=1> nX ;

  vector[nY] Y ;

  matrix[nY,nX] X ;

  int subjIndices[nSubj,2] ; #indices of start & finish of each subject in Y

}

transformed data{

  matrix[nSubj,nX] obs_coefs ;

  for(s in 1:nSubj){

    obs_coefs[s] = transpose(

      inverse(

        transpose( X[subjIndices[s,1]:subjIndices[s,2]] )

        * X[subjIndices[s,1]:subjIndices[s,2]]

      )

      * transpose( X[subjIndices[s,1]:subjIndices[s,2]] )

      * Y[subjIndices[s,1]:subjIndices[s,2]]

    ) ;

  }

}

So here would be the full model, with some scaling and unscaling of the coefficients for easy specification of weakly informed priors, but lacking any other optimizations:

data {

  int<lower=1> nSubj ;

  int<lower=1> nY ;

  int<lower=1> nX ;

  vector[nY] Y ;

  matrix[nY,nX] X ;

  int subjIndices[nSubj,2] ;

}

transformed data{

  matrix[nSubj,nX] obs_coefs ;

  matrix[nSubj,nX] scaled_coefs ;

  for(s in 1:nSubj){

    obs_coefs[s] = transpose(

      inverse(

        transpose( X[subjIndices[s,1]:subjIndices[s,2]] )

        * X[subjIndices[s,1]:subjIndices[s,2]]

      )

      * transpose( X[subjIndices[s,1]:subjIndices[s,2]] )

      * Y[subjIndices[s,1]:subjIndices[s,2]]

    );

  }

  #scale the coefs for easy default priors

  for(i in 1:nX){

    scaled_coefs[,i] = (obs_coefs[,i] - mean(obs_coefs[,i])) / sd(obs_coefs[,i]) ;

  }

}

parameters {

  vector[nX] scaled_coef_means ;

  vector<lower=0>[nX] scaled_coef_sds ;

  corr_matrix[nX] cor_mat ;

}

model {

  scaled_coef_means ~ normal(0,1) ;

  scaled_coef_sds ~ weibull(2,1) ;

  cor_mat ~ lkj_corr(1) ;

  for(s in 1:nSubj){

    scaled_coefs[s,] ~ multi_normal( scaled_coef_means , quad_form_diag(cor_mat, scaled_coef_sds) ) ;

  }

}

generated quantities{

  vector[nX] coef_means ;

  for(i in 1:nX){

    coef_means[i] = scaled_coef_means[i]

      * sd(obs_coefs[,i])

      + mean(obs_coefs[,i])

      ;

  }

}

[Mike Lawrence]

FYI, I've posted vids/data/code for this model:

lecture 1: https://www.youtube.com/watch?v=J15JHSeLbdg&index=17&list=PLu77iLvsj_GNmWqDdX-26kZ56dCZqkzBO

lecture 2: https://www.youtube.com/watch?v=tmFAXOuhloA&index=18&list=PLu77iLvsj_GNmWqDdX-26kZ56dCZqkzBO

Where I do more introductory stuff in the first and talk about the code from this thread (inc. application to missing data) in the second.

--
Mike Lawrence
Graduate Student
Department of Psychology & Neuroscience
Dalhousie University

~ Certainty is (possibly) folly ~

--
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+unsubscribe@googlegroups.com.
To post to this group, send email to stan-...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Bob Carpenter]

Thanks for sharing. Are you using all these lectures for
a class? I saw them tagged as PSYO6001.

- Bob

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

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

[Mike Lawrence]

Yup, it's a stats class for first year grad students in psychology & neuroscience.