============================================================================================================
Not Yet Published.
The data set (long format) consists of four columns. Factor A (with levels \(a_i\)) and Factor B (with levels \(b_j\)) for \(i,j\in\{1,2\}\).
simData_Holds <- function(seed=277, N=50, nullB=F) {
#' Input -
#' seed: random seed (default is fixed)
#' N: the number of subjects (default is 50)
#' nullB: whether the effect B is not significant.
#'
#' Output - the data frame for a 2x2 repeated measures design (compound symmetry holds)
set.seed(seed)
if(nullB) {
## 2x2 repeated measures data simulation (no effect B)
mat <- MASS::mvrnorm(N,rep(100,4),diag(c(2,2,2,2)))
} else { #nullB==F
## 2x2 repeated measures data simulation
mat <- MASS::mvrnorm(N,c(50,52,52,54),diag(c(20,20,20,20)))
}
mat <- sweep(mat, 1, rnorm(N,sd=10), "+") #rho=10^2 / (10^2 + 20)=0.83
data.frame("subj"=paste0("s",1:N),
"score"=c(mat),
"A"=rep(rep(c("a1","a2"),2),each=N),
"B"=rep(c("b1","b2"),each=N*2),
stringsAsFactors=T)
}
df <- simData_Holds(seed=1, nullB=T)
dat <- matrix(df$score, ncol=4,
dimnames=list(paste0("s",1:length(unique(df$subj))),
c("a1b1","a2b1","a1b2","a2b2"))) #wide data format
datalist.fixed <- list("nS"=nrow(dat), "I"=diag(4), "Y"=dat,
"hA"=.5, "hB"=.5, "hAB"=.5,
"lwr"=mean(dat)-6*sd(dat),
"upr"=mean(dat)+6*sd(dat))
cov(dat) %>% #sample covariance matrix
kable(digits=3) %>% kable_styling(full_width=F)| a1b1 | a2b1 | a1b2 | a2b2 | |
|---|---|---|---|---|
| a1b1 | 119.317 | 120.747 | 119.125 | 118.294 |
| a2b1 | 120.747 | 125.045 | 122.249 | 121.861 |
| a1b2 | 119.125 | 122.249 | 122.938 | 120.576 |
| a2b2 | 118.294 | 121.861 | 120.576 | 121.600 |
model.full.fixed1 <- "
data {
int<lower=1> nS; // number of subjects
vector[4] Y[nS]; // responses
real<lower=0> hA; // scale parameter of the variance of the standardized projected treatment effect A
real<lower=0> hB; // scale parameter of the variance of the standardized projected treatment effect B
real<lower=0> hAB; // scale parameter of the variance of the standardized projected interaction effect AB
matrix[4,4] I; // scale matrix for the inverse-Wishart prior
real lwr; // lower bound for the grand mean
real upr; // upper bound for the grand mean
}
parameters {
real<lower=lwr, upper=upr> mu; // grand mean
cov_matrix[4] Sigma; // error covariance matrix
real<lower=0> gA; // variance of the standardized projected treatment effect A
real<lower=0> gB; // variance of the standardized projected treatment effect B
real<lower=0> gAB; // variance of the standardized projected interaction effect AB
real tAfix; // standardized projected treatment effect A
real tBfix; // standardized projected treatment effect B
real tABfix; // standardized projected interaction effect AB
}
transformed parameters {
real multiplier; // (square root) geometric mean of the determinant of the error covariance matrix
matrix[2,2] muS; // condition means
multiplier = pow(determinant(Sigma), 0.125);
for (i in 1:2) {
for (j in 1:2) {
muS[i,j] = mu + multiplier * (pow(-1,i+1) * tAfix / sqrt(2)
+ pow(-1,j+1) * tBfix / sqrt(2)
+ pow(-1,i+j) * tABfix / 2);
}
}
}
model {
target += uniform_lpdf(mu | lwr, upr);
target += scaled_inv_chi_square_lpdf(gA | 1, hA);
target += scaled_inv_chi_square_lpdf(gB | 1, hB);
target += scaled_inv_chi_square_lpdf(gAB | 1, hAB);
target += inv_wishart_lpdf(Sigma | 5, I);
target += normal_lpdf(tAfix | 0, sqrt(gA));
target += normal_lpdf(tBfix | 0, sqrt(gB));
target += normal_lpdf(tABfix | 0, sqrt(gAB));
for (k in 1:nS) {
target += multi_normal_lpdf(Y[k] | to_vector(muS), Sigma);
}
}
"
stanmodel.full.fixed1 <- stan_model(model_code=model.full.fixed1)
stanfit.full.fixed1 <- sampling(stanmodel.full.fixed1, data=datalist.fixed,
iter=10000, warmup=1000, chains=2, seed=277, refresh=0)
rstan::summary(stanfit.full.fixed1)[1] %>% kable(digits=4) %>% kable_classic(full_width=F)
|
model.full.fixed2 <- "
data {
int<lower=1> nS; // number of subjects
vector[4] Y[nS]; // responses
real<lower=0> hA; // scale parameter of the variance of the standardized projected treatment effect A
real<lower=0> hB; // scale parameter of the variance of the standardized projected treatment effect B
real<lower=0> hAB; // scale parameter of the variance of the standardized projected interaction effect AB
matrix[4,4] I; // scale matrix for the inverse-Wishart prior
real lwr; // lower bound for the grand mean
real upr; // upper bound for the grand mean
}
parameters {
real<lower=lwr, upper=upr> mu; // grand mean
cov_matrix[4] Sigma; // error covariance matrix
real<lower=0> gA; // variance of the standardized projected treatment effect A
real<lower=0> gB; // variance of the standardized projected treatment effect B
real<lower=0> gAB; // variance of the standardized projected interaction effect AB
real tAfix; // standardized projected treatment effect A
real tBfix; // standardized projected treatment effect B
real tABfix; // standardized projected interaction effect AB
}
transformed parameters {
real multiplier; // using the trace of the error covariance matrix
matrix[2,2] muS; // condition means
multiplier = sqrt(trace(Sigma))/2;
for (i in 1:2) {
for (j in 1:2) {
muS[i,j] = mu + multiplier * (pow(-1,i+1) * tAfix / sqrt(2)
+ pow(-1,j+1) * tBfix / sqrt(2)
+ pow(-1,i+j) * tABfix / 2);
}
}
}
model {
target += uniform_lpdf(mu | lwr, upr);
target += scaled_inv_chi_square_lpdf(gA | 1, hA);
target += scaled_inv_chi_square_lpdf(gB | 1, hB);
target += scaled_inv_chi_square_lpdf(gAB | 1, hAB);
target += inv_wishart_lpdf(Sigma | 5, I);
target += normal_lpdf(tAfix | 0, sqrt(gA));
target += normal_lpdf(tBfix | 0, sqrt(gB));
target += normal_lpdf(tABfix | 0, sqrt(gAB));
for (k in 1:nS) {
target += multi_normal_lpdf(Y[k] | to_vector(muS), Sigma);
}
}
"
stanmodel.full.fixed2 <- stan_model(model_code=model.full.fixed2)
stanfit.full.fixed2 <- sampling(stanmodel.full.fixed2, data=datalist.fixed,
iter=10000, warmup=1000, chains=2, seed=277, refresh=0)
rstan::summary(stanfit.full.fixed2)[1] %>% kable(digits=4) %>% kable_classic(full_width=F)
|
model.full.fixed3 <- "
data {
int<lower=1> nS; // number of subjects
vector[4] Y[nS]; // responses
real<lower=0> hA; // scale parameter of the variance of the standardized projected treatment effect A
real<lower=0> hB; // scale parameter of the variance of the standardized projected treatment effect B
real<lower=0> hAB; // scale parameter of the variance of the standardized projected interaction effect AB
matrix[4,4] I; // scale matrix for the inverse-Wishart prior
real lwr; // lower bound for the grand mean
real upr; // upper bound for the grand mean
}
parameters {
real<lower=lwr, upper=upr> mu; // grand mean
cov_matrix[4] Sigma; // error covariance matrix
real<lower=0> gA; // variance of the standardized projected treatment effect A
real<lower=0> gB; // variance of the standardized projected treatment effect B
real<lower=0> gAB; // variance of the standardized projected interaction effect AB
real tAfix; // standardized projected treatment effect A
real tBfix; // standardized projected treatment effect B
real tABfix; // standardized projected interaction effect AB
}
transformed parameters {
real multiplier; // using the Frobenius norm of the error covariance matrix
matrix[2,2] muS; // condition means
multiplier = pow(sum(square(to_vector(Sigma))), 0.25)/sqrt(2);
for (i in 1:2) {
for (j in 1:2) {
muS[i,j] = mu + multiplier * (pow(-1,i+1) * tAfix / sqrt(2)
+ pow(-1,j+1) * tBfix / sqrt(2)
+ pow(-1,i+j) * tABfix / 2);
}
}
}
model {
target += uniform_lpdf(mu | lwr, upr);
target += scaled_inv_chi_square_lpdf(gA | 1, hA);
target += scaled_inv_chi_square_lpdf(gB | 1, hB);
target += scaled_inv_chi_square_lpdf(gAB | 1, hAB);
target += inv_wishart_lpdf(Sigma | 5, I);
target += normal_lpdf(tAfix | 0, sqrt(gA));
target += normal_lpdf(tBfix | 0, sqrt(gB));
target += normal_lpdf(tABfix | 0, sqrt(gAB));
for (k in 1:nS) {
target += multi_normal_lpdf(Y[k] | to_vector(muS), Sigma);
}
}
"
stanmodel.full.fixed3 <- stan_model(model_code=model.full.fixed3)
stanfit.full.fixed3 <- sampling(stanmodel.full.fixed3, data=datalist.fixed,
iter=10000, warmup=1000, chains=2, seed=277, refresh=0)
rstan::summary(stanfit.full.fixed3)[1] %>% kable(digits=4) %>% kable_classic(full_width=F)
|