@ND - lab:zhang
https://nd.psychstat.org/
2026-04-15T03:58:49+00:00@ND
https://nd.psychstat.org/
https://nd.psychstat.org/_media/wiki/dokuwiki.svgtext/html2016-01-24T14:48:26+00:00Anonymous (anonymous@undisclosed.example.com)an_r_function_to_run_mplus_and_obtain_bc_intervals
https://nd.psychstat.org/lab/zhang/an_r_function_to_run_mplus_and_obtain_bc_intervals
## Program to run Mplus for the simulation
## function to calculate the BC for mplus
bc.ci<-function(res, est){
iB<-length(res)
## BC interval
alpha<-c(.025, .975)
z0<-qnorm(sum(res<est)/iB)
zalpha<-qnorm(alpha)
bcalpha<-pnorm(2*z0+zalpha)
BC<-quantile(res,bcalpha)
c(est, sd(res), BC)
}
## run the simulation
## run the bootstrap to get sd, ci for indirect effect
boot.mplus<-function(B, N, datafile, scriptfile, resfile){
temp.data<-read.table(datafile)
write.table(temp.da…text/html2016-01-24T14:48:26+00:00Anonymous (anonymous@undisclosed.example.com)bayesian_exploratory_factor_analysis_using_informative_priors
https://nd.psychstat.org/lab/zhang/bayesian_exploratory_factor_analysis_using_informative_priors
Bayesian exploratory factor analysis using informative priors
Data generation
model{
for (i in 1:N){
f[i,1:2]~dmnorm(muf[1:2], pre.phi[1:2,1:2])
for (j1 in 1:6){
muy[i,j1]<-l1[j1]*f[i,1]+l2[j1]*f[i,2]
y[i,j1]~dnorm(muy[i,j1],pre.p[j1])
}
}
}
list(N=200, l1=c(.6,.6,.6,0,0,0), l2=c(0,0,0, .6,.6,.6),
muf=c(0,0),
pre.p=c(1.5625,1.5625, 1.5625, 2.78, 2.78, 2.78),
pre.phi=structure(.Data =c(1.33, -0.67, -.67, 1.33), .Dim = c(2,2)))text/html2016-01-24T14:48:26+00:00Anonymous (anonymous@undisclosed.example.com)bmem_simulation_setup_using_r
https://nd.psychstat.org/lab/zhang/bmem_simulation_setup_using_r
pat<-c('comp','mcar','mcar0','mcar1','mcar2','mar','mar0','mar1','mar2','mnar','mnar0','mnar1','mnar2','mnar0-1','mnar1-1','mnar2-1')
## create folders
for (mis in pat){
system(paste('mkdir ',mis, sep=''))
for (i in 1:3){
system(paste('mkdir ',mis, '/', i, sep=''))
}
}
## copy bmem to the folders
for (mis in pat){
for (i in 1:3){
system(paste('cp bmem ',mis, '/', i,'/.', sep=''))
}
}
## generate submit script
pat1<-c('mcar','mar','mnar')
for (mis in pat1){
filename <…text/html2016-01-24T14:48:26+00:00Anonymous (anonymous@undisclosed.example.com)classes
https://nd.psychstat.org/lab/zhang/classes
* index
A nice blog on MCMC: <http://theoreticalecology.wordpress.com/2011/12/09/mcmc-chain-analysis-and-convergence-diagnostics-with-coda-in-r/>
A paper
<http://onlinelibrary.wiley.com/doi/10.1111/j.1461-0248.2011.01640.x/pdf>text/html2016-01-24T14:48:26+00:00Anonymous (anonymous@undisclosed.example.com)dynamic_mediation_analysis_with_sam
https://nd.psychstat.org/lab/zhang/dynamic_mediation_analysis_with_sam
Dynamic mediation analysis with Sam
Process the data for analysis 2010-01-25text/html2016-01-24T14:48:26+00:00Anonymous (anonymous@undisclosed.example.com)fiml-mnar-0.1
https://nd.psychstat.org/lab/zhang/fiml-mnar-0.1
1 -0.032 0.1935456 0.062 2 0.07588255 0.3975445 0.205 3 -0.009842634 0.2686099 0.105 4 0.021 0.2611601 0.109 5 -0.012 0.239 0.106 6 0.009 0.216 0.098 7 0.05 0.3628793 0.159 8 0.003 0.23304 0.087 9 0.07 0.3251475 0.164 10 -0.032 0.2302722 0.076text/html2016-01-24T14:48:26+00:00Anonymous (anonymous@undisclosed.example.com)fiml-mnar-0.2
https://nd.psychstat.org/lab/zhang/fiml-mnar-0.2
1 -0.033 0.2582820 0.067 2 0.055 0.404501 0.188 3 -0.02429877 0.2843474 0.1 4 0.028 0.3080503 0.115 5 -0.002234903 0.285089 0.117 6 -0.02 0.2295361 0.08 7 0.045 0.3777618 0.165 8 -0.001 0.2848373 0.088 9 0.069 0.3521732 0.181 10 -0.05345098 0.271 0.084text/html2016-01-24T14:48:26+00:00Anonymous (anonymous@undisclosed.example.com)fiml-mnar-0.4
https://nd.psychstat.org/lab/zhang/fiml-mnar-0.4
1 -0.06645884 0.3172354 0.056 2 -0.056 0.3131110 0.087 3 -0.07701762 0.3269593 0.081 4 -0.01032937 0.3527093 0.149 5 -0.071 0.4015737 0.117 6 -0.1919515 0.1029960 0.001 7 -0.025 0.4826391 0.096 8 -0.029 0.3291485 0.075 9 0.05379908 0.4856425 0.205text/html2016-01-24T14:48:26+00:00Anonymous (anonymous@undisclosed.example.com)generate_complete_data_for_linear_model
https://nd.psychstat.org/lab/zhang/generate_complete_data_for_linear_model
library(mvtnorm)
grm.linear<-function(N, T, R){
## Constants and parameters
mL<-50
vL<-100
mS<-1
vS<-25
vLS<-0
vE<-25
mu<-c(mL,mS)
sigma<-array(c(vL,vLS,vLS,vS), dim=c(2,2))
y<-array(NA, dim=c(N, T))
LS<-array(NA, dim=c(N, 2))
for (i in 1:N){
LS[i,]<-rmvnorm(1, mu, sigma)
for (j in 1:T){
y[i, j] <- LS[i,1] + LS[i,2]*(j-1)/(T-1) + rnorm(1, 0, sqrt(vE))
}
}
#data.file<-paste('H:\\Private\\LongPower\\type1\\data\\N-',N,'-T-',T,'-',1000+R,'.txt',sep='')
data.file<-paste('LS-N-',N,'-T-…text/html2016-01-24T14:48:26+00:00Anonymous (anonymous@undisclosed.example.com)growth_curve_modeling_with_mnar_data
https://nd.psychstat.org/lab/zhang/growth_curve_modeling_with_mnar_data
Growth curve modeling with MNAR data
An working example
#WinBUGS codes generated by BAUW: http://bauw.psychstat.org
#USE WITH CAUTION!
Model{
# Model specification for linear growth curve model
for (i in 1:N){
LS[i,1:2]~dmnorm(muLS[i,1:2], Inv_cov[1:2,1:2])
muLS[i,1]<-bL[1]
muLS[i,2]<-bS[1]
for (t in 1:5){
y[i, t] ~ dnorm(muY[i,t], Inv_Sig_e2)
muY[i,t]<-LS[i,1]+LS[i,2]*t
## missing data mechanism
m[i,t]~dbern(p[i,t])
p[i,t]<-phi(gam[t,1]+gam…text/html2016-01-24T14:48:26+00:00Anonymous (anonymous@undisclosed.example.com)individual_probability_in_mixture_model
https://nd.psychstat.org/lab/zhang/individual_probability_in_mixture_model
## Conjugate priors
model
{
for( i in 1 : N ) {
y[i] ~ dnorm(mu[Z[i]], pre.phi[Z[i]])
Z[i] ~ dcat(lambda[1:2])
logp1[i]<-log(lambda[1])-.5*log(2*3.1415926)-.5*log(phi[1])-(y[i]-mu[1])*(y[i]-mu[1])/(2*phi[1])
logp2[i]<-log(lambda[2])-.5*log(2*3.1415926)-.5*log(phi[2])-(y[i]-mu[2])*(y[i]-mu[2])/(2*phi[2])
p1[i]<-exp(logp1[i])/(exp(logp1[i])+exp(logp2[i]))
…text/html2016-01-24T14:48:26+00:00Anonymous (anonymous@undisclosed.example.com)missing_data_mediation_analysis
https://nd.psychstat.org/lab/zhang/missing_data_mediation_analysis
Simulation results
[All results in Excel table]
Data simulation
* Stars plot for result comparison
* Simulation codes on 2010-10-25
Listwise deletion
* SAS macros for listwise deletion
* R codes for running the listwise deletion simulation
Pairwise deletion
* SAS macros for pairwise deletion
* R codes for running the pairwise deletion simulation
BMEM
* Run BMEME in R for simulation
* BMEM simulation setup using R
* Process BMEM results for BMEM paper
* [Results from …text/html2016-01-24T14:48:26+00:00Anonymous (anonymous@undisclosed.example.com)paperreviewed
https://nd.psychstat.org/lab/zhang/paperreviewed
2013
Title Number Aging, Neuropsychology and Cognition 2text/html2016-01-24T14:48:26+00:00Anonymous (anonymous@undisclosed.example.com)plot_running_time_for_bmem
https://nd.psychstat.org/lab/zhang/plot_running_time_for_bmem
y<-array(c(2.05, 4.13, 15, 2.79, 5.03, 16.24, 8.41, 15.34, 50.13, 20.35, 33.37, 125.46), dim=c(3,4))
plot(c(100,200,500,1000),y[1,], ylim=c(0,130), xlab='Sample size', ylab='Time (seconds)')
lines(c(100,200,500,1000),y[1,])
points(c(100,200,500,1000),y[2,])
lines(c(100,200,500,1000),y[2,],lty=2)
points(c(100,200,500,1000),y[3,])
lines(c(100,200,500,1000),y[3,],lty=3)
## for 10%
text(100, 2.05, '2.05', pos=1)
text(200, 2.79, '2.79', pos=1)
text(500, 8.41, '8.41', pos=1)
text(1000, 20.35, '20.…text/html2016-01-24T14:48:26+00:00Anonymous (anonymous@undisclosed.example.com)power_analysis_with_timo_manuel
https://nd.psychstat.org/lab/zhang/power_analysis_with_timo_manuel
Power analysis
possible journals:
Methodology: <http://www.psycontent.com/content/120235/>
FIML
R-Mplus-MNAR-Simulation
FIML-MNAR-0.1
FIML-MNAR-0.2
FIML-MNAR-0.4
Scripts and codes
* R codes for plot two-sided Q-Q plots
* Generate complete data for linear model
Output parameter bias in a table
linear.res<-function(filename)
{
##process the results
# res.file<-paste('res-N',N,'-T',T,'-alpha',alpha,'-R',r,'-a.txt',sep='')
## filename<-'C:\\zzy\\research\\LongPower\\Results\\MAR\\l…text/html2016-01-24T14:48:26+00:00Anonymous (anonymous@undisclosed.example.com)process_bmem_results_for_bmem_paper
https://nd.psychstat.org/lab/zhang/process_bmem_results_for_bmem_paper
res<-read.table('comp.txt', fill=T)
res<-res[ !is.na(res[,65]) ,]
res<-res[, 2:65]
## population parameter values
pop<-c(0,0,.39, .39, 0, 1, 1, .1521)
n<-dim(res)[1]
for (i in 1:8){
temp<-res[, ((i-1)*8+1):(i*8)]
if (pop[i] ==0){
bias<-(mean(temp[,1])-pop[i])*100
}else{
bias<-(mean(temp[,1])-pop[i])*100/pop[i]
}
cvg1<-sum(temp[,3]<pop[i] & temp[,4]>pop[i])/n
pow1<-sum(temp[,3]>0 | temp[,4]<0)/n
cvg2<-sum(temp[,5]<pop[i] & temp[,6]>pop[i])/n
pow2<-sum(…text/html2016-01-24T14:48:26+00:00Anonymous (anonymous@undisclosed.example.com)process_the_data_for_analysis_2010-01-25
https://nd.psychstat.org/lab/zhang/process_the_data_for_analysis_2010-01-25
## add missing data for missing days
## read data
dset<-read.csv('dynamic-mediation-data-2010-01-25.csv', head=T)
uniqueid<-unique(dset$X.id)
freqid<-table(dset$X.id)
k<-length(uniqueid)
id<-rep(1:length(uniqueid), freqid)
dset<-cbind(id, dset)
m<-length(uniqueid)*50
## sort the data using id and survey
dset<-dset[order(dset$survey), ]
dset<-dset[order(dset$id), ]
dim(dset)
aug<-NULL
len<-1:50
for (i in 1:k){
temp<-dset$survey[dset$id==i]
temp1<-len[-temp]
temp2<-cbind(rep(i,length…text/html2016-01-24T14:48:26+00:00Anonymous (anonymous@undisclosed.example.com)r_codes_for_plot_of_relative_deviance_according_to_number_of_imputations
https://nd.psychstat.org/lab/zhang/r_codes_for_plot_of_relative_deviance_according_to_number_of_imputations
dset<-structure(list(V1 = c(100L, 90L, 80L, 70L, 60L, 50L, 40L, 30L,
20L, 10L), V2 = c(0, 0.003841788, -0.001492275, -0.000219173,
-0.002724556, 0.00336318, -0.012693059, -0.010003541, -0.006839342,
-0.020346107), V3 = c(0, 0.001393081, 0.002897879, 0.001972387,
0.000464792, 0.002782286, 0.002782286, -0.000697999, 0.004973973,
0.011605194), V4 = c(0, -0.004692156, -0.005776311, -0.027694902,
-0.012266019, -0.029263698, -0.052509217, -0.083641089, -0.110293602,
-0.011652035), V5 = c(0, -0.…text/html2016-01-24T14:48:26+00:00Anonymous (anonymous@undisclosed.example.com)r_codes_for_plot_two-sided_q-q_plots
https://nd.psychstat.org/lab/zhang/r_codes_for_plot_two-sided_q-q_plots
## qqplot for testing the distributions
## for linear model
qq.linear<-function(T, N, r){
T=4
N=300
r=0.1
prob<-seq(.01,.99,.01)
plotname<-paste('QQ-Linear-', T, '-N-', N, '-R-', r, '.pdf', sep='')
pdf(file=plotname, width=5, height=5)
par(mar=c(5, 4, 4, 5) + 0.1)
## complete data
file.comp<-paste('C:\\zzy\\LongPower\\Results\\qqplot\\linear-MCAR','T-',T,'-N',N,'-R',0,'.txt',sep='')
y.comp<-scan(file.comp)
plot(qchisq(prob, 1), quantile(y.comp, prob=prob), ylim=c(0,7), xlim=c(0,7),type='l…text/html2016-01-24T14:48:26+00:00Anonymous (anonymous@undisclosed.example.com)r_codes_for_processing_results_from_mplus
https://nd.psychstat.org/lab/zhang/r_codes_for_processing_results_from_mplus
## Process results from Mplus
proc.res<-function(pat, r){
filename<-paste(pat,'-',r,'.txt', sep='')
res<-read.table(filename)
res<-res[!is.na(res[,8]),]
iL<-dim(res)[1]
est<-mean(res[,6])
bias<-(est-.1521)/.1521
cvg<-sum( res[,5]<.1521 & res[,7]>.1521 )/iL
power<-sum(res[,5]>0 | res[,7]<0)/iL
cat(c(est,bias, cvg,power, iL), "\n")
}
mis.pat<-c('mcar0','mcar1','mcar2','mar0','mar1','mar2','mnar0','mnar1','mnar2','mnar0-1','mnar1-1','mnar2-1')
rate…text/html2016-01-24T14:48:26+00:00Anonymous (anonymous@undisclosed.example.com)r_codes_for_processing_results_of_mi_using_sas
https://nd.psychstat.org/lab/zhang/r_codes_for_processing_results_of_mi_using_sas
Process all results
## cp result files into one folder
pat<-c('mcar','mcar0','mcar1','mcar2','mar','mar0','mar1','mar2','mnar','mnar0','mnar1','mnar2','mnar0-1','mnar1-1','mnar2-1')
for (mis in pat){
for (i in 1:3){
cmd<-paste('cp ', mis,'/',i,'/*.txt results/.', sep='')
system(cmd)
}
}
#############################
## function to print the results
proc.res<-function(pat,r){
filename<-paste(pat,'-',r,'.txt', sep='')
res<-read.table(filename, fill=T)
res<-res[ !is.na…text/html2016-01-24T14:48:26+00:00Anonymous (anonymous@undisclosed.example.com)r_codes_for_running_mi_simulation
https://nd.psychstat.org/lab/zhang/r_codes_for_running_mi_simulation
## R codes for running simulation of pairwise deletion and listwise deletion
for (ab in c(0, .14, .39, .59)){
for (N in c(50,100,200)){
for (miss in c('MCAR','MNAR1', 'MNAR2')){
for (rate in c(.1,.2,.4) ){
for (av in c('av','noav')){
## folder name
folder.name<-paste('ab',ab,'N',N,miss,'M',rate,sep='')
## change the folder name in SAS
repNS<-paste("sed 's/simucond/",folder.name,"/g' mi.sas > temp.sas",sep='')
system(repNS)
repNS<-paste("…text/html2016-01-24T14:48:26+00:00Anonymous (anonymous@undisclosed.example.com)r_codes_for_running_mplus_for_simulation
https://nd.psychstat.org/lab/zhang/r_codes_for_running_mplus_for_simulation
* For data without auxiliary variables
## Program to run Mplus for the simulation
## write the model
gen.model<-function(datafile){
## generate mplus scripts for data analysis
## Mplus input scripts
mplus.null<-'model.inp'
## for null model
cat('TITLE: Model\n',file=mplus.null);
cat('DATA:\n',file=mplus.null, append=T);
cat(' FILE=', datafile, ';\n',file=mplus.null, append=T);
cat('VARIABLE: \n',file=mplus.null, append=T);
cat(' NAMES ARE X M Y;\n', file=mplus.null, append=T);
ca…text/html2016-01-24T14:48:26+00:00Anonymous (anonymous@undisclosed.example.com)r_codes_for_running_the_listwise_deletion_simulation
https://nd.psychstat.org/lab/zhang/r_codes_for_running_the_listwise_deletion_simulation
## R codes for running simulation of pairwise deletion and listwise deletion
for (ab in c(0, .14, .39, .59)){
for (N in c(50,100,200)){
for (miss in c('MCAR','MNAR1', 'MNAR2')){
for (rate in c(.1,.2,.4) ){
## folder name
folder.name<-paste('ab',ab,'N',N,miss,'M',rate,sep='')
## change the folder name in SAS
repNS<-paste("sed 's/simucond/",folder.name,"/g' listwise.sas > sasmodel.sas",sep='')
system(repNS)
## Run SAS for simulation
s…text/html2016-01-24T14:48:26+00:00Anonymous (anonymous@undisclosed.example.com)r_codes_for_running_the_pairwise_deletion_simulation
https://nd.psychstat.org/lab/zhang/r_codes_for_running_the_pairwise_deletion_simulation
## R codes for running simulation of pairwise deletion and listwise deletion
for (ab in c(0, .14, .39, .59)){
for (N in c(50,100,200)){
for (miss in c('MCAR','MNAR1', 'MNAR2')){
for (rate in c(.1,.2,.4) ){
## folder name
folder.name<-paste('ab',ab,'N',N,miss,'M',rate,sep='')
## change the folder name in SAS
repNS<-paste("sed 's/simucond/",folder.name,"/g' pairwise.sas > sasmodel.sas",sep='')
system(repNS)
## Run SAS for simulation
s…text/html2016-01-24T14:48:26+00:00Anonymous (anonymous@undisclosed.example.com)r-mplus-mnar-simulation
https://nd.psychstat.org/lab/zhang/r-mplus-mnar-simulation
## Program to run Mplus for the simulation
## write the model
gen.model<-function( ){
## generate mplus scripts for data analysis
## Mplus input scripts
mplus.null<-'model.inp'
## for null model
cat('TITLE: Model\n',file=mplus.null);
cat('DATA:\n',file=mplus.null, append=T);
cat(' FILE=data.txt;\n',file=mplus.null, append=T);
cat('VARIABLE: \n',file=mplus.null, append=T);
cat(' NAMES ARE X M Y U V W;\n', file=mplus.null, append=T);
cat(' USEVARIABLES ARE X M Y;\n', file=mplus.null,…text/html2016-01-24T14:48:26+00:00Anonymous (anonymous@undisclosed.example.com)robust_growth_curve_modeling_using_student-t_distribution
https://nd.psychstat.org/lab/zhang/robust_growth_curve_modeling_using_student-t_distribution
Data generation
#WinBUGS codes generated by BAUW: http://bauw.psychstat.org
#USE WITH CAUTION!
Model{
# Model specification for linear growth curve model
for (i in 1:N){
LS[i,1:2]~dmnorm(muLS[i,1:2], Inv_cov[1:2,1:2])
muLS[i,1]<-bL[1]
muLS[i,2]<-bS[1]
for (t in 1:5){
y[i, t] ~ dt(muY[i,t], Inv_Sig_e2, 4)
muY[i,t]<-LS[i,1]+LS[i,2]*t
}
}
}
# The (naive) starting values for model parameters.
list(bL=c(5), bS = c(2), Inv_Sig_e2 = 1, Inv_cov= structure(.Data…text/html2016-01-24T14:48:26+00:00Anonymous (anonymous@undisclosed.example.com)robust_mediation_analysis
https://nd.psychstat.org/lab/zhang/robust_mediation_analysis
Robust mediation analysis
WinBUGS codes
* [Normal data]
mean sd MC_error val2.5pc median val97.5pc start sample Para[1] -0.8407 0.1289 0.005785 -1.091 -0.8418 -0.5905 10001 110000 Para[2] 0.2783 0.01888 4.56E-04 0.2411 0.2783 0.3151text/html2016-01-24T14:48:26+00:00Anonymous (anonymous@undisclosed.example.com)run_bmeme_in_r_for_simulation
https://nd.psychstat.org/lab/zhang/run_bmeme_in_r_for_simulation
## run BMEM for mediation analysis
## generate batch file
gen.batch<-function(datafile){
filename<-'batch.txt'
cat("out.txt\n", file=filename)
cat(datafile, file=filename, append=T)
cat("\n", file=filename, append=T)
cat("99999\n", file=filename, append=T)
cat(".95\n", file=filename, append=T)
cat(runif(1), file=filename, append=T)
cat("\n", file=filename, append=T)
cat("1000\n", file=filename, append=T)
}
runone<-function(datafile, resfile){
for (i in 1001:2000){
file<…text/html2016-01-24T14:48:26+00:00Anonymous (anonymous@undisclosed.example.com)sas_macros_for_listwise_deletion
https://nd.psychstat.org/lab/zhang/sas_macros_for_listwise_deletion
/*** Setup the global parameters ***/
/*The parameters below should be changed accordingly*/
*%LET filename="c:\zzy\mnarmediation\active1.txt"; * data file directory and name;
%LET varname=x m y; *specify variable names in the data file. Please
use x for the input variable, m for the mediation variable, and y for
the output variable. a1 and a2 are two auxiliary variables in the data
file. You can use any names except for x, m, and y for naming the
auxiliary variables;
%LET missing=99999; *s…text/html2016-01-24T14:48:26+00:00Anonymous (anonymous@undisclosed.example.com)sas_macros_for_pairwise_deletion
https://nd.psychstat.org/lab/zhang/sas_macros_for_pairwise_deletion
/*** Setup the global parameters ***/
/*The parameters below should be changed accordingly*/
*%LET filename="c:\zzy\mnarmediation\active1.txt"; * data file directory and name;
%LET varname=x m y; *specify variable names in the data file. Please
use x for the input variable, m for the mediation variable, and y for
the output variable. a1 and a2 are two auxiliary variables in the data
file. You can use any names except for x, m, and y for naming the
auxiliary variables;
%LET missing=99999; *s…text/html2016-01-24T14:48:26+00:00Anonymous (anonymous@undisclosed.example.com)sas_macros_for_simulation
https://nd.psychstat.org/lab/zhang/sas_macros_for_simulation
/*** Setup the global parameters ***/
/*The parameters below should be changed accordingly*/
*%LET filename="c:\zzy\mnarmediation\active1.txt"; * data file directory and name;
%LET varname=var00list; *specify variable names in the data file. Please
use x for the input variable, m for the mediation variable, and y for
the output variable. a1 and a2 are two auxiliary variables in the data
file. You can use any names except for x, m, and y for naming the
auxiliary variables;
%LET missing=99999…text/html2016-01-24T14:48:26+00:00Anonymous (anonymous@undisclosed.example.com)simulation_codes_on_2010-10-25
https://nd.psychstat.org/lab/zhang/simulation_codes_on_2010-10-25
normal data
R<-1000
pc<-0
pmu1<-0
pmu2<-0
pS1<-1
pS2<-1
pSX<-1
rho<-.5
sm<-1.1521
sy<-1.175234
set.seed(20100402)
dir.create('compdata')
for (r in 1:R){
for (ab in c(0, .14, .39, .59)){
pa<-ab
pb<-ab
## Generate complete data
X<-rnorm(500)
M<-pmu1+pa*X+rnorm(500)
Y<-pmu2+pb*M+pc*X+rnorm(500)
rho<-.5
## auxiliary data
by<-sqrt(rho^2/(1-rho^2)/sy)
bm<-sqrt(rho^2/(1-rho^2)/sm)
V<-bm*M+rnorm(500)
W<-by*Y+rnorm(500)
## all data together
dset<-cbind(X, M, Y, …text/html2016-01-24T14:48:26+00:00Anonymous (anonymous@undisclosed.example.com)stars_plot_for_result_comparison
https://nd.psychstat.org/lab/zhang/stars_plot_for_result_comparison
facedata<-read.table('facedata.txt')
names(facedata)<-c('B10%', 'B20%','B40%','C10%','C20%','C40%','P10%','P20%','P40%')
stars(facedata,ncol=6, lwd=2, key.loc=c(4,4))
stars(facedata,ncol=6, lwd=2, key.loc=c(8,-1), axes=T, labels=labels)
lines(c(2,6), c(9,9))
text(4.1,9.2,'Without Auxiliary variables')
lines(c(8.5,13), c(9,9))
text(10.6,9.2,'With Auxiliary variables')
text(2,7.5,'MI')
labels<-c('MCAR/MI','MCAR/FIML','MCAR/EM','MCAR/MI','MCAR/FIML','MCAR/EM','MAR/MI','MAR/FIML','MAR/EM',…text/html2016-01-24T14:48:26+00:00Anonymous (anonymous@undisclosed.example.com)todo
https://nd.psychstat.org/lab/zhang/todo
Tasks to be finished
From Beliefs to Virtuous Behaviors
SUBMISSION WINDOW: AUGUST 1 - OCTOBER 15
* more information
* Process
Monday, Sep 6
<todo>Finish the religious data analysis</todo>
Friday, Sep 2
<todo #>Lab meeting</todo>
<todo #>Analyze religious data</todo>text/html2016-01-24T14:48:26+00:00Anonymous (anonymous@undisclosed.example.com)wainer_analysis
https://nd.psychstat.org/lab/zhang/wainer_analysis
#WinBUGS codes generated by BAUW: http://bauw.psychstat.org
#USE WITH CAUTION!
Model{
# Model specification for 1PL item response model
for (i in 1:N){
for (j in 1:20){
y[i, j] ~ dcat(p[i,j,1:4])
p[i,j, 1] <-q[i,j,1]
q[i,j, 1] <-phi(tau[j,1]- (theta[i] - beta[j])*alpha[j])
q[i,j, 2] <-phi(tau[j,2]- (theta[i] - beta[j])*alpha[j])
q[i,j, 3] <-phi(tau[j,3]- (theta[i] - beta[j])*alpha[j])
p[i,j,2]<-q[i,j,2]-q[i,j,1]
p[i,j,3]<-q[i,j,3]-q[i,j,2]
p[i,j,…