View source: R/MeasurementMYCat.R
MeasurementMYCat | R Documentation |
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
MeasurementMYCat(Data,Prior,R)
Data |
list(X, m_tilde, y_tilde) |
Prior |
list(A_M,A_Y) |
R |
number of MCMC iterations, default = 10000 |
M = \beta_{0M} + X\beta_1 + U_M
(eq.1)
Y = \beta_{0Y} + M\beta_2 + X\beta_3 + U_Y
(eq.2)
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})
\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.
Data = list(X, m_tilde, y_tilde)
treatment variable vector
mediator indicators' matrix
dependent variable indicators' matrix
Prior = list(A_M,A_Y)
[optional]vector of coefficients' prior variances of eq.1, default = rep(100,2)
vector of coefficients' prior variances of eq.2, default = c(100,100,1)
matrix of eq.1 coefficients' draws
matrix of eq.2 coefficients' draws
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.
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.
Matrix of mediator indicator equations' coefficients' error variance draws
Matrix of dependent variable indicator equations' coefficients' error variance draws
array of discretized mediator indicators' cutoff values.
array of discretized dependent variable indicators' cutoff values.
Matrix of the augmented latent mediator
Matrix of the augmented latent dependent variable
vector of means indexing MCMC draws of the direct effect (used in BFSD to compute Bayes factor)
vector of variance indexing MCMC draws of the direct effect (used in BFSD to compute Bayes factor)
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)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.