What's this model called? Regression with a PCA-rotated bivariate gaussian variable as outcome.

[matus.s...@privatdemail.net]

In my current modeling effort I came up with the following model. I wonder whether there is a literature that discusses a similar approach. Maybe, somebody here has already encountered something similar?

model='''
data{
    int<lower=0> N; //nr of subjects
    int<lower=0> S; //nr of sessions
    vector[2] y[N,S]; // bivariate gaussian outcome
    real x[S]; // predictor

}parameters{
    vector<lower=0>[2] sigma; // standard deviation
    corr_matrix[2] r;
    vector[2] offset;
    vector[2] slope;
}transformed parameters{
    vector[2] mu[S]; // mean for each session
    for (s in 1:S){
        // Here comes the crucial idea
        // 2d matrix multiplication in two lines follows
        // in matrix notation  mu[s]= sigma*U*sqrt(lambda)*(offset+slope*x[s])
        // where lambda has eigenvectors on diagonal and U is eigenvalue matrix of the correlation matrix r, ie r=U*lambda*U^T
        mu[s][1]<-sigma[1]*cos(pi()/4)*sqrt(1+r[1,2])*(offset[1]+slope[1]*x[s])-sigma[1]*sin(pi()/4)*sqrt(1-r[1,2])*(offset[2]+slope[2]*x[s]);
        mu[s][2]<-sigma[2]*sin(pi()/4)*sqrt(1+r[1,2])*(offset[1]+slope[1]*x[s])+sigma[2]*cos(pi()/4)*sqrt(1-r[1,2])*(offset[2]+slope[2]*x[s]);
    }
}model{
    matrix[2,2] tempSigma; // covariance matrix
    r~lkj_corr(1);
    tempSigma<-quad_form_diag(r,sigma);
    for (n in 1:N)
        for (s in 1:S)
            y[n,s]~multi_normal(mu[s],tempSigma);
}'''

In case you wonder, this is the subject-level part of a hierarchical model of response times and accuracy (psychology). At the bottom, trial-level I use log and logit link functions to obtain approximately bivariate normal subject level estimates based on the response time and hit indicator.

[andre.p...@googlemail.com]

Matus,

many options we can take a road. If I'm not wrong, you use a eigenvalue/eigenvector 

composition to your multi normal distribution. The common STAN options is to use

the cholesky decomposition. There's nothing wrong go either way, but cholesky 

is more calculation efficient. 

You might take the EV/EW later by D*L*(D*L)'   =: Sigma (Covariance Matrix)

and use eigen(Sigma) later. You also might avoid the need to take the sqrt(r ...

in your code. 

Andre

[matus.s...@privatdemail.net]

Hi Andre,

you are correct it's eigenvalue/eigenvector decomposition of multi-normal.

I'm not familiar with the relationship between cholesky and eigenvalue decomposition. Is D*L=U*sqrt(lambda) always true? Because if not, slope and offset - my variables of interest, have a different interpretation.

m.

[andre.p...@googlemail.com]

Matus,

have a look at:

http://stats.stackexchange.com/questions/61719/cholesky-versus-eigendecomposition-for-drawing-samples-from-a-multivariate-norma

Andre