Simple "within-subjects" design

[Mike Lawrence]

I just know I'm going to smack my head on seeing the solution when posted, but I'm having trouble coding a simple "within-subjects" design where simulated data might be generated as follows in R:

#generate some fake data
N = 20

#parameters
intercept_mean = 100
intercept_sd = 15
effect_mean = 10
effect_sd = 3
intercept_effect_correlation = 0


#covariance matrix
v = matrix(
    data = c(
        intercept_sd^2
        , intercept_sd*effect_sd*intercept_effect_correlation
        , intercept_sd*effect_sd*intercept_effect_correlation
        , effect_sd^2
    )
    , nrow = 2
    , ncol = 2
)

obs = MASS::mvrnorm(N,c(intercept_mean,effect_mean),v)

#combine to long data format
long_data = data.frame(
    S = rep(1:N,times=2) #subject labels
    , Y = c(
        obs[,1]-obs[,2]/2
        , obs[,1]+obs[,2]/2
    )
    , condition = factor(rep(c('A','B'),each=N))
)

Now, if one were interested estimating credible intervals on the parameters listed in the "parameters" section, I'd first try a model that has no correlation between intercepts and effects, as:

data {
    int nY ; # number of observations
    int nS ; # number of subjects
    vector[nY] Y ; # observations
    vector[nY] X ; # condition (-.5,.5)
    int S[nY] ; # subject label (1:nS) for each observation
}
parameters {
    real intercept ; #intercept
    real effect ; #effect
    real<lower=0> intercept_sd ; #between-Ss sd on intercept
    real<lower=0> effect_sd ; #between-Ss sd on effect
    vector[nS] subject_intercept ; #storage vector for subject-by-subject intercept values
    vector[nS] subject_effect ; #storage vector for subject-by-subject effect values
}
transformed parameters{
    vector[nY] y ; #storage vector for observation-by-observation predicted y values
    for(ny in 1:nY){
        mu[ny] <- subject_intercept[S[ny]] + X[ny]*(subject_effect[S[ny]]/2) ;
    }
}
model {

    #priors for parameters
    intercept ~ normal(100,5) ; #prior on intercept
    effect ~ normal(0,5) ; #prior on effect
    intercept_sd ~ cauchy(0,15) ; #prior on intercept_sd
    effect_sd ~ cauchy(0,15) ; #prior on effect_sd

    #asserting that subjects are sampled from a population
    subject_intercept ~ normal(intercept,intercept_sd) ;
    subject_effect ~ normal(effect,effect_sd) ;

    for(ny in 1:nY){
        Y[ny] ~ y[ny] ;
    }
}

But obviously the line "Y[ny] ~ y[ny]" isn't well-formed and I need a way of expressing the model so that there's a sampling distribution on the right-hand side of that equation rather than a value computed from sampling distributions. Help?

[Mike Lawrence]

(Deleted my previous post as it provided an incorrect solution; working on the proper solution using a multivariate normal...)

[Mike Lawrence]

Nope, still flummoxed...

[Guido Biele]

assuming your Y is normally distributed, you need a sampling statement like

Y[ny] ~normal(y[ny],sigma)

you will have to declare sigma in the parameter block. one recommendation for the variance prior used to be the half cauchy (i.e. declare sigma with a lower bound 0 and sample from a chauchy distribution in the model block).

note that in the transformed parameters block you declare y, but you fill the undeclared vector mu. 

I don't think it is necessary to use the transformed parameters block here. all those mu/y will create a mess when you do print(stanfit).

you could just do the following in the model block

real mu; 

...

for(ny in 1:nY){

        mu <- subject_intercept[S[ny]] + X[ny]*(subject_effect[S[ny]]/2);
        Y[ny] ~ normal(mu,sigma) ;
    }

guido

PS: note that I did not look at your R code or tried to check if you calculated mu/y correctly! But I am wondering why you are dividing by two in that expression. for a simple within effect i would simply have an intercept plus a dummy variable (0,1) to capture the difference between the two conditions.

[Mike Lawrence]

On Tuesday, August 4, 2015, Guido Biele <guido...@gmail.com> wrote:

PS: note that I did not look at your R code 

Thanks for your time, but you would have saved said time had you done so. My data are not as you assumed so none of your suggestions apply.

--

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

~ Certainty is (possibly) folly ~

[Bob Carpenter]

I'm afraid I don't know what a "within-subjects" design is
(or if the quotes are meant to pick out a technical term or
indicate a non-standard usage), and the R code's a bit too complex
for me to understand without a serious time investment. There
are lots of us on the list who aren't R experts or even R users.

If you can write down the data, parameters, and log probability
function for the model, we can probably help you code it.

- 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.
> For more options, visit https://groups.google.com/d/optout.

[Mike Lawrence]

Resurrecting this older thread because I think I can clarify my statement of the problem now. 

Fundamentally this is an instance of multivariate regression, but one where you're not observing samples from the multivariate, but their linear transform in combination with condition contrasts. I believe the most direct analogy from the manual would be the example explored by "6.12. Multivariate Priors for Hierarchical Models" with the exception that only one "individual" is observed per "group", so you can't estimate the parameter sigma.

To express my data-generating R code in words and using a concrete example, we observe test scores of a class before and after lunch. We expect that students vary from one another in their overall performance and also vary from one another in the degree to which lunch/time affects their performance, and of course these two sources of variability could correlate. A given student's pre-lunch score is their intercept minus half their lunch effect and their post lunch score is their intercept plus half their lunch effect. Of course, I could simply compute each student's mean score and the difference between pre-lunch/post-lunch scores, and there's my data for the multivariate normal, but I'd have thought that there would be a way to express this model using just the observed test scores. Also, when I have this sorted out, I'm aiming to apply this to a more complicated design where the students are additionally observed on multiple days with some missing data for some students on some days, so direct computation of each subject's intercept & effects isn't possible. 

Does this clarification help?

[Bob Carpenter]

Yes, that makes sense (sorry for 3 month delay in answering; I'm
spring-cleaning my email). Andrew just suggested
exactly the same coding of pre/post-lunch scores in a different
application area to me yesterday.

- Bob

> To post to this group, send email to stan-...@googlegroups.com.

[Mike Lawrence]

Several months later, and I seem to have solved this finally! At least, for the scenario where you have data in wide format. I'll post the long-data compatible version in a bit, but in the meantime, here is the wide version that also permits missing data:

data {

  # N: number of cases

  int<lower=1> N ;

  # Y: matrix of outcomes for each case & variable

  matrix[N,2] Y ;

  # n_missing: number of instances of missing data

  int n_missing ;

  # indices_of_missing: matrix of the row and column indices of missing data

  int indices_of_missing[n_missing,2] ;

}

parameters {

  # means: vector of population-level means for intercept & slope

  vector[2] means ;

  # sds: vector of population-level sds for intercept & slope

  vector<lower=0>[2] sds ;

  # cor: population-level correlation between intercept & slope

  corr_matrix[2] cor ;

  # imputed: vector of imputed values for missing observations

  vector[n_missing] imputed ;

}

transformed parameters{

  # cov: covariance matrix (derived from cor & sds)

  matrix[2,2] cov ; #used in getting correlation to work

  # Y_with_imputed: matrix of observations including imputed values

  matrix[N,2] Y_with_imputed ;

  # Z: matrix of latent subject-by-subject intercepts & effects

  matrix[N,2] Z ;

  #replace any missing values in Y

  Y_with_imputed = Y ;

  for(i in 1:n_missing){

    Y_with_imputed[indices_of_missing[i][1],indices_of_missing[i][2]] = imputed[i] ;

  }

  #express Z as intercept/slope representation of Y_with_imputed

  for(i in 1:N){

    Z[i,1] = (Y_with_imputed[i,1]+Y_with_imputed[i,2])/2 ;

    Z[i,2] = Y_with_imputed[i,1]-Y_with_imputed[i,2] ;

  }

  #compute covariance matrix from cor & sds

  cov = quad_form_diag(cor,sds) ;

}

model {

  #priors on population parameters; presumes data has been scale()ed

  means ~ normal(0,1) ;

  sds ~ lognormal(0,1) ;

  cor ~ lkj_corr(1) ;

  #assert sampling of case-level outcomes on each variable given

  #  multivariate normal population

  for(s in 1:N){

    Z[s,] ~ multi_normal(means,cov) ;

  }

}