Outline Outline Introduction Introduction Using R as a Wrapper - - PowerPoint PPT Presentation

Using R as a Wrapper in Using R as a Wrapper in Simulation Studies Simulation Studies

Tim F Liao Tim F Liao University of Illinois University of Illinois

Outline Outline

• Introduction

Introduction

• Procedure for External Calls

Procedure for External Calls

• R Code

R Code

• An Example

An Example

• Conclusion

Conclusion

Introduction Introduction

• Monte Carlo simulation is a useful method for

Monte Carlo simulation is a useful method for assessing statistical robustness assessing statistical robustness

• Various specialty statistical software packages

Various specialty statistical software packages are popular among data analysts are popular among data analysts

• These programs are often limited in their

These programs are often limited in their capacities to perform simulation capacities to perform simulation

• R as a general

R as a general-

• purpose program that is

purpose program that is approaching the status of lingua franca, has approaching the status of lingua franca, has much to offer much to offer

• R can be used as a wrapper to call external

R can be used as a wrapper to call external programs to carry out simulation programs to carry out simulation

Procedure Procedure

• A: Initial call to run an external program

A: Initial call to run an external program

• Run a LCA model to simulate data

Run a LCA model to simulate data

• Estimate a model of simulated data

Estimate a model of simulated data

• B: Collect necessary output

B: Collect necessary output

• Check if output read is indeed output wanted

Check if output read is indeed output wanted

• Collect output in a single data matrix

Collect output in a single data matrix

• C: Monte Carlo simulation; repeat A & B a large

C: Monte Carlo simulation; repeat A & B a large number of times number of times

• Conduct post

Conduct post-

• simulation analysis of the output

simulation analysis of the output

• Draw violin plots of parameter estimates

Draw violin plots of parameter estimates

• Any other analyses

Any other analyses

R Code R Code

#test run of 1000 simulations #test run of 1000 simulations T < T <-

• 1000

1000 C < C <-

• 30

30 zM zM< <-

• matrix(rep(0,C*

matrix(rep(0,C*T),nrow T),nrow=T) =T) for(i for(i in (1:T)) { in (1:T)) { #the #the zM zM simulations simulations system("c:/temp/simlem/lem system("c:/temp/simlem/lem c:/temp/simlem/SimM1a.inp c:/temp/simlem/SimM1a.inp c:/temp/simlem/SimM.out c:/temp/simlem/SimM.out") ") system("c:/temp/simlem/lem system("c:/temp/simlem/lem c:/temp/simlem/Sim5E1.inp c:/temp/simlem/Sim5E1.inp c:/temp/simlem/Sim5E.out") c:/temp/simlem/Sim5E.out") a1< a1<-

• read.table("c:/temp/simlem/simla.out")[4,]

read.table("c:/temp/simlem/simla.out")[4,] … … … … b4< b4<-

• read.table("c:/temp/simlem/simla.out")[17,]

read.table("c:/temp/simlem/simla.out")[17,] b5< b5<-

• read.table("c:/temp/simlem/simla.out")[21,]

read.table("c:/temp/simlem/simla.out")[21,] a < a <-

• a1+a2+a3+a4+a5

a1+a2+a3+a4+a5 b < b <-

• b1+b2+b3+b4+b5

b1+b2+b3+b4+b5 if (a<b) { if (a<b) { zM[i,1]< zM[i,1]<-

• read.table("c:/temp/simlem/simla.out")[2,]

read.table("c:/temp/simlem/simla.out")[2,] } } else { else { zM[i,1]< zM[i,1]<-

• read.table("c:/temp/simlem/simla.out")[3,]

read.table("c:/temp/simlem/simla.out")[3,] } } … … … … } }

Example: Latent Class Analysis Example: Latent Class Analysis

• Assuming for the data the model structure:

Assuming for the data the model structure:

• Varying

Varying

• presence of non

presence of non-

• na

naï ïve MNAR association ve MNAR association

• true latent class proportions

true latent class proportions

• sample size

sample size

• Using the model to estimate

Using the model to estimate

( )

V V u u u u u u u u u u u u u m

r x XE xe XD xd XC xc XB xb XA xa j R j E e D d C c B b A a X x xabcder

− + + + + + + + + + + + + + =

∑

* * *

log β

π=0.1,β=0 π=0.1,β=1 π=0.2,β=0 π=0.2,β=1 π=0.3,β=0 π=0.3,β=1 π=0.4,β=0 π=0.4,β=1 π=0.5,β=0 π=0.5,β=1 40 60 80 100 120 140

Violin Plots of Likelihood Ratio of 1,000 Simulations (N=200) by True Latent Class Proportion and Association Parameter

Likelihood Ratio Statistic π=0.1,β=0 π=0.1,β=1 π=0.2,β=0 π=0.2,β=1 π=0.3,β=0 π=0.3,β=1 π=0.4,β=0 π=0.4,β=1 π=0.5,β=0 π=0.5,β=1 −560 −540 −520 −500 −480

Violin Plots of BIC of 1,000 Simulations (N=200) by True Latent Class Proportion and Association Parameter

BIC

π=0.1,β=0 π=0.1,β=1 π=0.2,β=0 π=0.2,β=1 π=0.3,β=0 π=0.3,β=1 π=0.4,β=0 π=0.4,β=1 π=0.5,β=0 π=0.5,β=1 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7

Violin Plots of Latent Probability of 1,000 Simulations (N=200) by True Latent Class Proportion and Association Parameter

Latent Class Probability

Conclusion Conclusion

• R is the most flexible in facilitating

R is the most flexible in facilitating simulation using existing specialty simulation using existing specialty software software

• The collection of necessary output can

The collection of necessary output can

• nly be ad hoc
• nly be ad hoc
• The procedure can apply to other popular

The procedure can apply to other popular statistical software packages such as statistical software packages such as MPlus MPlus, which can also run in dos mode , which can also run in dos mode