# n: (sub)sample size 
# cp: vector of class probabilities
# id: indices of the (sub)sample 
# n.draws: the number of pseudo-class draws
# sim.class: the matrix of pseudo-class draw results; 
## the jth column corresponds to the jth pseudo-class draw
# missum: the vactor of the number of missing data for each subject
# df: the vector of number of latent variables in each class

################################################################################
### Q-Q norm plot of pseudo-class adjusted standardized residuals at level 2 ###
################################# begin ########################################
for ( k in 1:n.class)
	{
 for ( u in ind.r[[k]])
   {
	 len_0
	 for(j in 1:n.draws)
	 {
	idin_id[sim.class[,j]==k]
	len_len+length(idin)
    }
    q_1:(len/n.draws)/(len/n.draws+1)
    qx_matrix(0,n.draws,length(q)) 
    for(j in 1:n.draws)
	{
	idin_id[sim.class[,j]==k]
    qx[j,]_ quantile(sd.res2[[k]][idin,u],q,na.rm=T)
    }

plot(qnorm(q), apply(qx,2,mean),ylim=c(-3.5,3.5),xlim=c(-3.5,3.5),ylab=paste("latent",u))
abline(0,1)
qxm_apply(qx,2,mean)

# add the 95% confidence intervals that do not cover the expected quantiles  
for( v in 1: length(q))
   {
	 if( (qnorm(q[v])>(qxm[v]+1.96*sqrt(q[v]*(1-q[v]))/sqrt(n*cp[k])/dnorm(qnorm(q[v]))))|(qnorm(q[v])<(qxm[v]-1.96*sqrt(q[v]*(1-q[v]))/sqrt(n*cp[k])/dnorm(qnorm(q[v])))) )
    segments(qnorm(q[v]),qxm[v]+1.96*sqrt(q[v]*(1-q[v]))/sqrt(n*cp[k])/dnorm(qnorm(q[v])),qnorm(q[v]),qxm[v]-1.96*sqrt(q[v]*(1-q[v]))/sqrt(n*cp[k])/dnorm(qnorm(q[v])))
    }
title(paste("LTC ",k))
   }
  }
#######################################end #######################################

################################################################################
### Q-Q plot of paired class standardized Ma residuals at level 1 ##############
################################# begin ########################################
for ( k1 in 1:(n.class-1))
	 {
	for ( k2 in (k1+1):n.class)
		{
	    len1_len2_0
       for(j in 1:n.draws)
	       { 
	        idin1_id[sim.class[,1]==k1]
	        idin2_id[sim.class[,1]==k2]
           len1_len1+length(idin1) 
           len2_len2+length(idin2)
	       }
	    qx_qy_matrix(0,n.draws,floor(len1/n.draws))
    	 q_1:floor(len1/n.draws)/(floor(len1/n.draws)+1)
	  
	    for(j in 1:n.draws)
	       { 
	        idin1_id[sim.class[,1]==k1]
	        idin2_id[sim.class[,1]==k2]
           qx[j,]_quantile(wils.hilf(ma.res1[[k1]][idin1],9-missum[idin1]),q,na.rm=T)
           qy[j,]_quantile(wils.hilf(ma.res1[[k2]][idin2],9-missum[idin2]),q,na.rm=T)
          }
       plot(apply(qx,2,mean),apply(qy,2,mean),ylab=paste("LTC",k2),xlab=paste("LTC",k1),xlim=c(-3,3),ylim=c(-3,3))
       abline(0,1,lty=2)

       qxm_apply(qx,2,mean)
       qym_apply(qy,2,mean)

      # add the 95% confidence intervals that do not cover the expected quantiles  
       for ( v in 1:length(q))
           {  	
           if ((qym[v]>(qxm[v]+1.96*sqrt(q[v]*(1-q[v]))/sqrt(n*cp[k2])/dnorm(qnorm(q[v]))))|( qym[v]<(qxm[v]-1.96*sqrt(q[v]*(1-q[v]))/sqrt(n*cp[k2])/dnorm(qnorm(q[v])))) )    
           segments(qxm[v],qym[v]+1.96*sqrt(q[v]*(1-q[v]))/sqrt(n*cp[k2])/dnorm(qnorm(q[v])),qxm[v],qym[v]-1.96*sqrt(q[v]*(1-q[v]))/sqrt(n*cp[k2])/dnorm(qnorm(q[v])))
           }
       for ( v in 1:length(q))
           {  
            if ((qxm[v]>(qym[v]+1.96*sqrt(q[v]*(1-q[v]))/sqrt(n*cp[k1])/dnorm(qnorm(q[v]))))|( qxm[v]<(qym[v]-1.96*sqrt(q[v]*(1-q[v]))/sqrt(n*cp[k1])/dnorm(qnorm(q[v])))) )
            segments(qxm[v]+1.96*sqrt(q[v]*(1-q[v]))/sqrt(n*cp[k1])/dnorm(qnorm(q[v])),qym[v],qxm[v]-1.96*sqrt(q[v]*(1-q[v]))/sqrt(n*cp[k1])/dnorm(qnorm(q[v])),qym[v])
           }

     }
   }
###################################### End ######################################


################################################################################
### Q-Q plot of paired class standardized Ma residuals at level 2 ##############
################################# begin ########################################

def_c(2,2,1,1)
for ( k1 in 1: (n.class-1)){
for( k2 in (k1+1):n.class){
    len1_len2_0
    for(j in 1:n.draws)
	    { 
	     idin1_id[sim.class[,1]==k1]
	     idin2_id[sim.class[,1]==k2]
        len1_len1+length(idin1) 
        len2_len2+length(idin2)
	    }
	 qx_qy_matrix(0,n.draws,floor(len1/n.draws))
	 q_1:floor(len1/n.draws)/(floor(len1/n.draws)+1)
	  
	 for(j in 1:n.draws)
	    { 
	     idin1_id[sim.class[,1]==k1]
	     idin2_id[sim.class[,1]==k2]
        qx[j,]_quantile(wils.hilf(ma.res2[[k1]][idin1],def[k1]),q,na.rm=T)
        qy[j,]_quantile(wils.hilf(ma.res2[[k2]][idin2],def[k2]),q,na.rm=T)
        }
     plot(apply(qx,2,mean),apply(qy,2,mean),ylab=paste("LTC",k2),xlab=paste("LTC",k1),xlim=c(-3,3),ylim=c(-3,3))
     abline(0,1)

     qxm_apply(qx,2,mean)
     qym_apply(qy,2,mean)

     # add the 95% confidence intervals that do not cover the expected quantiles  
     for( v in 1:length(q))
        {  	
        if( (qym[v]>(qxm[v]+1.96*sqrt(q[v]*(1-q[v]))/sqrt(n*cp[k2])/dnorm(qnorm(q[v]))))|( qym[v]<(qxm[v]-1.96*sqrt(q[v]*(1-q[v]))/sqrt(n*cp[k2])/dnorm(qnorm(q[v])))) )
          segments(qxm[v],qym[v]+1.96*sqrt(q[v]*(1-q[v]))/sqrt(n*cp[k2])/dnorm(qnorm(q[v])),qxm[v],qym[v]-1.96*sqrt(q[v]*(1-q[v]))/sqrt(n*cp[k2])/dnorm(qnorm(q[v])))
        }

     for( v in 1:length(q))
        {  
	      if( (qxm[v]>(qym[v]+1.96*sqrt(q[v]*(1-q[v]))/sqrt(n*cp[k1])/dnorm(qnorm(q[v]))))|( qxm[v]<(qym[v]-1.96*sqrt(q[v]*(1-q[v]))/sqrt(n*cp[k1])/dnorm(qnorm(q[v])))) )
         segments(qxm[v]+1.96*sqrt(q[v]*(1-q[v]))/sqrt(n*cp[k1])/dnorm(qnorm(q[v])),qym[v],qxm[v]-1.96*sqrt(q[v]*(1-q[v]))/sqrt(n*cp[k1])/dnorm(qnorm(q[v])),qym[v])
        }
     }
   }
###################################### End ######################################