
powzerodiff <- 
#sample size for group difference with null equal to exact diff
#dfa =      Degrees of freedom for model A
#dfb =      Degrees of freedom for model B
#alpha =    Significance level
#power =    Desired power
#rmseaa =   RMSEA for model A under H1
#rmseab =   RMSEA for model B under H1
#g =        Number of groups
function(dfa, dfb, alpha, power, rmseaa, rmseab, g)
{
ddiff <- dfa-dfb                        
fa <- (dfa*rmseaa**2)/sqrt(g)           
fb <- (dfb*rmseab**2)/sqrt(g)           
#initialize values
powa <- 0.0
n <- 0
#begin loop for finding initial level of n
while(powa < power) {
  n <- n + 100                          
  ncp <- (n-1)*(fa-fb)                  
  cval <- qchisq(1-alpha,df=ddiff,ncp=0)
  powa <- 1-pchisq(cval,df=ddiff,ncp=ncp)
}
#begin loop for interval halving
dir <- -1
newn <- n
intv <- 200
powdiff <- powa-power
while(powdiff > .001) {
  intv <- intv*.5
  newn <- newn + dir*intv*.5
  #compute new power
  ncp <- (newn-1)*(fa-fb)
  #compute power
  cval <- qchisq(1-alpha,df=ddiff,ncp=0)
  powa <- 1-pchisq(cval,df=ddiff,ncp=ncp)
  powdiff <- abs(powa-power)
  if (powa < power) dir <- 1 else dir <- -1
}
minn <- round(newn)
return(minn)
}


powsmalldiff <- 
#sample size for group difference with null equal to small diff
#dfa =      Degrees of freedom for model A
#dfb =      Degrees of freedom for model B
#alpha =    Significance level
#power =    Desired power
#rmsea0a =  RMSEA for model A (constrained) under H0
#rmsea1a =  RMSEA for model A (unconstrained) under H1
#rmsea0b =  RMSEA for model B (constrained) under H0
#rmsea1b =  RMSEA for model B (unconstrained) under H1
#g =        Number of groups
function(dfa, dfb, alpha, power, rmsea0a, rmsea1a, rmsea0b, rmsea1b, g)
{
ddiff <- dfa-dfb                         
f0a <- (dfa*rmsea0a**2)/sqrt(g)          
f0b <- (dfb*rmsea0b**2)/sqrt(g)          
f1a <- (dfa*rmsea1a**2)/sqrt(g)          
f1b <- (dfb*rmsea1b**2)/sqrt(g)          
#initialize values
powa <- 0.0
n <- 0
#begin loop for finding initial level of n
while(powa < power) {
  n <- n + 100                             
  ncp0 <- (n-1)*(f0a-f0b)                  
  ncp1 <- (n-1)*(f1a-f1b)                  
  cval <- qchisq(1-alpha,df=ddiff,ncp=ncp0)
  powa <- 1-pchisq(cval,df=ddiff,ncp=ncp1)
}
#begin loop for interval halving
dir <- -1
newn <- n
intv <- 200
powdiff <- powa-power
while(powdiff > .001) {
  intv <- intv*.5
  newn <- newn + dir*intv*.5
  #compute new power
  ncp0 <- (newn-1)*(f0a-f0b)                   
  ncp1 <- (newn-1)*(f1a-f1b)                   
  #compute power
  cval <- qchisq(1-alpha,df=ddiff,ncp=ncp0)
  powa <- 1-pchisq(cval,df=ddiff,ncp=ncp1)
  powdiff <- abs(powa-power)
  if (powa < power) dir <- 1 else dir <- -1
}
minn <- round(newn)
return(minn)
}



ci<-function (mu, Sigma, quant, alpha = 0.05, type = "MC", 
    plot = FALSE, plotCI = FALSE, n.mc = 1e+06, H0 = FALSE, mu0 = NULL, 
    Sigma0 = NULL, ...) 
{
    if (missing(mu) | is.null(mu)) 
        stop(paste("argument", sQuote("mu"), "must be specified"))
    if (missing(quant) | is.null(quant)) 
        stop(paste("argument", sQuote("quant"), "must be specified"))
    if (!type %in% c("MC", "mc", "asymp", "Asymp", 
        "all")) 
        stop(paste("Please enter a ccorrect", sQuote("type"), 
            "argument"))
    if (class(mu) != "lavaan") {
        if (is.null(Sigma) | missing(Sigma)) 
            stop(paste("argument", sQuote("Sigma"), 
                "cannot be a NULL value"))
        if (!is.matrix(Sigma)) {
            if (length(mu) != (sqrt(1 + 8 * length(Sigma)) - 
                1)/2) 
                stop(paste("Please check the length of", 
                  sQuote("Sigma"), "and", sQuote("mu"), 
                  ". If the length(dimension) of the", 
                  sQuote("mu"), "vector (", length(mu), 
                  ") is correct, the stacked lower triangle matrix", 
                  sQuote("Sigma"), "must have ", 
                  ((2 * length(mu) + 1)^2 - 1)/8, "elements, instead of", 
                  length(Sigma)))
            Sigma <- lav_matrix_vech_reverse(Sigma)
        }
        if (is.null(names(mu))) 
            names(mu) <- paste("b", 1:length(mu), sep = "")
        if (!all(all.vars(quant) %in% names(mu))) 
            stop(paste("The parameters names in formula", 
                sQuote("quant"), "must match the parameters names provided in", 
                sQuote("mu"), "."))
    }
    else {
        fm1 <- mu
        if (!fm1@Options$do.fit) 
            fm1 <- update(mu, do.fit = TRUE)
        pEstM1 <- coef(fm1)
        name1 <- all.vars(quant)
        if (!all(name1 %in% names(pEstM1))) 
            stop(paste("The parameters names in formula", 
                sQuote("quant"), "must match the parameters names provided in lavaan object."))
        mu <- pEstM1[name1]
        covM1 <- (vcov(fm1))
        Sigma <- covM1[name1, name1]
    }
    if (length(mu) * n.mc > .Machine$integer.max) {
        n.mc <- 1e+06
        warning(paste("n.mc is too large. It is reset to", 
            n.mc))
    }
    if (plot) {
        op <- par(mar = c(7, 4, 7, 4) + 0.1, xpd = TRUE, ask = FALSE)
    }
    if (type %in% c("mc", "MC")) 
        res <- confintMC(mu = mu, Sigma = Sigma, quant = quant, 
            alpha = alpha, plot = plot, plotCI = plotCI, n.mc = n.mc, 
            H0 = H0, mu0 = mu0, Sigma0 = Sigma0, ...)
    if (type %in% c("asymp", "Asymp")) 
        res <- confintAsymp(mu = mu, Sigma = Sigma, quant = quant, 
            alpha = alpha, plot = plot, plotCI = plotCI)
    if (type %in% "all") {
        res1 <- confintMC(mu = mu, Sigma = Sigma, quant = quant, 
            type = type, alpha = alpha, plot = plot, plotCI = plotCI, 
            n.mc = n.mc)
        res2 <- confintAsymp(mu = mu, Sigma = Sigma, quant = quant, 
            type = type, alpha = alpha, plot = plot, plotCI = plotCI)
        res <- list(MC = res1, Asymptotic = res2)
    }
    if (plot) {
        par(op)
    }
    return(res)
}


confintMC <- function(mu, Sigma, quant=NULL, alpha=0.05, type="MC", plot=FALSE, plotCI=FALSE, n.mc = 1e+06, H0=FALSE, mu0, Sigma0, ...){
  q1 <- quant
  quant <- parse(text=sub("~","",quant))
  df <- data.frame(mvrnorm(n.mc,mu,Sigma))
  colnames(df) <-names(mu)
  quant.vec <- eval(quant,df) # MC Vector
  CI <- quantile(quant.vec,c(alpha/2,1-alpha/2))
  names(CI) <- c(paste((alpha/2*100),"%"),paste((1-alpha/2)*100,"%"))
  quantMean <- mean(quant.vec)
  quantSE <- sd(quant.vec)
  quantError <- quantSE/n.mc
  pMinH1 <- mean(quant.vec > 0) 
  pMinH1 <- 2*(min(pMinH1, 1-pMinH1))
  # This is to determine min and max of MC Samples--Do Not Delete this
 # if(H0) xrange <- c(min(-4*quantSE,quantMean-4*quantSE), max(4*quantSE ,quantMean+4*quantSE) )
 # else xrange <- c(quantMean-4*quantSE, quantMean+4*quantSE)
  xrange <- c(quantMean-4*quantSE, quantMean+4*quantSE)
  max1 <- xrange[2]
  min1 <- xrange[1]
  

  # Added for min null test 6-14-14
  # Calculates Acceptance Region
  if(H0)
  {
    if(missing(Sigma0)| is.null(Sigma0) ) Sigma0 <- Sigma
    if(!is.matrix(Sigma0)){
      if(length(mu)!= (sqrt(1 + 8 * length(Sigma0)) - 1)/2) stop(
        paste("Please check the length of", sQuote("Sigma0"),"and",sQuote("mu"),". If the length(dimension) of the", sQuote("mu"),"vector (",length(mu),") is correct, the stacked lower triangle matrix", sQuote("Sigma0"), "must have ",((2*length(mu)+1)^2-1)/8, "elements, instead of", length(Sigma0)) 
      )
      
      Sigma0 <- lav_matrix_vech_reverse(Sigma0) #converts to a symmetric matrix
    }
    
    #If mu0 is not specified, we use conservative min approach
    if(is.null(mu0) | missing(mu0) ){
      mu0 <- mu
      mu0s <- mu0/sqrt(diag(Sigma0)); #Srandardized
      mu0[which(mu0s==min(mu0s))] <- 0 # setting the smallest z value to 0
    }
    names(mu0) <- names(mu)
    df0 <- data.frame(mvrnorm(n.mc,mu0,Sigma0))
    colnames(df0) <-names(mu0)
    H0quant.vec <- eval(quant,df0)
    H0Mean <- mean(H0quant.vec)
    H0SE <- sd(H0quant.vec)
    H0CI <- quantile(H0quant.vec, c(alpha/2,1-alpha/2) )
    names(H0CI) <- c(paste((alpha/2*100),"%"),paste((1-alpha/2)*100,"%"))
    #yciH0<-par("usr")[3] 
    pMinH0 <- mean(H0quant.vec > quantMean) 
    pMinH0 <- 2*(min(pMinH0, 1-pMinH0))
    H0xy <- H0quant.vec + quantMean
    H0xy <- H0xy[H0xy>min1 & H0xy<max1] #This is used to make the plot prettier
#    H0xy <- H0xy+quantMean # We add mean to make it comparable with the CI
    H0xyDens <- density(H0xy)
    H0CI <- H0CI + quantMean # We add mean to make it comparable with the CI
    H0res <- list(CI=H0CI, Estimate= H0Mean,SE= H0SE,p= pMinH0, mu0) #Results
    #names(H0res) <- c( paste( (1-alpha)*100, "% ", "AC",sep="" ) , "Mean", "SE", "p", "mu0")   
    #attr(H0res,"quant")  <- q1
  }

  ###########   Plot ##########################
  
  
  if (plot){
    
    outer <- FALSE #outer position
    mcex <- .8

    if(type=="all")
    {
      res.asymp <- confintAsymp(mu=mu, Sigma=Sigma, quant=quant, type=type, alpha=alpha)
      range.asymp <- c(res.asymp$Estimate-5*res.asymp$SE, res.asymp$Estimate+5*res.asymp$SE)
      max1 <- max(max1,range.asymp)
      min1 <- min(min1,range.asymp)      
    }
    xy <- quant.vec[quant.vec>min1 & quant.vec<max1]
    xyDens <- density(xy)
    #smidge <- par("cin")*abs(par("tcl"))
    # Added for min null test 6-14-14
    # To get a more reasonable y range
    if(H0){ 
      yrange <- range(quantile(xyDens$y,c(alpha/10,1-alpha/10)),quantile(H0xyDens$y,c(alpha/10,1-alpha/10)))
    }
    else{
      yrange <- quantile(xyDens$y,c(alpha/10,1-alpha/10))
    }
    
    plot(xyDens,xlab="", ylab="", axes=FALSE, xlim=xrange, ylim=yrange, main="", lwd=2)
    mtext(quant,1,5)
    xrange<-pretty(xrange,n=9)
    axis(1,xrange,xrange, line=2.5);
    axis(2,line=1.1)
    
    # This adds legends, bars etc
    if(H0)
    {
      lines(H0xyDens, lty=2, col="blue", lwd=2) #Reference Dist
      #arrows(H0CI[1],yciH0,H0CI[2],yciH0,length=0,angle=90,code=3,cex=1.5,lwd=5,lty=1, col="blue")
      segments(H0CI[1],par("usr")[3],H0CI[1], par("usr")[4]/2, cex=1.5,lwd=2,lty=2, col="blue")
      segments(H0CI[2],par("usr")[3],H0CI[2], par("usr")[4]/2, cex=1.5,lwd=2,lty=2, col="blue")
      
      mtext(paste("P Value=",round(pMinH0,4) ), side=3, line=-1 , outer=outer, at=max1-1*(max1-min1)/9, cex=mcex, adj=0, col="blue")
      
      mtext(paste("Kurtosis=",round(kurtosis(H0quant.vec, type=2),3)), side=3, line=0, outer=outer, at=max1-1*(max1-min1)/9, cex=mcex, adj=0, col="blue")
      mtext(paste("Skewness=",round(skewness(H0quant.vec, type=2),3)), side=3, line=1, outer=outer, at=max1-1*(max1-min1)/9, cex=mcex, adj=0, col="blue")
      mtext(paste("Critical Value=", round(H0CI[1],3)), side=3,line=2, outer=outer, at=max1-1*(max1-min1)/9, cex=mcex, adj=0, col="blue")
      mtext(paste("Critical Value=", round(H0CI[2],3)), side=3,line=3, outer=outer, at=max1-1*(max1-min1)/9, cex=mcex, adj=0, col="blue")
      mtext(paste("H0:",quant,"=0"), side=3,line=4, outer=outer, at=max1-1*(max1-min1)/9, cex=mcex, adj=0, col="blue",font=2)
      text(H0CI[1], par("usr")[4]/2, labels="Lower Critical Value",adj=c(.5,0), cex=mcex )
      text(H0CI[2], par("usr")[4]/2, labels="Upper Critical Value", adj=c(.5,0), cex=mcex)
    }
    
    
    if(type %in% c("mc", "MC")){
      #New- 1/24/14-DT
      mtext(paste("P Value=",round(pMinH1,4)), side=3, line=-1, outer=outer, at=max1-3*(max1-min1)/9, cex=mcex,adj=0)
      mtext(paste("Kurtosis=",round(kurtosis(quant.vec, type=2),3)), side=3, line=0, outer=outer, at=max1-3*(max1-min1)/9, cex=mcex, adj=0)
      mtext(paste("Skewness=",round(skewness(quant.vec, type=2),3)), side=3, line=1, outer=outer, at=max1-3*(max1-min1)/9, cex=mcex, adj=0)      
      mtext(paste("LL=",round(CI[1],3)),side=3,line=2,outer=outer, at=max1-3*(max1-min1)/9, cex=mcex, adj=0)
      mtext(paste("UL=",round(CI[2],3)),side=3,line=3,outer=outer,at=max1-3*(max1-min1)/9, cex=mcex, adj=0)
      mtext("Monte Carlo CI",side=3,line=4,outer=outer,at=max1-3*(max1-min1)/9, cex=mcex, font=2, adj=0)
    }
    
    if(type=="all"){
      #New- 1/24/14-DT
      mtext(paste("Kurtosis=",round(kurtosis(quant.vec, type=2),3)), side=3,line=1,outer=outer, at=max1-3*(max1-min1)/9, cex=mcex, adj=0)
      mtext(paste("Skewness=",round(skewness(quant.vec, type=2) ,3)) ,side=3,line=2,outer=outer, at=max1-3*(max1-min1)/9, cex=mcex, adj=0)
      mtext(paste("LL=", round(CI[1],3)), side=3,line=3, outer=outer, at=max1-3*(max1-min1)/9, cex=mcex, adj=0)
      mtext(paste("UL=", round(CI[2],3)), side=3,line=4, outer=outer, at=max1-3*(max1-min1)/9, cex=mcex, adj=0)
      mtext("Monte Carlo", side=3,line=5, outer=outer, at=max1-3*(max1-min1)/9, cex=mcex, adj=0)
      if (res.asymp$SE>40*.Machine$double.eps) legend(x=max1,y= par("usr")[4],c("Monte Carlo","Asymptotic Normal"),col=c("black","blue"),lty=c(1,2),lwd=c(2:2),bty="n",title="", cex=mcex, y.intersp=.5, xpd=FALSE, xjust=.5) 
    }
    
    if(plotCI){
      yci<-par("usr")[3]-1.2*diff(par("usr")[3:4])/25
      arrows(CI[1],yci,CI[2],yci,length=0,angle=90,code=3,cex=1.5,lwd=2)
      points(quantMean,yci,pch=19,cex=1.5)
    }
    
  } #end of if plot
  
  res <- list(CI,Estimate=quantMean,SE=quantSE,MCError=quantError,p= pMinH1) #Results
  attr(res,"quant")  <- q1
  #names(res) <- c( paste((1-alpha)*100,"% ", "CI",sep=""), "Estimate", "SE","MC Error", "p")
  if(H0) return(list(res,H0res))
  else return(res)
}



confintAsymp <- function(mu, Sigma, quant=NULL, alpha=0.05, type="asymp", plot=FALSE, plotCI=FALSE){
  #Sigma <- vech.reverse(Sigma)
  q1 <- quant
  fx <- deriv(quant,names(mu),func=TRUE)
  muList <- as.list(mu)
  grad1 <- do.call("fx",muList)
  grad <- as.vector(attr(grad1,"gradient"))
  quantSE <- sqrt(t(grad)%*%Sigma%*%grad)
  quant <- parse(text=sub("~","",quant))
  quantMean <- eval(quant,muList)
  q.a <- qnorm(alpha/2,0,1)
  CI <- c(quantMean+q.a*quantSE,quantMean-q.a*quantSE)
  res <- list(CI,quantMean,quantSE)
  names(res) <- c(paste((1-alpha/2)*100,"% ", "CI",sep=""), "Estimate", "SE")
  attr(res,"quant")  <- q1
  outer <- FALSE #outer position
  mcex <- .8
  
  if (plot & quantSE>40*.Machine$double.eps)
  {
    if(type=="all"){     
      f <- Vectorize(function(x) dnorm(x,quantMean,quantSE), "x")
      curve(f,lwd=2, col="blue", lty=2, add=TRUE)
      if(plotCI){
        max1 <- par("usr")[2]
        min1 <- par("usr")[1]
        yci<-par("usr")[3]-.1*diff(par("usr")[3:4])/25
        arrows(CI[1],yci,CI[2],yci,length=0,angle=90,code=3,cex=1.5,lwd=2, col="blue",lty=2)
        points(quantMean,yci,pch=19,cex=1.5, col="blue")
        #New
        mtext(paste("LL=",round(CI[1],3)),side=3,line=1,outer=outer,at=max1-1*(max1-min1)/9,col="blue", cex=mcex, adj=0)
        mtext(paste("UL=",round(CI[2],3)),side=3,line=2,outer=outer,at=max1-1*(max1-min1)/9, col="blue" , cex=mcex, adj=0)
        mtext("Asymptotic Normal",side=3,line=3,outer=outer, at=max1-1*(max1-min1)/9, col="blue", cex=mcex, adj=0)
      }                    
    }
    
    if(type %in% c("asymp", "Asymp")){
      f <- Vectorize(function(x) dnorm(x,quantMean,quantSE),"x")
    curve(f,from=quantMean-2.58*quantSE,to=quantMean+2.58*quantSE, lwd=2, xlab="", ylab="", axes=FALSE)
    range1 <- c(quantMean-2.58*quantSE, quantMean+2.58*quantSE)
    max1 <- range1[2]
    min1 <- range1[1]
    xrange<-seq(min1,max1,length=7)
    xrange<-pretty(xrange,n=9)
    mtext(quant,1,4.5)
    axis(1,xrange,xrange, line=2.5);
    axis(2,line=1.1)
    mtext(paste("LL=",round(CI[1],3)),side=3,line=1,outer=outer,at=max1-4*(max1-min1)/9,col="black", cex=mcex)
    mtext(paste("UL=",round(CI[2],3)),side=3,line=2,outer=outer,at=max1-4*(max1-min1)/9, col="black", cex=mcex)
    mtext("Asymptotic",side=3,line=3,outer=outer, at=max1-4*(max1-min1)/9, col="black", cex=mcex)

    if(plotCI){
      yci<-par("usr")[3]-1.2*diff(par("usr")[3:4])/25
      arrows(CI[1],yci,CI[2],yci,length=0,angle=90,code=3,cex=1.5,lwd=2)
      points(quantMean,yci,pch=19,cex=1.5)
    } #if
  } 
  } #if
  
  if (plot & quantSE<=40*.Machine$double.eps){
    if(type %in% c("asymp", "Asymp")){
#     range1 <- c(quantMean-2.58*quantSE, quantMean+2.58*quantSE)
#     max1 <- range1[2]
#     min1 <- range1[1]
#     xrange<-round(seq(min1,max1,length=7),1)
#     axis(1,c(xrange, round(c(quantMean,CI),2)));axis(2)
#     smidge <- par("cin")*abs(par("tcl"))
#     text(max1-0*(max1-min1)/9,(par("usr")[4]),adj=0, bquote(mu== .(round(quantMean,3)) ))
#     text(max1-0*(max1-min1)/9,(par("usr")[4]-1*par("cxy")[2]),adj=0, bquote(sigma== .(round(quantSE,3)) ))
    cat("
        \n ******************************************* 
Note: Because the SE =",quantSE,", the plot of the sampling distribution using aymptotic normal is not available.
        \n *******************************************\n")
    
    }
    
    if(type=="all"){
      max1 <- par("usr")[2]
      min1 <- par("usr")[1]
      yci<-par("usr")[3]+2.5*diff(par("usr")[3:4])/25
      #arrows(CI[1],yci,CI[2],yci,length=0,angle=90,code=3,cex=1.5,lwd=2, col="blue",lty=2)
      #points(quantMean,yci,pch=19,cex=1.5, col="blue")
#       text(max1+0*(max1-min1)/9,(par("usr")[4]-2.5*par("cxy")[2]),adj=0, paste("LL=",round(CI[1],3)), col="blue")
#       text(max1+0*(max1-min1)/9,(par("usr")[4]-3.5*par("cxy")[2]),adj=0, paste("UL=",round(CI[2],3)), col="blue")
      mtext(paste("LL=",round(CI[1],3)),side=3,line=2,outer=outer,at=max1-1*(max1-min1)/9,col="blue", cex=mcex)
      mtext(paste("UL=",round(CI[2],3)),side=3,line=3,outer=outer,at=max1-1*(max1-min1)/9, col="blue", cex=mcex)
      mtext("Asymptotic",side=3,line=4,outer=outer, at=max1-1*(max1-min1)/9, col="blue", cex=mcex)
      cat("Because the SE,",quantSE,", is near zero, the plot of the sampling distribution using aymptotic normal is not available.\n")
      
    }
    
    }
  
  return(res)
}


