npmmApproximation <- function(Hmin,Hmax,Kmin,Kmax,nL1vars,nL2vars,resvar_con=FALSE,Mplusoutput=NA,plotoutput=FALSE,
intslopes=NA, resvar=NA, mcwi=NA, mcws=NA, mcbi=NA)
{
approx.ubertable.fixed <- array(NA,c(1+nL1vars+nL2vars,1,as.numeric(paste(c(Hmax,Kmax),collapse=""))))
approx.ubertable.tau <- array(NA,c(1+nL1vars+nL2vars,1+nL1vars+nL2vars,as.numeric(paste(c(Hmax,Kmax),collapse=""))))
approx.ubertable.sigma <- array(NA,c(1+nL1vars+nL2vars,1+nL1vars+nL2vars,as.numeric(paste(c(Hmax,Kmax),collapse=""))))
fixedeffects.output <- list()
tau.output <-  list()
sigma.output <-  list()
for (z in Hmin:Hmax)
{
for (v in Kmin:Kmax)
{
##specify H and K
H <- z
K <- v
if (is.na(Mplusoutput)==FALSE)
{
mplusout <- NA
intslopes <- NA
resvar <- NA
mcwi <- NA
mcws <- NA
mcbi <- NA
##read in Mplus output
#total length of file
length <- H*K*(nL1vars+1)+H*K*(nL2vars+1)+H-1+K-1+(H-1)*(K-1)
if (resvar_con==TRUE)
{
length <- H*K*nL1vars+1+H*K*(nL2vars+1)+H-1+K-1+(H-1)*(K-1)
}
mplusout <- NA
intslopes <- NA
resvar <- NA
mcwi <- NA
mcws <- NA
mcbi <- NA
##read in Mplus output data
try(
Mplusout <- matrix(data=scan(paste(Mplusoutput,H,"by",K,".dat",sep="")),length,1)
)
if (is.na(Mplusout)==TRUE) next
##intercepts
intercepts <- Mplusout[H*K*(nL1vars+1)+1,]
if (H*K > 1)
{
for (i in 1:c(K*H-1))
{
intercepts <- c(intercepts,Mplusout[c(H*K*(nL1vars+1)+1+i*(nL2vars+1)),])
}
}
if (resvar_con==TRUE)
{
intercepts <- Mplusout[H*K*nL1vars+2,]
if (H*K > 1)
{
for (i in 1:c(K*H-1))
{
intercepts <- c(intercepts,Mplusout[c(H*K*nL1vars+2+i*(nL2vars+1)),])
}
}
}
##xslopes
xslopes <- Mplusout[1:nL1vars,]
if (H*K > 1)
{
for (i in 0:c(K*H-2))
{
xslopes <- c(xslopes,Mplusout[c(nL1vars+2+i*(nL1vars+1)):c(nL1vars+2+i*(nL1vars+1)+nL1vars-1),])
}
}
if (resvar_con==TRUE)
{
xslopes <- Mplusout[1:nL1vars,]
if (H*K > 1)
{
xslopes <- c(xslopes,Mplusout[c(nL1vars+2):c(H*K*nL1vars+1),])
}
}
##reorder
xslopes_temp <- list()
for (i in 1:nL1vars)
{
xslopes_temp[[i]] <- xslopes[seq(i, length(xslopes), nL1vars)]
if (i==1)
{
xslopes_temp2 <- xslopes_temp[[i]]
}
if (i>1)
{
xslopes_temp2 <- c(xslopes_temp2,xslopes_temp[[i]])
}
}
xslopes <- xslopes_temp2
##wslopes
wslopes <- Mplusout[c(H*K*(nL1vars+1)+2):c(H*K*(nL1vars+1)+2+nL2vars-1),]
if (H*K > 1)
{
for (i in 1:c(K*H-1))
{
wslopes <- c(wslopes,Mplusout[c(H*K*(nL1vars+1)+2+i*(nL2vars+1)):c(H*K*(nL1vars+1)+2+i*(nL2vars+1)+nL2vars-1),])
}
}
if (resvar_con==TRUE)
{
wslopes <- Mplusout[c(H*K*nL1vars+3):c(H*K*nL1vars+3+nL2vars-1),]
if (H*K > 1)
{
for (i in 1:c(K*H-1))
{
wslopes <- c(wslopes,Mplusout[c(H*K*nL1vars+3+i*(nL2vars+1)):c(H*K*nL1vars+3+i*(nL2vars+1)+nL2vars-1),])
}
}
}
##reorder
wslopes_temp <- list()
for (i in 1:nL1vars)
{
wslopes_temp[[i]] <- wslopes[seq(i, length(wslopes), nL1vars)]
if (i==1)
{
wslopes_temp2 <- wslopes_temp[[i]]
}
if (i>1)
{
wslopes_temp2 <- c(wslopes_temp2,wslopes_temp[[i]])
}
}
wslopes <- wslopes_temp2
##resvar
resvar=Mplusout[nL1vars+1,]
if (resvar_con==FALSE)
{
if (H*K > 1)
{
for (i in 1:c(K*H-1))
{
resvar <- c(resvar,Mplusout[c(nL1vars+1+i*(nL1vars+1)),])
}
}
}
##full ints and slopes
intslopes=c(intercepts,xslopes,wslopes)
##mcw
if (K>1)
{
mcwi=c(Mplusout[c(H*K*(nL1vars+1)+H*K*(nL2vars+1)+1+H-1):c(H*K*(nL1vars+1)+H*K*(nL2vars+1)+1+H-1+K-2)],0)
}
if (K==1)
{
mcwi=0
}
if (resvar_con==TRUE)
{
if (K>1)
{
mcwi=c(Mplusout[c(H*K*nL1vars+1+H*K*(nL2vars+1)+1+H-2+1):c(H*K*nL1vars+1+H*K*(nL2vars+1)+1+H-2+1+K-2)],0)
}
if (K==1)
{
mcwi=0
}
}
##ms
if (H>1 & K>1)
{
mcws_temp=Mplusout[c(H*K*(nL1vars+1)+H*K*(nL2vars+1)+1+H-1+K-1):c(H*K*(nL1vars+1)+H*K*(nL2vars+1)+1+H-1+K-1+(H-1)*(K-1)-1)]
mcws <- list()
for (i in 1:H)
{
mcws <- c(mcws,mcws_temp[c(1+(i-1)*(K-1)):c(i+i*(K-2))],0)
}
mcws[c((H-1)*K+1):c(H*K)] <- 0
}
if (H==1 | K==1)
{
mcws=matrix(0,H*K,1)
}
mcws <- as.numeric(mcws)
if (resvar_con==TRUE)
{
if (H>1 & K>1)
{
mcws_temp=Mplusout[c(H*K*nL1vars+1+H*K*(nL2vars+1)+1+H-2+1+K-2+1):c(H*K*nL1vars+1+H*K*(nL2vars+1)+1+H-2+1+K-2+1+(H-1)*(K-1)-1)]
mcws <- list()
for (i in 1:H)
{
mcws <- c(mcws,mcws_temp[c(1+(i-1)*(K-1)):c(i+i*(K-2))],0)
}
mcws[c((H-1)*K+1):c(H*K)] <- 0
}
if (H==1 | K==1)
{
mcws=matrix(0,H*K,1)
}
mcws <- as.numeric(mcws)
}
##mcb
if (H>1)
{
mcbi=c(Mplusout[c(H*K*(nL1vars+1)+H*K*(nL2vars+1)+1):c(H*K*(nL1vars+1)+H*K*(nL2vars+1)+1+H-2)],0)
}
if (H==1)
{
mcbi=0
}
if (resvar_con==TRUE)
{
if (H>1)
{
mcbi=c(Mplusout[c(H*K*nL1vars+1+H*K*(nL2vars+1)+1):c(H*K*nL1vars+1+H*K*(nL2vars+1)+1+H-2)],0)
}
if (H==1)
{
mcbi=0
}
}
}
##read in intercepts and slopes
int_slopes <- aperm(array(data=intslopes,dim=c(K,H,1+nL1vars+nL2vars)),perm=c(2,1,3))
##read in residual variances
residual_var <- matrix(data=resvar,H,K,byrow=TRUE)
##read in multinomials
#L1 class intercepts
mcw <- matrix(data=mcwi,1,K)
mcw[1,K] <- 0
#L2 class on L1 class slopes
ms <- matrix(data=mcws,H,K,byrow=TRUE)
ms[H,1:K] <- 0
ms[1:H,K] <- 0
#L2 class intercepts
mcb <- matrix(data=mcbi,H,1)
mcb[H,1] <- 0
##compute probabilities of class membership
#denominator for p_cbcw for each cbcw
den_cbcw <- matrix(NA,H,K)
for (i in 1:H)
{
for (j in 1:K)
{
den_cbcw[i,j] <- sum(exp(mcw[1,1:K]+ms[i,1:K]))
}
}
#prob of cw given cb
prob_cwgivencb <- matrix(NA,H,K)
for (i in 1:H)
{
for (j in 1:K)
{
prob_cwgivencb[i,j] <- exp(mcw[1,j]+ms[i,j]) / den_cbcw[i,j]
}
}
#marginal prob of cb
prob_cb <- matrix(NA,H,1)
for (i in 1:H)
{
prob_cb[i,1] <- exp(mcb[i,1])/ sum(exp(mcb[1:H,1]))
}
#class combination prob
prob_cbcw <- matrix(NA,H,K)
for (i in 1:H)
{
for (j in 1:K)
{
prob_cbcw[i,j] <- prob_cb[i,1]*prob_cwgivencb[i,j]
}
}
##compute marginal L2 class parameters
margL2params <- array(3,c(H,1,1+nL1vars+nL2vars))
for (i in 1:H)
{
for (j in 1:c(1+nL1vars+nL2vars))
{
margL2params[i,1,j] <- sum(prob_cwgivencb[i,1:K]*int_slopes[i,1:K,j])
}
}
##compute fixed effects
fixedeffects <- matrix(NA,1+nL1vars+nL2vars,1)
for (i in 1:c(1+nL1vars+nL2vars))
{
fixedeffects[i,1] <- sum(prob_cb*margL2params[,,i])
}
##compute implied tau matrix
impliedtau <- matrix(NA,1+nL1vars+nL2vars,1+nL1vars+nL2vars)
for (i in 1:c(1+nL1vars+nL2vars))
{
for (j in 1:c(1+nL1vars+nL2vars))
{
impliedtau[i,j] <- sum(prob_cb*margL2params[,1,i]*margL2params[,1,j])-sum(prob_cb*margL2params[,1,i])*sum(prob_cb*margL2params[,1,j])
}
}
##implied sigma equation
sigmamatrix <- matrix(NA,1+nL1vars+nL2vars,1+nL1vars+nL2vars)
for (i in 1:c(1+nL1vars+nL2vars))
{
for (j in 1:c(1+nL1vars+nL2vars))
{
sigmamatrix[i,j] <- sum(prob_cbcw*((int_slopes[,,i]-margL2params[,,i])*
(int_slopes[,,j]-margL2params[,,j])))
}
}
sigmamatrix[1,1] <- sigmamatrix[1,1]+sum(prob_cbcw*residual_var)
approx.ubertable.fixed[,,as.numeric(paste(c(z,v),collapse=""))] <-
fixedeffects
approx.ubertable.tau[,,as.numeric(paste(c(z,v),collapse=""))] <-
impliedtau
approx.ubertable.sigma[,,as.numeric(paste(c(z,v),collapse=""))] <-
sigmamatrix
#fixedeffects.output[[as.numeric(paste(c(z,v),collapse=""))]] <- rbind2(paste("# L1 classes = ",K,"; # L2 classes = ",H),fixedeffects)
#tau.output[[as.numeric(paste(c(z,v),collapse=""))]] <- list(paste("# L1 classes = ",K,"; # L2 classes = ",H),impliedtau)
#sigma.output[[as.numeric(paste(c(z,v),collapse=""))]] <- list(paste("# L1 classes = ",K,"; # L2 classes = ",H),sigmamatrix)
#names(fixedeffects.output[[as.numeric(paste(c(z,v),collapse=""))]]) <- paste("# L1 classes = ",K,"; # L2 classes = ",H)
fixedeffects.output[[as.numeric(paste(c(z,v),collapse=""))]] <- matrix(c(paste("K = ",K,"; H = ",H),fixedeffects),2+nL1vars+nL2vars,1)
rownames(fixedeffects.output[[as.numeric(paste(c(z,v),collapse=""))]]) <- c("",seq(from=1,to=1+nL1vars+nL2vars,by=1))
colnames(fixedeffects.output[[as.numeric(paste(c(z,v),collapse=""))]]) <- ""
blank <- replicate(nL1vars+nL2vars,"")
headers <- seq(from=1,to=1+nL1vars+nL2vars,by=1)
tau.output[[as.numeric(paste(c(z,v),collapse=""))]] <- matrix(c(paste("K = ",K,"; H = ",H),blank,headers,impliedtau),3+nL1vars+nL2vars,1+nL1vars+nL2vars,byrow=TRUE)
rownames(tau.output[[as.numeric(paste(c(z,v),collapse=""))]]) <- c("","",seq(from=1,to=1+nL1vars+nL2vars,by=1))
colnames(tau.output[[as.numeric(paste(c(z,v),collapse=""))]]) <- replicate(1+nL1vars+nL2vars,"")
sigma.output[[as.numeric(paste(c(z,v),collapse=""))]] <- matrix(c(paste("K = ",K,"; H = ",H),blank,headers,sigmamatrix),3+nL1vars+nL2vars,1+nL1vars+nL2vars,byrow=TRUE)
rownames(sigma.output[[as.numeric(paste(c(z,v),collapse=""))]]) <- c("","",seq(from=1,to=1+nL1vars+nL2vars,by=1))
colnames(sigma.output[[as.numeric(paste(c(z,v),collapse=""))]]) <- replicate(1+nL1vars+nL2vars,"")
}
}
##plots
if (plotoutput != FALSE)
{
pdf(file = paste(plotoutput,".pdf"))
classlist <- as.numeric(paste(c(Hmin,Kmin),collapse=""))
if (Hmax > Hmin)
{
for (j in c(Hmin+1):Hmax)
{
classlist <- c(classlist,as.numeric(paste(c(j,Kmin),collapse="")))
}
}
##fixed effects plots
for (i in 1:c(1+nL1vars+nL2vars))
{
plot(y=as.numeric(approx.ubertable.fixed[i,,c(classlist)]),x=c(Hmin:Hmax),type="b",xlim=c(Hmin,Hmax),ylim=c(min(as.numeric(approx.ubertable.fixed[i,,]),na.rm=TRUE),max(as.numeric(approx.ubertable.fixed[i,,]),na.rm=TRUE)),lwd=3,col=Kmin,xlab=expression(italic("H")),ylab=paste("Implied fixed effect ",i))
classlist_temp <- classlist
if (Kmax > Kmin)
{
for (l in c(Kmin+1):Kmax)
{
classlist_temp <- classlist_temp+1
lines(y=as.numeric(approx.ubertable.fixed[i,,classlist_temp]),x=c(Hmin:Hmax),type="b",col=l,lty=1,lwd=3)
}
legend("topright", legend = paste("K =",Kmin:Kmax), col=c(Kmin:Kmax), lty=c(1:1),lwd=3,seg.len=4)
}
}
##tau plots
for (i in 1:c(1+nL1vars+nL2vars))
{
for (j in 1:c(1+nL1vars+nL2vars))
{
plot(y=as.numeric(approx.ubertable.tau[i,j,c(classlist)]),x=c(Hmin:Hmax),type="b",ylim=c(min(as.numeric(approx.ubertable.tau[i,j,]),na.rm=TRUE),max(as.numeric(approx.ubertable.tau[i,j,]),na.rm=TRUE)),xlim=c(Hmin,Hmax),lwd=3,col=Kmin,xlab=expression(italic("H")),ylab=paste("Implied tau ",i,j))
classlist_temp <- classlist
if (Kmax > Kmin)
{
for (l in c(Kmin+1):Kmax)
{
classlist_temp <- classlist_temp+1
lines(y=as.numeric(approx.ubertable.tau[i,j,classlist_temp]),x=c(Hmin:Hmax),type="b",col=l,lty=1,lwd=3)
}
legend("topright", legend = paste("K =",Kmin:Kmax), col=c(Kmin:Kmax), lty=c(1:1),lwd=3,seg.len=4)
}
}
}
##sigma plots
classlist_sigma <- as.numeric(paste(c(Hmin,Kmin),collapse=""))
if (Kmax > Kmin)
{
for (j in c(Kmin+1):Kmax)
{
classlist_sigma <- c(classlist_sigma,as.numeric(paste(c(Hmin,j),collapse="")))
}
}
for (i in 1:c(1+nL1vars+nL2vars))
{
for (j in 1:c(1+nL1vars+nL2vars))
{
plot(y=as.numeric(approx.ubertable.sigma[i,j,c(classlist_sigma)]),x=c(Kmin:Kmax),type="b",ylim=c(min(as.numeric(approx.ubertable.sigma[i,j,]),na.rm=TRUE),max(as.numeric(approx.ubertable.sigma[i,j,]),na.rm=TRUE)),xlim=c(Kmin,Kmax),lwd=3,col=Kmin,xlab=expression(italic("K")),ylab=paste("Implied sigma ",i,j))
classlist_sigma_temp <- classlist_sigma
if (Hmax > Hmin)
{
for (l in c(Hmin+1):Hmax)
{
classlist_sigma_temp <- classlist_sigma_temp+10
lines(y=as.numeric(approx.ubertable.sigma[i,j,classlist_sigma_temp]),x=c(Kmin:Kmax),type="b",col=l,lty=1,lwd=3)
}
legend("topright", legend = paste("H =",Hmin:Hmax), col=c(Hmin:Hmax), lty=c(1:1),lwd=3,seg.len=4)
}
}
}
dev.off()
}
##Output
fixedeffects.output <- Filter(Negate(is.null), fixedeffects.output)
tau.output <- Filter(Negate(is.null), tau.output)
sigma.output <- Filter(Negate(is.null), sigma.output)
Output <- c(noquote(fixedeffects.output),noquote(tau.output),noquote(sigma.output))
names(Output) <- c(replicate(length(fixedeffects.output),'Implied Fixed 
Effects'),replicate(length(fixedeffects.output),"Implied Tau Matrix"),replicate(length(tau.output),"Implied Sigma Matrix"))
return(Output)
}