MeasurementMYCat: Sampler for Partial Mediation Model with Multiple Categorical...

View source: R/MeasurementMYCat.R

MeasurementMYCatR Documentation

Sampler for Partial Mediation Model with Multiple Categorical Indicator for the Mediator and DV

Description

Estimates a partial mediation model with multiple categorical indicator for the mediator and the dependent variable using a mixture of Metropolis-Hastings and Gibbs sampling

Usage

MeasurementMYCat(Data,Prior,R)

Arguments

Data

list(X, m_tilde, y_tilde)

Prior

list(A_M,A_Y)

R

number of MCMC iterations, default = 10000

Details

Model

M = \beta_{0M} + X\beta_1 + U_M (eq.1)
Y = \beta_{0Y} + M\beta_2 + X\beta_3 + U_Y (eq.2)

Measurement equations

m^*_1 = M + U_{m^*_1}
\tilde{m}_1 = OrdProbit(m^*_1 ,C_{m_1})
m^*_2 = \lambda_{20} + M + U_{m^*_2}
\tilde{m}_2 = OrdProbit(m^*_2 ,C_{m_2})
...
m^*_K = \lambda_{K0} + M + U_{m^*_K}
\tilde{m}_K = OrdProbit(m^*_K ,C_{m_K})
y^*_1 = Y + U_{y^*_1}
\tilde{y}_1 = OrdProbit(y^*_1 ,C_{y_1})
y^*_2 = \tau_{20} + Y + U_{y^*_2}
\tilde{y}_2 = OrdProbit(y^*_2 ,C_{y_2})
...
y^*_L = \tau_{L0} + Y + U_{y^*_L}
\tilde{y}_L = OrdProbit(y^*_L ,C_{y_L})

Prior specification:

\beta_{0M} \sim N(0,100), \beta_{0Y} \sim N(0,100)
\beta_1 \sim N(0,100), \beta_2 \sim N(0,100), \beta_3 \sim N(0,1)
\lambda_{20},...,\lambda_{K0} \sim N(0,100)
\sigma^2_{m*_1}, ..., \sigma^2_{m*_K} \sim Inv\chi^2(\nu,\nu S)
\tau_{20}, ..., \tau_{L0} \sim N(0,100)
\sigma^2_{y*_1}, ..., \sigma^2_{y*_L} \sim Inv\chi^2(\nu,\nu S)
C*_{m_1}, ..., C*_{m_K}, C*_{y_1}, ..., C*_{y_L} \sim N(0,I)

Note: C*_{m_1}, ..., C*_{m_K}, C*_{y_1}, ..., C*_{y_L} are untransformed cutoffs, which are then exponentially transformed to impose sign and order constraint on them. Subjective prior values for the error variances are \nu=1, S=3.

Argument Details

Data = list(X, m_tilde, y_tilde)

X(N x 1)

treatment variable vector

m_tilde(N x M_ind)

mediator indicators' matrix

y_tilde(N x Y_ind)

dependent variable indicators' matrix

Prior = list(A_M,A_Y) [optional]

A_M

vector of coefficients' prior variances of eq.1, default = rep(100,2)

A_Y

vector of coefficients' prior variances of eq.2, default = c(100,100,1)

Value

beta_M(R X 2)

matrix of eq.1 coefficients' draws

beta_Y(R X 3)

matrix of eq.2 coefficients' draws

lambda (M_ind X 2 X R)

array of mediator indicator coefficients' draws. Each slice is one draw, where rows represent the indicator equation and columns the coefficients. All Slope coefficients as well as intercept of the first equation are fixed to 1 and 0 respectively.

tau (Y_ind X 2 X R)

array of dependent variable indicator coefficients' draws. Each slice is one draw, where rows represent the indicator equation and columns the coefficients. All Slope coefficients as well as intercept of the first equation are fixed to 1 and 0 respectively.

ssq_m_star(R X M_ind)

Matrix of mediator indicator equations' coefficients' error variance draws

ssq_y_star(R X Y_ind)

Matrix of dependent variable indicator equations' coefficients' error variance draws

cutoff_M (M_ind X k_M X R)

array of discretized mediator indicators' cutoff values.

cutoff_Y (Y_ind X k_Y X R)

array of discretized dependent variable indicators' cutoff values.

Mdraw(R X N)

Matrix of the augmented latent mediator

Ydraw(R X N)

Matrix of the augmented latent dependent variable

mu_draw

vector of means indexing MCMC draws of the direct effect (used in BFSD to compute Bayes factor)

var_draw

vector of variance indexing MCMC draws of the direct effect (used in BFSD to compute Bayes factor)

Examples

set.seed(60)
SimMeasurementMYCat = function(X, beta_M, cutoff_M, beta_Y, cutoff_Y, M_ind, Y_ind,
                               lambda, tau, ssq_m_star, ssq_y_star){

  nobs = dim(X)[1]
  m_tilde = m_star = matrix(double(nobs*M_ind), ncol = M_ind)
  y_tilde = y_star = matrix(double(nobs*Y_ind), ncol = Y_ind)

  M = cbind(rep(1,nobs),X)%*%beta_M + rnorm(nobs)
  for(i in 1: M_ind){
    m_star[,i] = lambda[i] + M + sqrt(ssq_m_star[i])*rnorm(nobs);
    m_tilde[,i] = cut(m_star[,i], br = cutoff_M[i,], right=TRUE, include.lowest = TRUE, labels = FALSE)
  }

  Y = cbind(rep(1,nobs),cbind(M,X))%*%beta_Y + rnorm(nobs)
  for(i in 1: Y_ind){
    y_star[,i] = tau[i] + Y + sqrt(ssq_y_star[i])*rnorm(nobs);
    y_tilde[,i] = cut(y_star[,i], br = cutoff_Y[i,], right=TRUE, include.lowest = TRUE, labels = FALSE)
  }

  return(list(Y = Y, M = M, y_tilde = y_tilde, m_tilde = m_tilde, X = X,
              k_M=dim(cutoff_M)[2]-1, beta_M = beta_M, lambda = lambda,
              ssq_m_star = ssq_m_star, m_star = m_star, cutoff_M = cutoff_M,
              k_Y=dim(cutoff_Y)[2]-1, beta_Y = beta_Y, tau = tau, ssq_y_star = ssq_y_star,
              y_star = y_star, cutoff_Y = cutoff_Y, M_ind=M_ind, Y_ind=Y_ind))
}

M_ind = 2
Y_ind = 2
Mcut = Ycut = 8
nobs=1000
X=as.matrix(runif(nobs,min=0, max=1))
beta_M = c(.5,1)
beta_Y = c(.7, 1.5, 0)
ssq_m_star = c(.5,.3)
lambda = c(0,-.5) #the intercepts for the latent M indicators w. measurement error
ssq_y_star = c(.2,.2)
tau = c(0,-.5)
cutoff_M = matrix(c(-100, 0, 1.6, 2, 2.2, 3.3, 6,  100,
                    -100, 0, 1, 2, 3, 4, 5, 100) ,ncol= Mcut, byrow = TRUE)
cutoff_Y =  matrix(c(-100, 0, 1.6, 2, 2.2, 3.3, 6,  100,
                     -100, 0, 1, 2, 3, 4, 5, 100) ,ncol= Ycut, byrow = TRUE)
DataMYCat = SimMeasurementMYCat(X, beta_M, cutoff_M, beta_Y, cutoff_Y, M_ind,
                                Y_ind, lambda, tau, ssq_m_star, ssq_y_star)


#estimation
Mcut = max(DataMYCat$m_tilde) + 1
Ycut = max(DataMYCat$y_tilde) + 1
Data = list(X=cbind(rep(1,length(DataMYCat$X)),DataMYCat$X), m_tilde=as.matrix(DataMYCat$m_tilde),
            y_tilde=as.matrix(DataMYCat$y_tilde), k_M = Mcut-1, k_Y=Ycut-1,
            M_ind=dim(as.matrix(DataMYCat$m_tilde))[2], Y_ind=dim(as.matrix(DataMYCat$y_tilde))[2])
out = MeasurementMYCat(Data=Data,R = 10000)

#results
colMeans(out$beta_M)
colMeans(out$beta_Y)
apply(out$lambdadraw,MARGIN = c(1,2),FUN = mean)
apply(out$taudraw,MARGIN = c(1,2),FUN = mean)

apply(out$cutoff_M,c(1,2),FUN = mean)
apply(out$cutoff_Y,c(1,2),FUN = mean)

arashl1364/BFMediate documentation built on Oct. 11, 2023, 5:54 p.m.