#for making multiple histograms of all variables in a matrix 
multi.histogram <- function(x) {nvar <- dim(x)[2]  #number of variables
     nsize=trunc(sqrt(nvar))+1   #size of graphic
     old.par <- par(no.readonly = TRUE) # all par settings which can be changed
     par(mfrow=c(nsize,nsize))       #set new graphic parameters
     for (i in 1:nvar) {
     name=names(x)[i]                #get the names for the variables
     hist(x[,i],main=name,xlab=name) }  #draw the histograms for each variable
     on.exit(par(old.par))   #set the graphic parameters back to the original
     }



# for pairs program to plot histograms in diagonals for scatter matrix
panel.hist <- function(x, ...)
{
    usr <- par("usr"); on.exit(par(usr))
    par(usr = c(usr[1:2], 0, 1.5) )
    h <- hist(x, plot = FALSE)
    breaks <- h$breaks; nB <- length(breaks)
    y <- h$counts; y <- y/max(y)
    rect(breaks[-nB], 0, breaks[-1], y, col = "cyan", ...)
}

# for pairs program to write correlations in upper triangle for scatter matrix
 panel.cor <- function(x, y, digits=2, prefix="", cex.cor)
     {
         usr <- par("usr"); on.exit(par(usr))
         par(usr = c(0, 1, 0, 1))
         r = (cor(x, y,use='pairwise.complete.obs'))
         txt <- format(c(r, 0.123456789), digits=digits)[1]
         txt <- paste(prefix, txt, sep="")
         if(missing(cex.cor)) cex <- 1.1
         text(0.5, 0.5, txt, cex = cex )
     }



#spss skewness
skewG1=function(x) {
  w=length(x)
  m1=mean(x)
  m2=sum((x-m1)^2)
  m3=sum((x-m1)^3)
  m4=sum((x-m1)^4)
  s1=sd(x)
  skew=w*m3/(w-1)/(w-2)/s1^3
  sdskew=sqrt( 6*w*(w-1) / ((w-2)*(w+1)*(w+3)) )
  kurtosis=(w*(w+1)*m4 - 3*m2^2*(w-1)) / ((w-1)*(w-2)*(w-3)*s1^4)
  sdkurtosis=sqrt( 4*(w^2-1) * sdskew^2 / ((w-3)*(w+5)) )
  mat=matrix(c(skew,kurtosis, sdskew,sdkurtosis), 2,
        dimnames=list(c("skew","kurtosis"), c("estimate","se")))
  return(skew)
}

#spss kurtosis
kurtG2=function(x) {
  w=length(x)
  m1=mean(x)
  m2=sum((x-m1)^2)
  m3=sum((x-m1)^3)
  m4=sum((x-m1)^4)
  s1=sd(x)
  skew=w*m3/(w-1)/(w-2)/s1^3
  sdskew=sqrt( 6*w*(w-1) / ((w-2)*(w+1)*(w+3)) )
  kurtosis=(w*(w+1)*m4 - 3*m2^2*(w-1)) / ((w-1)*(w-2)*(w-3)*s1^4)
  sdkurtosis=sqrt( 4*(w^2-1) * sdskew^2 / ((w-3)*(w+5)) )
  mat=matrix(c(skew,kurtosis, sdskew,sdkurtosis), 2,
        dimnames=list(c("skew","kurtosis"), c("estimate","se")))
  return(kurtosis)
}


#moments package results
mskew<-
function(x) {
 n=length(x)
m3 <- sum((x - mean(x))^3)/length(x)
s3 <- sqrt(sum((x-mean(x))^2)/n)^3
j<-m3/s3
return(j)
}

#moments package results but minus 3
mkurtosis<-
function(x) {
 n=length(x)
m4 <- sum((x - mean(x))^4)/length(x)
s4 <- (sum((x-mean(x))^2)/n)^2
j<-m4/s4 - 3
return(j)
}


completeFun <- function(data, desiredCols) {
  completeVec <- complete.cases(data[, desiredCols])
  return(data[completeVec, ])
}



#p1 is growth variabiliy
#sigma is residual variability
p4hlm=function(n=112,m=3,beta10=.1,sigma2=.42,p1=.04/.42){
 a=.05
 tau1=sigma2*p1
 t0=rep(1,m+1)
 t1=0:m
 z1=t1-t(t1)%*%t0/sum(t0^2)*t0
 lambda=sqrt(n)*beta10/sqrt(tau1+sigma2/sum(z1^2))
 t0=qt(1-a/2,n-1)
 1-pt(t0,n-1,lambda)+pt(-t0,n-1,lambda)
 }


#ci for factor loadings
ciload <- function(L,se,alpha) {
      alpha<- qnorm(1-alpha/2)
         x<-.5*(log((1+L)/(1-L)))
         a<-(1+L)/(1-L)
         b<-exp(alpha*((2*se)/(1-L*L))) 
        xlow<-(a/b-1)/(a/b+1)
       xhigh<-(a*b-1)/(a*b+1)
  list(xlow,xhigh) 
        }  


#ci for factor loadings with no transformation
ciload2 <- function(L,se,alpha) {
      cv<- qnorm(1-alpha/2)
       xlow<-L-se*cv
       xhigh<-L+se*cv
  list(xlow,xhigh) 
        }  


seBeta<-function (X = NULL, y = NULL, cov.x = NULL, cov.xy = NULL, var.y = NULL, 
    Nobs = NULL, alpha = 0.05, estimator = "ADF", digits = 3) 
{
    adfCOV <- function(X, y) {
        dev <- scale(cbind(X, y), scale = FALSE)
        nvar <- ncol(dev)
        N <- nrow(dev)
        ue <- nvar * (nvar + 1)/2
        s <- vector(length = ue, mode = "character")
        z <- 0
        for (i in 1:nvar) {
            for (j in i:nvar) {
                z <- z + 1
                s[z] <- paste(i, j, sep = "")
            }
        }
        v <- expand.grid(s, s)
        V <- paste(v[, 1], v[, 2], sep = "")
        id.mat <- matrix(0, nrow = ue^2, 4)
        for (i in 1:4) id.mat[, i] <- as.numeric(sapply(V, substr, 
            i, i))
        M <- matrix(1:ue^2, ue, ue, byrow = TRUE)
        r <- M[lower.tri(M, diag = TRUE)]
        ids <- id.mat[r, ]
        adfCovMat <- matrix(0, ue, ue)
        covs <- matrix(0, nrow(ids), 1)
        for (i in 1:nrow(ids)) {
            w_ij <- crossprod(dev[, ids[i, 1]], dev[, ids[i, 
                2]])/N
            w_ik <- crossprod(dev[, ids[i, 1]], dev[, ids[i, 
                3]])/N
            w_il <- crossprod(dev[, ids[i, 1]], dev[, ids[i, 
                4]])/N
            w_jk <- crossprod(dev[, ids[i, 2]], dev[, ids[i, 
                3]])/N
            w_jl <- crossprod(dev[, ids[i, 2]], dev[, ids[i, 
                4]])/N
            w_kl <- crossprod(dev[, ids[i, 3]], dev[, ids[i, 
                4]])/N
            w_ijkl <- (t(dev[, ids[i, 1]] * dev[, ids[i, 2]]) %*% 
                (dev[, ids[i, 3]] * dev[, ids[i, 4]])/N)
            covs[i] <- (N * (N - 1) * (1/((N - 2) * (N - 3))) * 
                (w_ijkl - w_ij * w_kl) - N * (1/((N - 2) * (N - 
                3))) * (w_ik * w_jl + w_il * w_jk - (2/(N - 1)) * 
                w_ij * w_kl))
        }
        adfCovMat[lower.tri(adfCovMat, diag = T)] <- covs
        vars <- diag(adfCovMat)
        adfCovMat <- adfCovMat + t(adfCovMat) - diag(vars)
        adfCovMat
    }
    vech <- function(x) t(x[!upper.tri(x)])
    Dn <- function(x) {
        mat <- diag(x)
        index <- seq(x * (x + 1)/2)
        mat[lower.tri(mat, TRUE)] <- index
        mat[upper.tri(mat)] <- t(mat)[upper.tri(mat)]
        outer(c(mat), index, function(x, y) ifelse(x == y, 1, 
            0))
    }
    DIAG <- function(x = 1, nrow, ncol) {
        if (length(x) == 1) 
            x <- as.matrix(x)
        if (is.matrix(x)) 
            return(diag(x))
        else return(diag(x, nrow, ncol))
    }
    if ((is.null(X) | is.null(y)) & estimator == "ADF") 
        stop("\nYou need to supply both X and y for ADF Estimation\n")
    if (is.null(X) & !is.null(y)) 
        stop("\n y is not defined\n Need to specify both X and y\n")
    if (!is.null(X) & is.null(y)) 
        stop("\n X is not defined\n Need to specify both X and y\n")
    if (is.null(X) & is.null(y)) {
        if (is.null(cov.x) | is.null(cov.xy) | is.null(var.y) | 
            is.null(Nobs)) 
            stop("\nYou need to specify covariances and sample size\n")
        scov <- rbind(cbind(cov.x, cov.xy), c(cov.xy, var.y))
        N <- Nobs
        p <- nrow(cov.x)
    }
    else {
        X <- as.matrix(X)
        y <- as.matrix(y)
        scov <- cov(cbind(X, y))
        N <- length(y)
        p <- ncol(X)
    }
    if (estimator == "ADF") {
        cov.cov <- adfCOV(X, y)
    }
    else {
        Kp.lft <- solve(t(Dn(p + 1)) %*% Dn(p + 1)) %*% t(Dn(p + 
            1))
        cov.cov <- 2 * Kp.lft %*% (scov %x% scov) %*% t(Kp.lft)
    }
    param <- c(vech(scov))
    ncovs <- length(param)
    if (p == 1) {
        v.x.pl <- 1
    }
    else {
        v.x.pl <- c(1, rep(0, p - 1))
        for (i in 2:p) v.x.pl[i] <- v.x.pl[i - 1] + p - (i - 
            2)
    }
    cx <- scov[1:p, 1:p]
    cxy <- scov[1:p, p + 1]
    vy <- scov[p + 1, p + 1]
    sx <- sqrt(DIAG(cx))
    sy <- sqrt(vy)
    bu <- solve(cx) %*% cxy
    ncx <- length(vech(cx))
    db <- matrix(0, p, ncovs)
    V <- matrix(0, p, ncx)
    V[as.matrix(cbind(1:p, v.x.pl))] <- 1
    db[, 1:ncx] <- (DIAG(c(solve(DIAG(2 * sx * sy)) %*% bu)) %*% 
        V - DIAG(sx/sy) %*% (t(bu) %x% solve(cx)) %*% Dn(p))
    db[, (ncx + 1):(ncx + p)] <- DIAG(sx/sy) %*% solve(cx)
    db[, ncovs] <- -DIAG(sx/(2 * sy^3)) %*% bu
    cx.nms <- matrix(0, p, p)
    cxy.nms <- c(rep(0, p), "var_y")
    for (i in 1:p) for (j in 1:p) cx.nms[i, j] <- paste("cov_x", 
        i, "x", j, sep = "")
    for (i in 1:p) cxy.nms[i] <- paste("cov_x", i, "y", sep = "")
    old.ord <- c(vech(cx.nms), cxy.nms)
    new.ord <- vech(rbind(cbind(cx.nms, cxy.nms[1:p]), c(cxy.nms)))
    db <- db[, match(new.ord, old.ord)]
    if (p == 1) 
        DEL.cmat <- t(db) %*% cov.cov %*% db/N
    else DEL.cmat <- db %*% cov.cov %*% t(db)/N
    b.nms <- NULL
    for (i in 1:p) b.nms[i] <- paste("beta_", i, sep = "")
    rownames(DEL.cmat) <- colnames(DEL.cmat) <- b.nms
    DELse <- sqrt(DIAG(DEL.cmat))
    CIs <- as.data.frame(matrix(0, p, 3))
    colnames(CIs) <- c("lbound", "estimate", "ubound")
    for (i in 1:p) rownames(CIs)[i] <- paste("beta_", i, sep = "")
    tc <- qt(alpha/2, N - p - 1, lower.tail = F)
    beta <- DIAG(sx) %*% bu * sy^-1
    for (i in 1:p) {
        CIs[i, ] <- c(beta[i] - tc * DELse[i], beta[i], beta[i] + 
            tc * DELse[i])
    }
    cat("\n", 100 * (1 - alpha), "% CIs for Standardized Regression Coefficients:\n\n", 
        sep = "")
    # print(round(CIs, digits))
    invisible(list(cov.mat = DEL.cmat, SEs = DELse, alpha = alpha, 
        CIs = CIs, estimator = estimator))
}


#intraclass correlations

ICCest <-function (x, y, data = NULL, alpha = 0.05, CI.type = c("THD", 
    "Smith")) 
{
    square <- function(z) {
        z^2
    }
    icall <- list(y = substitute(y), x = substitute(x))
    if (is.character(icall$y)) {
        warning("passing a character string to 'y' is deprecated since ICC vesion 2.3.0 and will not be supported in future versions. The argument to 'y' should either be an unquoted column name of 'data' or an object")
        if (missing(data)) 
            stop("Supply either the unquoted name of the object containing 'y' or supply both 'data' and then 'y' as an unquoted column name to 'data'")
        icall$y <- eval(as.name(y), data, parent.frame())
    }
    if (is.name(icall$y)) 
        icall$y <- eval(icall$y, data, parent.frame())
    if (is.call(icall$y)) 
        icall$y <- eval(icall$y, data, parent.frame())
    if (is.character(icall$y)) 
        icall$y <- eval(as.name(icall$y), data, parent.frame())
    if (is.character(icall$x)) {
        warning("passing a character string to 'x' is deprecated since ICC vesion 2.3.0 and will not be supported in future versions. The argument to 'x' should either be an unquoted column name of 'data' or an object")
        if (missing(data)) 
            stop("Supply either the unquoted name of the object containing 'x' or supply both 'data' and then 'x' as an unquoted column name to 'data'")
        icall$x <- eval(as.name(x), data, parent.frame())
    }
    if (is.name(icall$x)) 
        icall$x <- eval(icall$x, data, parent.frame())
    if (is.call(icall$x)) 
        icall$x <- eval(icall$x, data, parent.frame())
    if (is.character(icall$x) && length(icall$x) == 1) 
        icall$x <- eval(as.name(icall$x), data, parent.frame())
    tdata <- data.frame(icall)
    tdata <- na.omit(tdata)
    a <- length(unique(tdata$x))
    if (!is.null(attributes(tdata)$na.action)) {
        warning(cat("NAs removed from rows:\n", unclass(attributes(tdata)$na.action), 
            "\n"))
    }
    if (!is.factor(tdata$x)) {
        warning("'x' has been coerced to a factor")
        tdata$x <- as.factor(tdata$x)
    }
    else {
        if (length(levels(tdata$x)) > a) {
            tdata$x <- factor(as.character(tdata$x), levels = unique(tdata$x))
            warning("Missing levels of 'x' have been removed")
        }
    }
    tmpbb <- anova(aov(y ~ x, data = tdata))
    num.df <- tmpbb$Df[1]
    denom.df <- tmpbb$Df[2]
    MSa <- tmpbb$"Mean Sq"[1]
    MSw <- var.w <- tmpbb$"Mean Sq"[2]
    tmp.outj <- aggregate(y ~ x, data = tdata, FUN = length)$y
    k <- (1/(a - 1)) * (sum(tmp.outj) - (sum(square(tmp.outj))/sum(tmp.outj)))
    var.a <- (MSa - MSw)/k
    r <- var.a/(var.w + var.a)
    low.F <- qf(alpha/2, num.df, denom.df, lower.tail = FALSE)
    N <- nrow(tdata)
    n.bar <- N/a
    n.not <- n.bar - sum(square(tmp.outj - n.bar)/((a - 1) * 
        N))
    type <- match.arg(CI.type)
    if (type == "THD") {
        up.F <- qf(alpha/2, denom.df, num.df, lower.tail = FALSE)
        FL <- (MSa/MSw)/low.F
        FU <- (MSa/MSw) * up.F
        low.CI <- (FL - 1)/(FL + n.not - 1)
        up.CI <- (FU - 1)/(FU + n.not - 1)
    }
    if (type == "Smith") {
        z.not <- qnorm(alpha/2)
        Vr <- (2 * square(1 - r)/square(n.not)) * ((square((1 + 
            r * (n.not - 1)))/(N - a)) + ((a - 1) * (1 - r) * 
            (1 + r * (2 * n.not - 1)) + square(r) * (sum(square(tmp.outj)) - 
            2 * (1/N) * sum((tmp.outj^3)) + (1/square(N)) * square(sum(square(tmp.outj)))))/square(a - 
            1))
        low.CI <- r + z.not * sqrt(Vr)
        up.CI <- r - z.not * sqrt(Vr)
    }
    list(ICC = r, LowerCI = low.CI, UpperCI = up.CI, N = a, k = k, 
        varw = var.w, vara = var.a)
}


#bootstrap an interpolated median
imedianCI <- function(x, w=1, na.rm=TRUE, conf.level=0.95, bootn=1000) {
#  library(psych)
  bootmed = apply(X=matrix(sample(x, rep=TRUE, bootn*length(x)), nrow=bootn), MARGIN = 1, FUN = interp.median, w=w, na.rm=TRUE);
  CI <- quantile(bootmed, c((1-conf.level)/2, conf.level+(1-conf.level)/2));
  interpmed <- interp.median(x=x, w=w, na.rm=na.rm);
  return(list("ci"=CI));
}


#calculate MAD using interpolated median
imad <- function(x,c=1) {
#library(psych)
med <- interp.median(x=x, w=1)
jj<-abs(x-med)
jj<-(interp.median(jj))*c
  return(jj);
}


#khb recaling adjustments in logistic/probit modeling
paste0 <- function(..., collapse = NULL) {
    paste(..., sep = "", collapse = collapse)
}
.getNamespace <- function(name) .Internal(getRegisteredNamespace(as.name(name)))
khb <- function(reduced, full, corrected.se=FALSE, med.sandwich=NULL, glm.sandwich=FALSE,  savefs=FALSE){
	## Variables in reduced model
	keyvar <- attr(reduced$terms, "term.labels")
	## Variables in full model
	vf <- attr(full$terms, "term.labels")
	## Define mediators as variable present only in full model
	mediators <- vf[!(vf %in% keyvar)]
	message(" [>] Mediators: ", paste(mediators, collapse=", "))
	datamm <- model.matrix(full)
	if(inherits(full, "clm")){
		datamm <- datamm$X
	}
	datammdf <- as.data.frame(datamm)
	## Some mediators may have several z variable (ie factors)
	getVarList <- function(kv){
		## position of variables in term assign list
		ikv <- match(kv, vf)
		kvList <- list()
		for(i in seq_along(kv)){
			kvList[[kv[i]]] <- colnames(datamm)[attr(datamm, "assign") == ikv[i]]
		}
		return(kvList)
	}
	keyvarList <- getVarList(keyvar)
	mediatorList <- getVarList(mediators)
	allkeyvar <- unlist(keyvarList, use.names=F) ##colnames(datamm)[attr(datamm, "assign") %in% ikeyvar]
	allmediators <- unlist(mediatorList, use.names=F) ##colnames(datamm)[attr(datamm, "assign") %in% imediators]
	nz <- length(allmediators)
	## Choosing covariance estimator
	if(is.null(med.sandwich)){
		med.sandwich <- any(sapply(datammdf[, allmediators], function(x){return(length(unique(x))==2)}))
	}
	stata.sandwich <- function(...){
		return(sandwich(..., adjust=TRUE))
	}
	vcovfunc <- if(med.sandwich){stata.sandwich} else{ vcov}
	vcovglm <- if(glm.sandwich){stata.sandwich} else{vcov}
	
	## Return object
	ret <- list()
	## Storing summary informations.
	ff <- list()
	for(z in allmediators){
		ff[[z]] <- as.formula(paste0("`",z,"`~`",paste(c(allkeyvar), collapse="`+`"),"`"))
	}
	##getAdjustedModel <- function(reduced){
		## Try to get the same as stata
	fitsur <- systemfit(ff, data=datammdf, method="SUR", control=systemfit.control(maxiter = 2, methodResidCov = "noDfCor"))
	## Compute residuals 
	zresid <- resid(fitsur)
	names(zresid) <- paste("KHB__KHB__", names(zresid), sep="")
	## Updating the reduced model to include the 
	adjusted <- update(reduced, formula = as.formula(paste0(".~.+", paste0("`KHB__KHB__", allmediators,"`", collapse=" + "))), data = data.frame(model.frame(reduced), zresid))
	
	##	return(adjusted)
	##}
	##adjusted <- getAdjustedModel(reduced)
	## Needed for statistical signif
	## sgrzr <- summary(grzr)
	## sgf <- summary(gf)
	zcovmat <- matrix(0, ncol=2*nz, nrow=2*nz)
	zcovmat[1:nz, 1:nz] <- vcovglm(full)[allmediators, allmediators]
	if(savefs){
		estimation <- list()
	}
	adjustedCoef  <- coeftest(adjusted, vcov.=vcovglm)
	fullCoef <- coeftest(full, vcov.=vcovglm)
	reducedCoef <- coef(reduced)
	ckgr <- coef(reduced)[allkeyvar]
	ckgrzr <- adjustedCoef[allkeyvar, 1]
	ckgf <- fullCoef[allkeyvar, 1]
	sum_conf <- data.frame(Ratio=ckgrzr/ckgf, Percentage=100*(ckgrzr-ckgf)/ckgrzr, 
							Rescaling=ckgrzr/ckgr)
	rownames(sum_conf) <- allkeyvar
	
	for(kk in allkeyvar){
		ret[[kk]] <- list()
		medSUR <- match(paste(allmediators, kk, sep="_"), names(coef(fitsur)))
		zcoef <- coef(fitsur)[medSUR]
		names(zcoef) <- allmediators
		zcovmat[nz+1:nz, nz+1:nz] <- vcovfunc(fitsur)[medSUR, medSUR]
		a <- c(zcoef, fullCoef[allmediators, 1])
		sdconf <-  sqrt(t(a) %*% zcovmat %*% a)
		if(savefs){
			estimation[[kk]] <- list(a=a, sigma=zcovmat, se=sdconf) 
		}

		res <- matrix(0, ncol=4, nrow=3)
		colnames(res) <- c("Estimate", "Std. Error", "z value", "Pr(>|z|)")
		rownames(res) <- c("Reduced", "Full", "Diff")
		
		res[1, ] <- adjustedCoef[kk, ]
		if(corrected.se){
			SIGMA <- matrix(0, nrow=2*nz+1, ncol=2*nz+1)
			intvar <- c(kk, allmediators)
			aa <- c(zcoef, fullCoef[intvar, 1])
			SIGMA[1:(1+nz), 1:(1+nz)] <- vcovglm(full)[intvar, intvar]
			SIGMA[nz+1+1:nz,nz+1+1:nz] <- zcovmat[nz+1:nz, nz+1:nz]
			res[1, 2] <- sqrt(t(aa) %*% SIGMA %*% aa)
			res[1, 3] <- res[1, 1]/res[1, 2]
			res[1,4] <- 2 * pnorm(-abs(res[1, 3]))
			
		}
		res[2, ] <- fullCoef[kk, ]
		res[3, ] <- res[1,]-res[2, ]
		res[3, 2] <- sdconf
		res[3, 3] <- res[3, 1]/sdconf
		res[3,4] <- 2 * pnorm(-abs(res[3, 3]))
		ret[[kk]]$table <- res
	
		res <- matrix(0, ncol=6, nrow=1+nz)
		colnames(res) <- c("% diff", "% effect", "Estimate", "Std. Error", "z value", "Pr(>|z|)")
		rownames(res) <- c("Grand Indirect", allmediators)
		res["Grand Indirect", ] <- c(100, 100*sum_conf[kk, "Percentage"], ret[[kk]]$table[3, 1:4])
		for(zz in allmediators){
			izz <- match(zz, allmediators)
			res[zz, 3] <- fullCoef[zz, 1]*zcoef[zz]
			izzpos <- c(izz, nz+izz)
			res[zz, 4] <- sqrt(t(a[izzpos]) %*% zcovmat[izzpos,izzpos] %*% a[izzpos])
			res[zz, 5] <- res[zz, 3]/res[zz, 4]
			res[zz, 6] <- 2 * pnorm(-abs(res[zz, 5]))
			res[zz, 1] <- 100*res[zz, 3]/ res[1, 3]
			res[zz, 2] <- 100*res[zz, 3]/ ret[[kk]]$table[1,1]
		}
		ret[[kk]]$detail <- res
		##ret[[kk]]$ape["Diff"] <- ret$ape["Reduced"] - ret$ape["Full"]
	}
	finalobj <-list(sum_conf=sum_conf, key=ret, reduced=reduced, adjusted=adjusted, full=full,
		keyvar=keyvarList, mediators=mediatorList)
	if(savefs){
		estimation$fitsur <- fitsur
		finalobj$estimation <- estimation
	}
	class(finalobj) <- "khb"
	return(finalobj)

}
print.khb <- function(x, type="summary", keyvar=NULL, disentangle=FALSE, ...){
	if(type=="models"){
		args <- c(list(Reduced=x$reduced, Adjusted=x$adjusted, Full=x$full), list(...))
		print(do.call(regTable, args), quote=FALSE)
		return(invisible())
	}
	
	cat("KHB method\n")
	cat("Model type:",  class(x$reduced), ifelse(inherits(x$reduced, "glm"),paste0("(", family(x$reduced)$link,")"), "") , "\n")
	varList <- function(vl){
		ret <- character(length(vl))
		names(ret) <- names(vl)
		for(nn in names(vl)){
			vlnn <- vl[[nn]]
			if(length(vlnn)==1 && vlnn==nn){
				ret[nn] <- nn
			}else{
				ret[nn] <- paste0(nn, " (", paste(vlnn, collapse=", "), ")")
			}
		}
		return(paste(ret, collapse=", "))
	}
	cat("Variables of interest: ", varList(x$keyvar), fill=TRUE, sep="")
	cat("Z variables (mediators): ", varList(x$mediators), fill=TRUE, sep="")
	cat("\nSummary of confounding\n")
	printCoefmat(x$sum_conf)
	#cat("\n")
	
	if(is.null(keyvar)){
		keyvar <- names(x$key)
	}else if(length(keyvar)==1 && keyvar==FALSE){
		return(invisible())
	} else{
		cond <- keyvar %in% names(x$key)
		cond2 <- keyvar %in% names(x$keyvar)
		if(any((!cond) & (!cond2))){
			stop(" [!] The following keyvar were not found: ", paste(keyvar[(!cond) & (!cond2)], collapse=", "))
		}
		keyvar <- c(keyvar[cond], unlist(x$keyvar[keyvar[cond2&(!cond)]]))
	}
	for(k in keyvar){
		cat(rep("--", 22),"\n", sep="-")
		cat(k,":\n")
		printCoefmat(x$key[[k]]$table, cs.ind=1:2, tst.ind=3, signif.legend=FALSE)
		if(nrow(x$key[[k]]$detail)>2 && disentangle){
			cat("\n\nDecomposition of the difference:\n")
			printCoefmat(x$key[[k]]$detail[-1,], cs.ind=3:4, tst.ind=5, zap.ind=1:2)
		}
	}
	return(invisible())
}

#more khb
apecoef <- function(x){
	if(family(x)$link=="probit"){
		dform <- dnorm
	}else{
		dform <- dlogis
	}
	return(mean(dform(predict(x, type="link")))*coef(x))
}


getRatioByType <- function(gr, gf, ratio.type, coefFunc=coef){
	cr <- coefFunc(gr)
	cf <-coefFunc(gf)
	uc <- names(cr)[names(cr) %in% names(cf)]
	cr <- cr[uc]
	cf <-cf[uc]
	if(ratio.type=="perc"){
		arc <- 100*(cr-cf)/cr
	} else if (ratio.type=="diff"){
		arc <- cf-cr
	} else{
		arc <- cr/cf
	}
	return(arc)
}


coefChange <- function(reduced, full, method="ystand", ratio.type="perc", R=1){
	method <- tolower(method)
	if(method=="ape"){
		coefFunc <- apecoef
	} else if(method=="naive"){
		coefFunc <- coef
	}else if(method=="ystand"){
		coefFunc <- ystandcoef
	} else if(method=="khb"){
		x <- khb(reduced, full)
		reduced <- x$adjusted
		coefFunc <- coef
	}
	
	if(R>1){
		coefFunc <- cmpfun(coefFunc)
		coefChangeInternal <- cmpfun(function(dat, ww){
			ggr <- update(reduced, weights=ww)
			ggf <- update(full, weights=ww)
			return(getRatioByType(ggr, ggf, ratio.type=ratio.type, coefFunc=coefFunc))
		})
		bts <- boot(1:nrow(reduced$model), statistic=coefChangeInternal, sim="ordinary", R=R, stype="f")
		return(bts)
	}

	return(getRatioByType(reduced, full, ratio.type=ratio.type, coefFunc=coefFunc))
}


ystandcoef <- function(x){
	const <- (pi^2)/3
	if(family(x)$link=="probit"){
		const <- 1
	}
	sdy <- sqrt(var(predict(x, type="link"))+const)
	xc <- coef(x)/sdy
	attr(xc, "sdy") <- sdy
	return(xc)
}

compareModels <- function(reduced, full, method="ystand", ratio.type="perc"){
	method <- tolower(method)
	if(method=="ape"){
		coefFunc <- apecoef
	} else if(method=="naive"){
		coefFunc <- coef
	}else if(method=="ystand"){
		coefFunc <- ystandcoef
	} else if(method=="khb"){
		x <- khb(reduced, full)
		reduced <- x$adjusted
		coefFunc <- coef
	}
	rtable <- regTable(Reduced=reduced, Full=full, coefFunc=coefFunc, 
				signif.stars=FALSE, se=FALSE, ratio.type=ratio.type, 
				redum="Reduced", fullm="Full")
	return(rtable)
}

regTable <- function(..., digits=max(3L, getOption("digits") - 2L), signif.stars = getOption("show.signif.stars"), 
    signif.legend = signif.stars, se=TRUE, coefFunc=coef, ratio.type=NULL, redum=1, fullm=2, khb.resid=TRUE){
	sigfunc <- function(pv){
		return(symnum(pv, corr = FALSE, na = FALSE, 
                  cutpoints = c(0, 0.001, 0.01, 0.05, 0.1, 1), 
                  symbols = c("***", "**", "*", ".", " ")))
	}
	roundedDigit <- function(val){
		val <- zapsmall(val, digits)
		acs <- abs(val)
		if (any(ia <- is.finite(acs))) {
			digmin <- 1 + if (length(acs <- acs[ia & acs != 0])) 
				floor(log10(range(acs[acs != 0], finite = TRUE)))
			else 0
			xx <- format(round(val, max(1L, digits - 
				digmin)), digits = digits)
		}else{
			xx <- format(val, digits = digits)
		}
		return(xx)
	}
	coeffSignif <- function(rl){
		if(signif.stars){
			allcof <- summary(rl)$coefficients
			cof <- allcof[, 1]
			significance <- allcof[, 4]
			mycoeff <- paste(
					roundedDigit(cof),
					sigfunc(significance)
					)
			names(mycoeff) <- names(cof)
		}
		else{
			mycoeff <- roundedDigit(coefFunc(rl))
		}
		if(se){
			nc <- length(mycoeff)
			withse <- character(2*nc)
			withse[(1:nc)*2-1] <- mycoeff
			names(withse)[(1:nc)*2-1] <- names(mycoeff)
			withse[(1:nc)*2] <- paste0("(", roundedDigit(allcof[,2]), ")")
			names(withse)[(1:nc)*2] <- paste("SE_", names(mycoeff), sep="")
			return(withse)
		}else{
			return(mycoeff)
		}
		
		
	}
	allreg <- list(...)
	allregraw <- allreg
	if(is.null(names(allreg))){
		names(allreg) <- paste("(", 1:length(allreg), ")", sep="")
	}
	allregcoef <- list()
	diffcoef <- NULL
	for(rr in names(allreg)){
		allregcoef[[rr]] <- coeffSignif(allreg[[rr]])
		diffcoef <- c(diffcoef, names(allregcoef[[rr]]))
	}
	diffcoef <- unique(diffcoef)
	regmat <- matrix("", nrow=length(diffcoef), ncol=length(allreg)+!is.null(ratio.type))
	rownames(regmat) <- diffcoef
	##colnames(regmat) <- names(allreg)
	for(rr in 1:length(allreg)){
		arc <- allregcoef[[names(allreg)[rr]]]
		regmat[names(arc), rr] <- arc
	}
	if(!is.null(ratio.type)){
		arc <- getRatioByType(allreg[[redum]], allreg[[fullm]], ratio.type=ratio.type, coefFunc=coefFunc)
		regmat[names(arc), length(allreg)+1] <- roundedDigit(arc)
	}
	if(se){
		rownames(regmat)[(1:(nrow(regmat)/2))*2] <- ""
	}
	
	rowcond <- grepl("KHB__KHB__", rownames(regmat), fixed=TRUE)
	rownames(regmat)[rowcond] <- gsub("\\`?KHB__KHB__(.*)\\`?$", "Resid(\\1)", rownames(regmat)[rowcond])
	if(!khb.resid){
		regmat <- regmat[!rowcond, ]
	}
	qualmat <- matrix("", ncol=length(allreg), nrow=7)
	rownames(qualmat) <- c("Pseudo R2", "Dev.", "Null", "Chisq", "Sig", "Dl", "BIC")
	qualmat["Pseudo R2", ] <- sapply(allreg, function(x) roundedDigit(NagelkerkeR2(x)$R2))
	qualmat["Dev.", ] <- sapply(allreg, function(x) roundedDigit(deviance(x)))
	qualmat["Null", ] <- sapply(allreg, function(x) roundedDigit(x$null.deviance))
	qualmat["Chisq", ] <- sapply(allreg, function(x) {sx <- summary(x); roundedDigit(sx$null.deviance - sx$deviance)})
	qualmat["Sig", ] <- sapply(allreg, function(x) {sx <- summary(x);
											dl <- sx$df.null - sx$df.residual;
											sigfunc(pchisq(sx$null.deviance - sx$deviance, dl, lower.tail=FALSE))
											}
	)
	qualmat["Dl", ] <- sapply(allreg, function(x) {sx <- summary(x);  paste(sx$df.null - sx$df.residual, sep="")})
	qualmat["BIC", ] <- sapply(allreg, function(x) {sx <- summary(x); roundedDigit(sx$deviance+(sx$df.null - sx$df.residual)*log(nobs(x)))})
	
	finalmat <- rbind(rep("", length(allreg)), qualmat)
	if(!is.null(ratio.type)) finalmat <- cbind(finalmat, rep("", nrow(finalmat)))
	finalmat <- rbind(regmat, finalmat)
	colnames(finalmat) <- c(names(allreg), ratio.type)
	return(finalmat)
}

NagelkerkeR2<- function (rr) 
{
    n <- nrow(rr$model)
    R2 <- (1 - exp((rr$dev - rr$null)/n))/(1 - exp(-rr$null/n))
    RVAL <- list(N = n, R2 = R2)
    return(RVAL)
}



#Millsap sem selection stuff
pnormmix <- function(q, mean1 = 0, sd1 = 1, mean2 = 0, sd2 = 1, pmix1 = 0.5,
lower.tail = TRUE) {
# Distribution function (pdf) of a mixture of two normal distributions.
#
# Args:
# q: vector of quantiles.
# mean1: mean of first normal distribution.
# sd1: standard deviation of first normal distribution.
# mean2: mean of second normal distribution.
# sd2: standard deviation of second normal distribution.
# pmix1: mixing proportion for the first distribution. Should be a
# number in the range (0, 1).
# lower.tail: logical; if TRUE(default), probabilities are P[X <= x];
# otherwise, P[X > x].
#
# Returns:
# Cumulative probability of q or 1 minus it on the mixture normal
# distribution.
# Error handling
stopifnot(pmix1 > 0, pmix1 < 1)
as.vector(c(pmix1, 1 - pmix1) %*%
sapply(q, pnorm, mean = c(mean1, mean2), sd = c(sd1, sd2),
lower.tail = lower.tail))
}
# test
# pnormmix(1, 0, 3.1, 1.7, 3.1, lower.tail = FALSE)
qnormmix <- function(p, mean1 = 0, sd1 = 1, mean2 = 0, sd2 = 1, pmix1 = 0.5,
lower.tail = TRUE) {
# Quantile function of a mixture of two normal distributions.
#
# Args:
# p: vector of probabilities.
# mean1: mean of first normal distribution.
# sd1: standard deviation of first normal distribution.
# mean2: mean of second normal distribution.
# sd2: standard deviation of second normal distribution.
# pmix1: mixing proportion for the first distribution. Should be a
# number in the range (0, 1).
# lower.tail: logical; if TRUE(default), probabilities are P[X <= x];
# otherwise, P[X > x].
#
# Returns:
# Quantile corresponding to p or 1 - p on the mixture normal
# distribution.
# Error handling
stopifnot(pmix1 > 0, pmix1 < 1, p >= 0, p <= 1)
f <- function(x) (pnormmix(x, mean1, sd1, mean2, sd2, pmix1,
lower.tail) - p)^2
start <- as.vector(c(pmix1, 1 - pmix1) %*%
sapply(p, qnorm, c(mean1, mean2), c(sd1, sd2),
lower.tail = lower.tail))
nlminb(start, f)$par
}
.bvnorm_kernel <- function(x, y, mu_x = 0, mu_y = 0, sd_x = 1, sd_y = 1,
cov_xy = 0) {
# Helper funcction for computing the kernel for bivariate normal density
cor <- cov_xy / sd_x / sd_y
numer <- (x - mu_x)^2 / sd_x^2 + (y - mu_y)^2 / sd_y^2 -
2 * cor * (x - mu_x) * (y - mu_y) / sd_x / sd_y
numer / (1 - cor^2)
}
contour_bvnorm <- function(mean1 = 0, sd1 = 1, mean2 = 0, sd2 = 1,
cor12 = 0, cov12 = NULL, density = .95, length_out = 101, bty = "L",...)
 {
# Plot contour for a bivariate normal distribution
#
# Args:
# mean1: mean of first normal distribution (on x-axis).
# sd1: standard deviation of first normal distribution.
# mean2: mean of second normal distribution (on y-axis).
# sd2: standard deviation of second normal distribution.
# cor12: correlation in the bivariate normal.
# cov12: covariance in the bivariate normal. If not input, compute the
# covariance using the correlation and the standard deviations.
# density: density level, i.e., probability enclosed by the ellipse.
# length_out: number of values on the x-axis and on the y-axis to be
# evaluated; default to 101.
# bty: argument passed to the `contour` function.
# ...: other arguments passed to the `countour` funcction.
#
# Returns:
# a plot showing the contour of the bivariate normal distribution on
# a two-dimensional space.
# Error handling
stopifnot(cor12 >= -1, cor12 <= 1)
if (is.null(cov12)) cov12 <- cor12 * sd1 * sd2s
x_seq <- mean1 + seq(-3, 3, length.out = length_out) * sd1
y_seq <- mean2 + seq(-3, 3, length.out = length_out) * sd2
z <- outer(x_seq, y_seq, .bvnorm_kernel, mu_x = mean1, mu_y = mean2,
sd_x = sd1, sd_y = sd2, cov_xy = cov12)
contour(x_seq, y_seq, z, levels = qchisq(density, 2), drawlabels = FALSE,
bty = bty, ...)
}
.partit_bvnorm <- function(cut1, cut2, mean1 = 0, sd1 = 1, mean2 = 0, sd2 = 1,
cor12 = 0, cov12 = cor12 * sd1 * sd2) {
# Helper function for computing summary statistics from a selection approach
Sigma <- matrix(c(sd1^2, cov12, cov12, sd2^2), nrow = 2)
C <- pmnorm(c(cut1, cut2), c(mean1, mean2), Sigma)
B <- pnorm(cut1, mean1, sd1) - C
D <- pnorm(cut2, mean2, sd2) - C
A <- 1 - B - C - D
propsel <- A + B
success_ratio <- A / propsel
sensitivity <- A / (A + D)
specificity <- C / (C + B)
c(A, B, C, D, propsel, success_ratio, sensitivity, specificity)
}
PartInv <- function(propsel, cut_z = NULL, kappa_r, kappa_f = kappa_r,
phi_r, phi_f = phi_r, lambda_r, lambda_f = lambda_r,
Theta_r, Theta_f = Theta_r, tau_r, tau_f = tau_r,
pmix_ref = 0.5, plot_contour = TRUE, ...) {
# Evaluate partial measurement invariance using Millsap & Kwok's (2004)
# approach
#
# Args:
# propsel: proportion of selection. If missing, computed using `cut_z`.
# cut_z: prespecified cutoff score on the observed composite. This argument
# is ignored when `propsel` has input.
# kappa_r: latent factor mean for the reference group.
# kappa_f: (optional) latent factor mean for the focal group;
# if no input, set equal to kappa_r.
# phi_r: latent factor variance for the reference group.
# phi_f: (optional) latent factor variance for the focal group;
# if no input, set equal to phi_r.
# lambda_r: a vector of factor loadings for the reference group.
# lambda_f: (optional) a vector of factor loadings for the focal group;
# if no input, set equal to lambda_r.
# tau_r: a vector of measurement intercepts for the reference group.
# tau_f: (optional) a vector of measurement intercepts for the focal group;
# if no input, set equal to tau_r.
# Theta_r: a matrix of the unique factor variances and covariances
# for the reference group.
# Theta_f: (optional) a matrix of the unique factor variances and
# covariances for the focal group;
# if no input, set equal to Theta_r.
# pmix_ref: Proportion of the reference group;
# default to 0.5 (i.e., two populations have equal size).
# plot_contour: logical; whether the contour of the two populations
# should be plotted; default to TRUE.
# ...: other arguments passed to the `countour` funcction.
#
# Returns:
# a list of four elements and a plot if plot_contour == TRUE:
# - propsel: echo the same argument as input.
# - cutpt_xi: cut point on the latent scale (xi).
# - cutpt_z: cut point on the observed scale (Z).
# - summary: A 8 x 2 table, with columns representing the reference
# and the focal groups, and the rows represent
# probabilities of true positive (A), false positive (B),
# true negative (C), false negative (D); proportion selected,
# success ratio, sensitivity, and specificity.
# Error handling
stopifnot(length(kappa_r) == 1, length(kappa_f) == 1, length(phi_r) == 1,
length(phi_f) == 1)
if (length(Theta_r) == length(lambda_r)) Theta_r <- diag(Theta_r)
if (length(Theta_f) == length(lambda_f)) Theta_f <- diag(Theta_f)
library(mnormt) # load `mnormt` package
mean_zr <- sum(tau_r) + sum(lambda_r) * kappa_r
mean_zf <- sum(tau_f) + sum(lambda_f) * kappa_f
sd_zr <- sqrt(sum(lambda_r)^2 * phi_r + sum(Theta_r))
sd_zf <- sqrt(sum(lambda_f)^2 * phi_f + sum(Theta_f))
cov_z_xir <- sum(lambda_r) * phi_r
cov_z_xif <- sum(lambda_f) * phi_f
sd_xir <- sqrt(phi_r)
sd_xif <- sqrt(phi_f)
if (!missing(propsel)) {
if (!is.null(cut_z)) {
warning("Input to `cut_z` is ignored.")
}
cut_z <- qnormmix(propsel, mean_zr, sd_zr, mean_zf, sd_zf,
pmix_ref, lower.tail = FALSE)
} else if (!is.null(cut_z) & missing(propsel)) {
propsel <- pnormmix(cut_z, mean_zr, sd_zr, mean_zf, sd_zf,
pmix_ref, lower.tail = FALSE)
}
cut_xi <- qnormmix(propsel, kappa_r, sd_xir, kappa_f, sd_xif,
pmix_ref, lower.tail = FALSE)
partit_1 <- .partit_bvnorm(cut_xi, cut_z, kappa_r, sd_xir, mean_zr, sd_zr,
cov12 = cov_z_xir)
partit_2 <- .partit_bvnorm(cut_xi, cut_z, kappa_f, sd_xif, mean_zf, sd_zf,
cov12 = cov_z_xif)
dat <- data.frame("Reference" = partit_1, "Focal" = partit_2,
row.names = c("A (true positive)", "B (false positive)",
"C (true negative)", "D (false negative)",
"Proportion selected", "Success ratio",
"Sensitivity", "Specificity"))
p <- NULL
if (plot_contour) {
x_lim <- range(c(kappa_r + c(-3, 3) * sd_xir,
kappa_f + c(-3, 3) * sd_xif))
y_lim <- range(c(mean_zr + c(-3, 3) * sd_zr,
mean_zf + c(-3, 3) * sd_zf))
contour_bvnorm(kappa_r, sd_xir, mean_zr, sd_zr, cov12 = cov_z_xir,
xlab = bquote("Latent Score" ~ (xi)),
ylab = bquote("Observed Composite" ~ (italic(Z))),
lwd = 2, col = "red", xlim = x_lim, ylim = y_lim,...)
contour_bvnorm(kappa_f, sd_xif, mean_zf, sd_zf, cov12 = cov_z_xif,
add = TRUE, lty = "dashed", lwd = 2, col = "blue", ...)
legend("topleft", c("Reference group", "Focal group"),
lty = c("solid", "dashed"), col = c("red", "blue"))
abline(h = cut_z, v = cut_xi)
x_cord <- rep(cut_xi + c(.25, -.25) * sd_xir, 2)
y_cord <- rep(cut_z + c(.25, -.25) * sd_zr, each = 2)
text(x_cord, y_cord, c("A", "B", "D", "C"))
p <- recordPlot()
}
list(propsel = propsel, cutpt_xi = cut_xi, cutpt_z = cut_z,
summary = round(dat, 3), p = p)
}


#score method of ci for dependent proportions
#b is off diagonal count
#cis other off diagonal count
#n is total sample size
#confdence coefficient = .95

scoreci<-
function (b, c, n, conf.level) 
{
    pa = 2 * n
    z = qnorm(1 - (1 - conf.level)/2)
    if (c == n) {
        ul = 1
    }
    else {
        proot = (c - b)/n
        dp = 1 - proot
        niter = 1
        while (niter <= 50) {
            dp = 0.5 * dp
            up2 = proot + dp
            pb = -b - c + (2 * n - c + b) * up2
            pc = -b * up2 * (1 - up2)
            q21 = (sqrt(pb^2 - 4 * pa * pc) - pb)/(2 * pa)
            score = (c - b - n * up2)/sqrt(n * (2 * q21 + up2 * 
                (1 - up2)))
            if (abs(score) < z) {
                proot = up2
            }
            niter = niter + 1
            if ((dp < 1e-07) || (abs(z - score) < 1e-06)) {
                niter = 51
                ul = up2
            }
        }
    }
    if (b == n) {
        ll = -1
    }
    else {
        proot = (c - b)/n
        dp = 1 + proot
        niter = 1
        while (niter <= 50) {
            dp = 0.5 * dp
            low2 = proot - dp
            pb = -b - c + (2 * n - c + b) * low2
            pc = -b * low2 * (1 - low2)
            q21 = (sqrt(pb^2 - 4 * pa * pc) - pb)/(2 * pa)
            score = (c - b - n * low2)/sqrt(n * (2 * q21 + low2 * 
                (1 - low2)))
            if (abs(score) < z) {
                proot = low2
            }
            niter = niter + 1
            if ((dp < 1e-07) || (abs(z - score) < 1e-06)) {
                ll = low2
                niter = 51
            }
        }
    }
    cint <- c(ll, ul)
    attr(cint, "conf.level") <- conf.level
    rval <- list(conf.int = cint)
    class(rval) <- "htest"
    return(rval)
}

