Dear all,
I’m trying to make a model with human data (31 subjects, each ~ 180 trials).
There are two steps:
- For each person predict their behaviour on each trial of the task as a function of different factors (regressors) in a logistic regression.
- Measure the correlation between one of the regression weights (in the example here number 2) and an independent questionnaire measure (that is read in as ‘data’)
I have maybe found a solution, but I’m not certain that it is correct. What I’m wondering about specifically is how I’m using the regression weight (regressor 2) for each person to correlate it with the questionnaire data; whether it is wrong to put the ‘data’ together with an estimated parameter as a ‘transformed parameter’ to put it into the multi_normal? In the model fit, maybe not surprisingly, there are only 2 ‘samples’ for the part of the transformed parameter that are given as data (and some rhats for this are ‘nan’). Everything else has rhats of about 1.0 and many samples. [I’m also uploading the .stan file because while the formatting looks good in the screen where I enter it, it doesn’t look good in the preview).
data {
int<lower=1> N; //number of data points
int<lower=1> Nsub; //number of participants
int<lower=1> Npreds; //number of regressors/predictors (count excludes the constant)
int subI[N]; //which person a trial bleongs to
real QuesScores[Nsub]; //independent questionnaire measurement for each person
int<lower=0,upper=1> respMat[N]; //participants' responses(0 or 1 - independent variable)
row_vector[Npreds+1] datMat[N]; //predictors for regression (also has a constant)
}
parameters {
// Group level
vector[Npreds+1] beta_mu;
vector<lower=0>[Npreds+1] beta_sd;
// Individual people
vector[Npreds+1] beta_pr[Nsub];
// For correlation
real regQuesMu;
vector<lower=0>[1] sigma;//[2] sigma;
real<lower=-1,upper=1> r; // correlation
}
transformed parameters{
vector[2] combRegQues[Nsub]; // ??putting together the estimated parameter and the questionnaire score (to later use with multinormal distribution)
vector[2] combRegQuesMu;
vector[Npreds+1] beta[Nsub];
cov_matrix[2] T;
// transform regression weights (non-centered distr. to speed up)
for (is in 1:Nsub){
beta[is] = beta_pr[is].*beta_sd+beta_mu;
}
// Add in the questionnaires??
combRegQuesMu[1] = beta_mu[2];
combRegQuesMu[2] = regQuesMu;
for (is in 1:Nsub){
combRegQues[is,1] = beta[is,2];
combRegQues[is,2] = QuesScores[is];
}
// fill the covariance matrix
T[1,1] = sigma[1]*sigma[1];
T[1,2] = r * sigma[1] * beta_sd[2];//sigma[2];
T[2,1] = T[1,2];
T[2,2] = beta_sd[2]*beta_sd[2];//sigma[2]*sigma[2];
}
// This block runs the actual model
model {
real mu[N];
// Define the priors
beta_mu ~ cauchy(0,5);
beta_mu[1] ~ cauchy(0,20);
beta_sd ~ cauchy(0,5);
regQuesMu ~ cauchy(0,5); // i have zscored the input
for (is in 1:Nsub){
beta_pr[is] ~ normal(0,1);//(beta_mu,beta_sd);
}
// For each trial, combine all the predictors and outcomes
for (i in 1:N){
mu[i] = datMat[i]*beta[subI[i]];
}
// Describe the regression model (i.e. how data and parameters relate)
respMat ~ bernoulli_logit(mu);
// Correlation between regression weight and questionnaire score
combRegQues ~ multi_normal(combRegQuesMu,T);
}
temp.data.R (833.5 KB)
ToyExamplCorrRegressorQuestionnaireScoreV2.stan (2.4 KB)
ToyExamplCorrRegressorQuestionnaireScoreV2.stan (2.4 KB)
[edit: escaped code]