@ND - lab:projects:02robust_growth_curve_modeling_using_student-t_distribution
https://nd.psychstat.org/
2026-05-08T22:50:36+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)a_working_simulation_example
https://nd.psychstat.org/lab/projects/02robust_growth_curve_modeling_using_student-t_distribution/a_working_simulation_example
Analyze data as student-t data
mean sd MC_error val2.5pc median val97.5pc start sample Inv_Sig_e2 0.9794 0.05691 0.002632 0.8701 0.9795 1.093 500 2511 Inv_cov[11] 1.135 0.227 0.02002 0.7734 1.103 1.663 500 2511 Inv_cov[12]text/html2016-01-24T14:48:26+00:00Anonymous (anonymous@undisclosed.example.com)alternative_specification_using_auxiliary_variables
https://nd.psychstat.org/lab/projects/02robust_growth_curve_modeling_using_student-t_distribution/alternative_specification_using_auxiliary_variables
#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]
Inv_Sig_e2k[i]<-Inv_Sig_e2*w[i]
for (t in 1:5){
y[i, t] ~ dnorm(muY[i,t], Inv_Sig_e2k[i])
muY[i,t]<-LS[i,1]+LS[i,2]*t
}
w[i]~dgamma(k2, k2)
}
k2<-k/2
k<-1/m
m~dunif(0,.5)
#Priors for model parameter
for …text/html2016-01-24T14:48:26+00:00Anonymous (anonymous@undisclosed.example.com)bugs_codes_for_data_simulation
https://nd.psychstat.org/lab/projects/02robust_growth_curve_modeling_using_student-t_distribution/bugs_codes_for_data_simulation
Data generation
Model{
# Model specification for linear growth curve model
for (i in 1:300){
LS[i,1:2]~dmnorm(muLS[i,1:2], Inv_cov[1:2,1:2])
muLS[i,1]<-6
muLS[i,2]<-1
for (t in 1:5){
y[i, t] ~ dnorm(muY[i,t], Inv_Sig_e2[i])
muY[i,t]<-LS[i,1]+LS[i,2]*(t-1)/4
}
w[i]~dgamma(a, a)
Inv_Sig_e2[i]<-w[i]*Pre_e2
}
a<-1.5
Pre_e2<- 2
}
# The (naive) starting values for model parameters.
list(Inv_cov= structure(.Data = c(.5,0,0,2),.Dim=c(2,2)))
Model{
…text/html2016-01-24T14:48:26+00:00Anonymous (anonymous@undisclosed.example.com)gcm_with_wishart_prior
https://nd.psychstat.org/lab/projects/02robust_growth_curve_modeling_using_student-t_distribution/gcm_with_wishart_prior
#WinBUGS codes generated by BAUW: http://bauw.psychstat.org
#USE WITH CAUTION!
Model{
# Model specification for linear growth curve model
for (i in 1:NS){
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[i])
muY[i,t]<-LS[i,1]+LS[i,2]*A[t]
# muY[i,t]<-LS[i,1]+LS[i,2]*t
}
w[i]~dgamma(a, a)
Inv_Sig_e2[i]<-w[i]*Pre_e2
}
A[1]<-0
A[2]~dnorm(0, 1.0E-2)
A[3]~dnorm(0,…text/html2016-01-24T14:48:26+00:00Anonymous (anonymous@undisclosed.example.com)generate_student_t_random_numbers_using_auxiliary_variable_method_in_r
https://nd.psychstat.org/lab/projects/02robust_growth_curve_modeling_using_student-t_distribution/generate_student_t_random_numbers_using_auxiliary_variable_method_in_r
## student t random number generator
## using data augmentation method
n<-10000
df = 5
mu =0
sig = 1
w<-rep(NA, n)
t.y<-rep(NA, n)
for (i in 1:n){
#w[i]<-rchisq(1,df)
w[i]<-rgamma(1, df/2, df/2)
t.y[i]<-rnorm(1, mu, sig/sqrt(w[i]))
}
## Q-Q plot
prob<-seq(.01, .99, .01)
plot(qt(prob, df), quantile(t.y, prob))
qqplot(t.y, rt(n, df))text/html2016-01-24T14:48:26+00:00Anonymous (anonymous@undisclosed.example.com)index
https://nd.psychstat.org/lab/projects/02robust_growth_curve_modeling_using_student-t_distribution/index
Manuscripts
* Version 1.
Read data analysis
Data
* NLSY real data
* R codes related to real data analysis
* [NLSY digit span variable (???)]
Simulations
* BUGS codes for data simulation
* R codes to process simulated data from BUGS
* R codes to generate simulation submission scripts
* R codes to process the coda files
* R codes for final results process
* R function for processing DIC
Using SAS to run simulation
* SAS codes
* Latent basis GCM in bugs
* Initial…text/html2016-01-24T14:48:26+00:00Anonymous (anonymous@undisclosed.example.com)latent_basis_gcm_in_bugs
https://nd.psychstat.org/lab/projects/02robust_growth_curve_modeling_using_student-t_distribution/latent_basis_gcm_in_bugs
#WinBUGS codes generated by BAUW: http://bauw.psychstat.org
#USE WITH CAUTION!
Model{
# Model specification for linear growth curve model
for (i in 1:NS){
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[i])
muY[i,t]<-LS[i,1]+LS[i,2]*A[t]
# muY[i,t]<-LS[i,1]+LS[i,2]*t
}
w[i]~dgamma(a, a)
Inv_Sig_e2[i]<-w[i]*Pre_e2
}
A[1]<-0
A[2]~dnorm(0, 1.0E-2)
A[3]~dnorm(0,…text/html2016-01-24T14:48:26+00:00Anonymous (anonymous@undisclosed.example.com)nlsy_real_data
https://nd.psychstat.org/lab/projects/02robust_growth_curve_modeling_using_student-t_distribution/nlsy_real_data
Mplus format
6.1 5 4.9 5 .
5.9 7.5 6.9 7 7.6
4.9 5.2 5.4 6.2 5.9
5.4 4.1 5.9 6.8 3.8
6.5 5.9 6.6 7.3 .
6.1 5.4 5.6 5.5 6.5
6.3 7.2 6.4 . 6.9
5.8 4.9 5.4 3.9 5
5.2 5.6 5.6 10 7.6
2.4 3 2.7 2.7 2
9.3 9.3 9.5 9.7 10
6.3 5.3 9 8.5 8.7
8 8.1 7.6 9.9 8.8
6.2 6.9 . 4.9 .
5.4 3.8 5.3 . 5.9
9.2 9.3 9.7 9.7 10
7.4 6.4 7.8 8 7.3
5.6 . 4.3 6.8 8
4.2 5 5.7 5.7 6.2
8 8.5 8.9 9.2 9.7
5.1 5.7 . . 6
7.5 7.8 9.2 9.4 9.4
6 3.7 4.2 5.2 5.3
8 6.5 8.5 8.6 8.6
6.5 6.9 . 7.9 8.3
5.6 5.1 6.1 6 5.5
4.3 7.8 7 9 6.9
6 6.…text/html2016-01-24T14:48:26+00:00Anonymous (anonymous@undisclosed.example.com)openbugs_codes_for_normal_error_analysis
https://nd.psychstat.org/lab/projects/02robust_growth_curve_modeling_using_student-t_distribution/openbugs_codes_for_normal_error_analysis
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
}
}
#Priors for model parameter
for (i in 1:1){
bL[i] ~ dnorm(0, 1.0E-6)
bS[i] ~ dnorm(0, 1.0E-6)
}
Inv_cov[1:2,1:2]~dwish(R[1:2,1:2], 2)
R[1,1]<-1
R[2,2]<-1
R[2,1]<-R[1,2]
R[1,2]<-0
Inv_Sig_e2 ~ …text/html2016-01-24T14:48:26+00:00Anonymous (anonymous@undisclosed.example.com)openbugs_codes_for_student_t_error_analysis_with_exponential_prior_lambda_.1._01_.001
https://nd.psychstat.org/lab/projects/02robust_growth_curve_modeling_using_student-t_distribution/openbugs_codes_for_student_t_error_analysis_with_exponential_prior_lambda_.1._01_.001
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, k)
muY[i,t]<-LS[i,1]+LS[i,2]*t
}
}
k~dexp(.1)
#Priors for model parameter
for (i in 1:1){
bL[i] ~ dnorm(0, 1.0E-6)
bS[i] ~ dnorm(0, 1.0E-6)
}
Inv_cov[1:2,1:2]~dwish(R[1:2,1:2], 2)
R[1,1]<-1
R[2,2]<-1
R[2,1]<-R[1,2]
R[1,2]<-0
…text/html2016-01-24T14:48:26+00:00Anonymous (anonymous@undisclosed.example.com)openbugs_codes_for_student_t_error_analysis_with_uniform_prior
https://nd.psychstat.org/lab/projects/02robust_growth_curve_modeling_using_student-t_distribution/openbugs_codes_for_student_t_error_analysis_with_uniform_prior
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, k)
muY[i,t]<-LS[i,1]+LS[i,2]*t
}
}
k~dunif(2, 30)
#Priors for model parameter
for (i in 1:1){
bL[i] ~ dnorm(0, 1.0E-6)
bS[i] ~ dnorm(0, 1.0E-6)
}
Inv_cov[1:2,1:2]~dwish(R[1:2,1:2], 2)
R[1,1]<-1
R[2,2]<-1
R[2,1]<-R[1,2]
R[1,2]<…text/html2016-01-24T14:48:26+00:00Anonymous (anonymous@undisclosed.example.com)r_codes_for_final_results_process
https://nd.psychstat.org/lab/projects/02robust_growth_curve_modeling_using_student-t_distribution/r_codes_for_final_results_process
## Calculate bias, cvg, type1 error, and power from simulation results
## For student t analysis
proc.res<-function(resname){
res<-read.table(resname,header=T,row.names=NULL)
N<-dim(res)[1]
par.name<-c('mL','vL','mS','vS','vLS','vE','k')
par<-c(6,3,2,1,0,1,3)
for (i in 0:6){
ind<-c(i+2, i+9, i+16, i+23, 2*i+30, 2*i+31)
temp<-res[,ind]
avg<-mean(temp[,1])
if (par[i+1]==0){
bias<-(avg-par[i+1])*100
}else{
bias<-(avg-par[i+1])/par[i+1]*100
}
esd<-sd(temp[,1]…text/html2016-01-24T14:48:26+00:00Anonymous (anonymous@undisclosed.example.com)r_codes_related_to_real_data_analysis
https://nd.psychstat.org/lab/projects/02robust_growth_curve_modeling_using_student-t_distribution/r_codes_related_to_real_data_analysis
## coda analysis
student<-read.coda('coda-uniform.txt','index-uniform.txt')
geweke.diag(student,.2,.5)
HPDinterval(student)
## history plot
## plot coda
coda<-read.table('/Volumes/~zzhang4/Private/zzy/research/Bayesian analysis of longitudinal data/Student t distribution/real-data-analysis/latent basis coda chain.txt')
par<-matrix(c(coda[,2]),ncol=11)
## history plot
par(mfrow=c(5,2), mar=c(3, 5, 1, 1) + 0.1)
plot.ts(par[,1], ylab='A[2]')
plot.ts(par[,2], ylab='A[3]')
plot.ts(par[,3], ylab='…text/html2016-01-24T14:48:26+00:00Anonymous (anonymous@undisclosed.example.com)r_codes_to_generate_simulation_submission_scripts
https://nd.psychstat.org/lab/projects/02robust_growth_curve_modeling_using_student-t_distribution/r_codes_to_generate_simulation_submission_scripts
Uniform priors
## for the analysis with uniform priors
## The results are saved in the folder of *R
gen.script<-function(DATA, N){
folder<-paste(DATA,'R', sep='')
folder1<-paste(DATA, sep='')
NSS<-paste('N',N,sep='')
## create a folder for a given sample size
mkNSS<-paste('mkdir ../',folder,'/',NSS, sep='')
system(mkNSS)
## replace "NS" in model.txt by the sample size
repNS<-paste("sed 's/NS/",N,"/g' model.txt > tempmodel.txt",sep='')
system(repNS)
for (i in 1:10){
## create a folder
mkf…text/html2016-01-24T14:48:26+00:00Anonymous (anonymous@undisclosed.example.com)r_codes_to_process_simulated_data_from_bugs
https://nd.psychstat.org/lab/projects/02robust_growth_curve_modeling_using_student-t_distribution/r_codes_to_process_simulated_data_from_bugs
## Simulated t data
y<-read.table('data-generated-student-coda.txt')
ymat<-array(y[, 2], dim=c(1000,5, 500))
## write data into txt files
for (i in 1001:2000){
t(ymat[i-1000,,])->temp
filename<-paste('./Tdata/TdataN500I',i,'.txt', sep='')
write.table(temp, filename, row.names=F, col.names=F)
}
source("http://nd.psychstat.org/_media/classdata/rtobugs.r")
for (i in 1001:2000){
t(ymat[i-1000,,])->temp
filename<-paste('./Tdata/TdataN500I',i,'bugs.txt', sep='')
rtobugs(list(N=500, y=temp),…text/html2016-01-24T14:48:26+00:00Anonymous (anonymous@undisclosed.example.com)r_codes_to_process_the_coda_files
https://nd.psychstat.org/lab/projects/02robust_growth_curve_modeling_using_student-t_distribution/r_codes_to_process_the_coda_files
2010 Sep 20
geweke<-function(x, frac1 = 0.2, frac2 = 0.5){
xstart <- c(start(x), end(x) - frac2 * (end(x) - start(x)))[c(1,3)]
xend <- c(start(x) + frac1 * (end(x) - start(x)), end(x))[c(1,3)]
y.variance <- y.mean <- vector("list", 2)
for (i in 1:2) {
y <- window(x, start = xstart[i], end = xend[i])
y.mean[[i]] <- apply(as.matrix(y), 2, mean)
y.variance[[i]] <- apply(as.matrix(y), 2, var)
}
z <- (y.mean[[1]] - y.mean[[2]])/sqrt(y.variance[[1]] + y.va…text/html2016-01-24T14:48:26+00:00Anonymous (anonymous@undisclosed.example.com)r_function_for_processing_dic
https://nd.psychstat.org/lab/projects/02robust_growth_curve_modeling_using_student-t_distribution/r_function_for_processing_dic
proc.coda<-function(stem,N, I){
res<-NULL
for (i in 1001:2000){
codafile<-paste('../', stem,I,"/",stem,"CODAN",N,"I",i,"dic.txt",sep='')
file.res<-paste(stem,'-', I,'-',N,'DIC.txt', sep='')
if (file.exists(codafile)){
DIC<-suppressWarnings(readLines(codafile))
DIC<-DIC[substr(DIC,1,1)=="y"]
cat(c(DIC, i), file=file.res, append=T)
cat('\n', file=file.res, append=T)
}
}
}
proc.coda('Tdata',100,'Linear')
proc.coda('Tdata',200,'Linear')
proc.coda(…text/html2016-01-24T14:48:26+00:00Anonymous (anonymous@undisclosed.example.com)real_data_analysis_using_auxiliary_variables
https://nd.psychstat.org/lab/projects/02robust_growth_curve_modeling_using_student-t_distribution/real_data_analysis_using_auxiliary_variables
#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[i])
muY[i,t]<-LS[i,1]+LS[i,2]*t
}
w[i]~dgamma(a, a)
Inv_Sig_e2[i]<-w[i]*Pre_e2
}
#Priors for model parameter
#k~dexp(0.1)
k~dunif(0,100)
a<-k/2
for (…text/html2016-01-24T14:48:26+00:00Anonymous (anonymous@undisclosed.example.com)sas_codes
https://nd.psychstat.org/lab/projects/02robust_growth_curve_modeling_using_student-t_distribution/sas_codes
*copy data to the running directory;
%MACRO runsim(i);
data temp;
number = &i;
datastr=put(number, 4.);
datastrcomp=compress(datastr,'');
commd ='cp /pscratch/zzhang4/Private/student-grm/data/folder/DATAN500I'||datastrcomp||'bugs.txt data.txt';
run;
data _null_;
set temp;
file 'cp.sh';
put commd;
run;
data _null_;
X 'chmod 755 cp.sh';
X './cp.sh';
X 'rm -f cp.sh';
run;
quit;
*copy coda to the results locations;
data temp1;
r1 = &i;
rstr=put(r1, 4.);
codar=compress(rstr,'');
commd1 =…text/html2016-01-24T14:48:26+00:00Anonymous (anonymous@undisclosed.example.com)simulate_data_and_analyze_in_mplus
https://nd.psychstat.org/lab/projects/02robust_growth_curve_modeling_using_student-t_distribution/simulate_data_and_analyze_in_mplus
## Simulated t data
y<-read.table('data-generated-normal-coda.txt')
ymat<-array(y[, 2], dim=c(1000,5, 500))
t(ymat[2,,])->temp
plot(1:5, temp[1,1:5], ylim=c(30, 150), type='l')
for (i in 2:200){
lines(1:5, temp[i,1:5])
}
## write data into txt files
for (i in 1001:2000){
t(ymat[i-1000,,])->temp
filename<-paste('./Tdata/TdataN500I',i,'.txt', sep='')
write.table(temp, filename, row.names=F, col.names=F)
}
source("http://nd.psychstat.org/_media/classdata/rtobugs.r")
for (i in 1001:2000)…text/html2016-01-24T14:48:26+00:00Anonymous (anonymous@undisclosed.example.com)tabulated_results
https://nd.psychstat.org/lab/projects/02robust_growth_curve_modeling_using_student-t_distribution/tabulated_results
Parameters include
* $\beta_L=60$
* $\beta_S=3$
* $\sigma_L^2=200$
* $\sigma_S^2=3$
* $\sigma_{LS}=0$
* $\sigma_e^2=30$
* $k=3$
T data using T model
point est bias s.e. cvg1 power1 cvg2 power2 $\beta_L=60$ 59.95971 -0.000671522 0.7178075 0.953 1 0.951 1 $\beta_S=3$ 2.993026 -0.002324607 0.1281185 0.952 1 0.951 1 $\sigma^{2}_L=200$ 198.5271 $\sigma^{2}_s=3$$\sigma_{LS}=0$$\sigma^{2}_e=30$$k=3$$\beta_L=60$$\beta_S=3$$\sigma^{2}_L=200$$\sigma^{2}_s=3$$\si…text/html2016-01-24T14:48:26+00:00Anonymous (anonymous@undisclosed.example.com)using_the_uniform_prior_for_real_data
https://nd.psychstat.org/lab/projects/02robust_growth_curve_modeling_using_student-t_distribution/using_the_uniform_prior_for_real_data
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, k)
muY[i,t]<-LS[i,1]+LS[i,2]*t
}
}
k~dunif(2, 30)
#Priors for model parameter
for (i in 1:1){
bL[i] ~ dnorm(0, 1.0E-6)
bS[i] ~ dnorm(0, 1.0E-6)
}
Inv_cov[1:2,1:2]~dwish(R[1:2,1:2], 2)
R[1,1]<-1
R[2,2]<-1
R[2,1]<-R[1,2]
R[1,2]<…