Learn R Programming

macc (version 1.0.1)

sim.data.multi: Generate two/three-level simulation data

Description

This function generates a two/three-level dataset with given parameters.

Usage

sim.data.multi(Z.list, N, K = 1, Theta, Sigma, 
    Psi = diag(rep(1, 3)), Lambda = diag(rep(1, 3)))

Arguments

Z.list

a list of data. When K = 1 (a two-level dataset), each list is a vector containing the treatment/exposure assignment of the trials for each subject; When K > 1 (a three-level dataset), each list is a list of length K, and each contains a vector of treatment/exposure assignment.

N

an integer, indicates the number of subjects.

K

an integer, indicates the number of sessions of each subject.

Theta

a 2 by 2 matrix, containing the population level model coefficients.

Sigma

a 2 by 2 matrix, is the covariance matrix of the model errors in the single-level model.

Psi

the covariance matrix of the random effects in the mixed effects model of the model coefficients. This is used only when K > 1. Default is a 3-dimensional identity matrix.

Lambda

the covariance matrix of the model errors in the mixed effects model if K > 1 or the linear model if K = 1 of the model coefficients.

Value

data

a list of data. When K = 1, each list is a contains a dataframe; when K > 1, each list is a list of length K, and within each list is a dataframe.

A

the value of \(A\)s. When K = 1, it is a vector of length N; when K > 1, it is a N by K matrix.

B

the value of \(B\)s. When K = 1, it is a vector of length N; when K > 1, it is a N by K matrix.

C

the value of \(C\)s. When K = 1, it is a vector of length N; when K > 1, it is a N by K matrix.

type

a character indicates the type of the dataset. When K = 1, type = twolevel; when K > 1, type = multilevel

Details

When K > 1 (three-level data), for the \(n_{ik}\) trials in each session of each subject, the single level mediation model is $$M_{ik}=Z_{ik}A_{ik}+E_{1ik},$$ $$R_{ik}=Z_{ik}C_{ik}+M_{ik}B_{ik}+E_{2ik},$$ where \(Z_{ik}\), \(M_{ik}\), \(R_{ik}\), \(E_{1ik}\), and \(E_{2ik}\) are vectors of length \(n_{ik}\). Sigma is the covariance matrix of \((E_{1ik},E_{2ik})\) (for simplicity, Sigma is the same across sessions and subjects). For coefficients \(A_{ik}\), \(B_{ik}\) and \(C_{ik}\), we assume a mixed effects model. The random effects are from a trivariate normal distribution with mean zero and covariance Psi; and the model errors are from a trivariate normal distribution with mean zero and covariance Lambda. For the fixed effects \(A\), \(B\) and \(C\), the values are specified in Theta with Theta[1,1] = A, Theta[1,2] = C and Theta[2,2] = B.

When K = 1 (two-level data), the single-level model coefficients \(A_{i}\), \(B_{i}\) and \(C_{i}\) are assumed to follow a trivariate linear regression model, where the population level coefficients are specified in Theta and the model errors are from a trivariate normal distribution with mean zero and covariance Lambda.

See Section 5.2 of the reference for details.

References

Zhao, Y., & Luo, X. (2014). Estimating Mediation Effects under Correlated Errors with an Application to fMRI. arXiv preprint arXiv:1410.7217.

Examples

Run this code
# NOT RUN {
   ###################################################
    # Generate a two-level dataset

    # covariance matrix of errors
    delta<-0.5
    Sigma<-matrix(c(1,2*delta,2*delta,4),2,2)

    # model coefficients
    A0<-0.5
    B0<--1
    C0<-0.5

    Theta<-matrix(c(A0,0,C0,B0),2,2)

    # number of subjects, and trials of each
    set.seed(2000)
    N<-50
    K<-1
    n<-matrix(NA,N,K)
    for(i in 1:N)
    {
      n0<-rpois(1,100)
      n[i,]<-rpois(K,n0)
    }

    # treatment assignment list
    set.seed(100000)
    Z.list<-list()
    for(i in 1:N)
    {
      Z.list[[i]]<-rbinom(n[i,1],size=1,prob=0.5)
    }

    # Lambda
    rho.AB=rho.AC=rho.BC<-0
    Lambda<-matrix(0,3,3)
    lambda2.alpha=lambda2.beta=lambda2.gamma<-0.5
    Lambda[1,2]=Lambda[2,1]<-rho.AB*sqrt(lambda2.alpha*lambda2.beta)
    Lambda[1,3]=Lambda[3,1]<-rho.AC*sqrt(lambda2.alpha*lambda2.gamma)
    Lambda[2,3]=Lambda[3,2]<-rho.BC*sqrt(lambda2.beta*lambda2.gamma)
    diag(Lambda)<-c(lambda2.alpha,lambda2.beta,lambda2.gamma)

    # Data
    set.seed(5000)
    re.dat<-sim.data.multi(Z.list=Z.list,N=N,K=K,Theta=Theta,Sigma=Sigma,Lambda=Lambda)
    data2<-re.dat$data
    ###################################################

    ###################################################
    # Generate a three-level dataset

    # covariance matrix of errors
    delta<-0.5
    Sigma<-matrix(c(1,2*delta,2*delta,4),2,2)

    # model coefficients
    A0<-0.5
    B0<--1
    C0<-0.5

    Theta<-matrix(c(A0,0,C0,B0),2,2)

    # number of subjects, and trials of each
    set.seed(2000)
    N<-50
    K<-4
    n<-matrix(NA,N,K)
    for(i in 1:N)
    {
      n0<-rpois(1,100)
      n[i,]<-rpois(K,n0)
    }

    # treatment assignment list
    set.seed(100000)
    Z.list<-list()
    for(i in 1:N)
    {
      Z.list[[i]]<-list()
      for(j in 1:K)
      {
        Z.list[[i]][[j]]<-rbinom(n[i,j],size=1,prob=0.5)
      }
    }

    # Psi and Lambda
    sigma2.alpha=sigma2.beta=sigma2.gamma<-0.5
    theta2.alpha=theta2.beta=theta2.gamma<-0.5
    Psi<-diag(c(sigma2.alpha,sigma2.beta,sigma2.gamma))
    Lambda<-diag(c(theta2.alpha,theta2.beta,theta2.gamma))

    # Data
    set.seed(5000)
    re.dat<-sim.data.multi(Z.list,N,K,Theta,Sigma,Psi,Lambda)
    data3<-re.dat$data
    ###################################################
# }

Run the code above in your browser using DataLab