# Some helpful functions for twin modeling with OpenMx/umx/twinflex


# 1. maketwinvars
###############################################################################
maketwinvars <- function(var, sep = "_", sib = FALSE) {
  if (!is.character(sep)) {
    stop("The argument sep needs a character.")
  }
  if (length(sep) > 1) {
    stop("The length of the argument sep needs to be 1.")    
  }
  var1 <- paste0(var,sep,1)
  var2 <- paste0(var,sep,2)
  acevar <- paste(rep(c("A","C","E"), each = length(var)), var, sep = "")
  acevar1 <- paste(rep(c("A","C","E"), each = length(var1)), var1, sep = "")
  acevar2 <- paste(rep(c("A","C","E"), each = length(var2)), var2, sep = "")
  varswide <- c(var1,var2)
  varstotal <- c(varswide,acevar1,acevar2)
  twinvars <- list(var = var, acevar = acevar, var1 = var1, var2 = var2, acevar1 = acevar1, acevar2 = acevar2, varswide = varswide, varstotal = varstotal)
  if (sib == TRUE) {
    var3 <- paste0(var,sep,3)
    acevar <- paste(rep(c("A","C","T","E"), each = length(var)), var, sep = "")
    acevar1 <- paste(rep(c("A","C","T","E"), each = length(var1)), var1, sep = "")
    acevar2 <- paste(rep(c("A","C","T","E"), each = length(var2)), var2, sep = "")
    acevar3 <- paste(rep(c("A","C","T","E"), each = length(var3)), var2, sep = "")
    varswide <- c(var1,var2,var3)
    varstotal <- c(varswide,acevar1,acevar2,acevar3)
    twinvars <- list(var = var, acevar = acevar, var1 = var1, var2 = var2, var3 = var3, acevar1 = acevar1, acevar2 = acevar2, acevar3 = acevar3, varswide = varswide, varstotal = varstotal)
  }
  return(twinvars)
}


# 2. makesvmean: make starting values for mean vector
###############################################################################

makesvmean <- function(var, df_long, ordinal = NULL) {
  svmean <- 1
  for (i in 1:length(var)) {
    if(is.null(ordinal) | !(var[i] %in% ordinal)) {
      svmean[i] <- mean(df_long[,var[i]], na.rm = TRUE)
    }
    if(var[i] %in% ordinal) {
      svmean[i] <- 0
    }
  }
  result <- c(svmean,svmean,rep(0,2*3*length(var)))
  return(result)
}
# 3. make start values for clpm model
###############################################################################
pathclpm <- function(mat, v = 2, synchr = FALSE, horder = FALSE, ACE = FALSE) {
  mat2 <- mat
  dim <- 1:dim(mat)[1]
  col <- dim[seq(1, length(dim), v)]
  print(col)
  if (synchr == FALSE) {
    for (i in 1:dim(mat)[1]) {
      if ((i+(v-1)) <= dim(mat)[1] & i %in% col) {
        mat[1:(i+(v-1)),i:(i+(v-1))] <- 0    
        
      }
      if ((i+(v-1)) > dim(mat)[1]) {
        mat[,i] <- 0
      }
    }
  }
  if (synchr == TRUE) {
    for (i in 1:dim(mat)[1]) {
      for (j in col[2:length(col)]) {
        if (i >= j) {
          mat[1:(j-1),i] <- 0
        }
      }
    }
  }
  if (horder == FALSE) {
    for (i in 1:dim(mat)[1]) {
      if ((i+2*v) < dim(mat)[1]  & i %in% col) {
        mat[(i+2*v):(dim(mat)[1]),i:(i+(v-1))] <- 0
      }
      if ((i+2*v) == dim(mat)[1]) {
        mat[((dim(mat)[1]-1):dim(mat)[1]),i] <- 0
      }
    }
  } 
  if (ACE == TRUE) {
    diag(mat) <- diag(mat2)
  }
  return(mat)
}

# 4. label ACE matrices
###############################################################################
acelabel <- function(acevar, var, mod = FALSE, nmod = NULL, sib = FALSE) {
  acelab <- list(A = NULL,C = NULL,E = NULL)
  if (sib == TRUE) {
    acelab <- list(A = NULL,C = NULL,T = NULL, E = NULL)
  }
  nv <- length(var)
  acelab$Adim <- acevar[1:nv]
  acelab$Cdim <- acevar[(nv+1):(nv*2)]
  if (sib == TRUE) {
    acelab$Edim <- acevar[((nv*3)+1):(nv*4)]
    acelab$Tdim <- acevar[((nv*2)+1):(nv*3)]
  } else {
    acelab$Edim <- acevar[((nv*2)+1):(nv*3)]
  }
  acelab$A <- matrix(paste(rep(paste0("b_",var), each = nv), acelab$Adim, sep = "_"),nrow = nv, ncol = nv, byrow = TRUE, dimnames = list(acelab$Adim,acelab$Adim))
  acelab$C <- matrix(paste(rep(paste0("b_",var), each = nv), acelab$Cdim, sep = "_"),nrow = nv, ncol = nv, byrow = TRUE, dimnames = list(acelab$Cdim,acelab$Cdim))
  acelab$T <- matrix(paste(rep(paste0("b_",var), each = nv), acelab$Tdim, sep = "_"),nrow = nv, ncol = nv, byrow = TRUE, dimnames = list(acelab$Tdim,acelab$Tdim))
  acelab$E <- matrix(paste(rep(paste0("b_",var), each = nv), acelab$Edim, sep = "_"),nrow = nv, ncol = nv, byrow = TRUE, dimnames = list(acelab$Edim,acelab$Edim))
  
  if (mod == TRUE) {
    if (is.null(nmod)) {
      stop("When mod = TRUE, you need to specify the number of moderators with nmod!")
    }
    pathMod1 <- list(A = NULL, C = NULL, E = NULL)
    pathMod2 <- list()
    for (i in 1:nmod) {
      blab <- paste0("b",i)
      pathMod1$A <- matrix(paste(rep(paste0(blab,"_",var), each = nv), acelab$Adim, sep = "_"),nrow = nv, ncol = nv, byrow = TRUE, dimnames = list(acelab$Adim,acelab$Adim))
      pathMod1$C <- matrix(paste(rep(paste0(blab,"_",var), each = nv), acelab$Cdim, sep = "_"),nrow = nv, ncol = nv, byrow = TRUE, dimnames = list(acelab$Cdim,acelab$Cdim))
      pathMod1$E <- matrix(paste(rep(paste0(blab,"_",var), each = nv), acelab$Edim, sep = "_"),nrow = nv, ncol = nv, byrow = TRUE, dimnames = list(acelab$Edim,acelab$Edim)) 
      if (sib == TRUE) {
        pathMod1$T <- matrix(paste(rep(paste0(blab,"_",var), each = nv), acelab$Tdim, sep = "_"),nrow = nv, ncol = nv, byrow = TRUE, dimnames = list(acelab$Tdim,acelab$Tdim)) 
        
      }
      name <- paste0("Mod",i)
      pathMod2[[name]] <- pathMod1
    }
    acelab$Mod <- pathMod2
  }
  acelab$Acov <- matrix(paste(rep(paste0("cov_",acelab$Adim), each = nv), acelab$Adim, sep = "_"),nrow = nv, ncol = nv, byrow = TRUE, dimnames = list(acelab$Adim,acelab$Adim))
  acelab$Acov[upper.tri(acelab$Acov)] <- t(acelab$Acov)[upper.tri(acelab$Acov)]
  acelab$Ccov <- matrix(paste(rep(paste0("cov_",acelab$Cdim), each = nv), acelab$Cdim, sep = "_"),nrow = nv, ncol = nv, byrow = TRUE, dimnames = list(acelab$Cdim,acelab$Cdim))
  acelab$Ccov[upper.tri(acelab$Ccov)] <- t(acelab$Ccov)[upper.tri(acelab$Ccov)]
  acelab$Ecov <- matrix(paste(rep(paste0("cov_",acelab$Edim), each = nv), acelab$Edim, sep = "_"),nrow = nv, ncol = nv, byrow = TRUE, dimnames = list(acelab$Edim,acelab$Edim))
  acelab$Ecov[upper.tri(acelab$Ecov)] <- t(acelab$Ecov)[upper.tri(acelab$Ecov)]
  if (sib == TRUE) {
    acelab$Tcov <- matrix(paste(rep(paste0("cov_",acelab$Tdim), each = nv), acelab$Tdim, sep = "_"),nrow = nv, ncol = nv, byrow = TRUE, dimnames = list(acelab$Tdim,acelab$Tdim))
    acelab$Tcov[upper.tri(acelab$Tcov)] <- t(acelab$Tcov)[upper.tri(acelab$Tcov)]
  }
  return(acelab)
}


# 5. make A matrix
###############################################################################
makeAmat <- function(b, pZ, pCov, a, c, e, Bottom,dimname) {
  values <- rbind(cbind(b$values,pZ$values,pCov$values,a$values,c$values,e$values,pZ$values,pZ$values,pZ$values),
                  cbind(pZ$values,b$values,pCov$values,pZ$values,pZ$values,pZ$values,a$values,c$values,e$values),
                  Bottom$values)
  
  dimnames(values) <- list(dimname, dimname)
  values <- as.data.frame(values)
  labels <- rbind(cbind(b$labels,pZ$labels,pCov$labels,a$labels,c$labels,e$labels,pZ$labels,pZ$labels,pZ$labels),
                  cbind(pZ$labels,b$labels,pCov$labels,pZ$labels,pZ$labels,pZ$labels,a$labels,c$labels,e$labels),
                  Bottom$labels)
  dimnames(labels) <- list(dimname, dimname)
  labels <- as.data.frame(labels)
  free <- rbind(cbind(b$free,pZ$free,pCov$free,a$free,c$free,e$free,pZ$free,pZ$free,pZ$free),
                cbind(pZ$free,b$free,pCov$free,pZ$free,pZ$free,pZ$free,a$free,c$free,e$free),
                Bottom$free)
  dimnames(free) <- list(dimname, dimname)
  free <- as.data.frame(free)
  matrixA <- list(values = values, labels = labels, free = free)
  return(matrixA)
}

# 6. make S matrix
makeSmat <- function(cMan, cZ, cACEM, cA, cC, cE, dimname) {
  mzvalues <- cbind(
    rbind(cbind(cMan$values,cZ$values),
          cbind(cZ$values,cMan$values),
          t(cACEM$values)),
    rbind(cACEM$values,
          cbind(cA$values,cZ$values,cZ$values,cA$values,cZ$values,cZ$values),
          cbind(cZ$values,cC$values,cZ$values,cZ$values,cC$values,cZ$values),
          cbind(cZ$values,cZ$values,cE$values,cZ$values,cZ$values,cZ$values),
          cbind(cA$values,cZ$values,cZ$values,cA$values,cZ$values,cZ$values),
          cbind(cZ$values,cC$values,cZ$values,cZ$values,cC$values,cZ$values),
          cbind(cZ$values,cZ$values,cZ$values,cZ$values,cZ$values,cE$values)))
  dimnames(mzvalues) <- list(dimname, dimname)
  Alab <- matrix(paste0(".5*",cA$labels), ncol = dim(cA$labels)[1], nrow = dim(cA$labels)[1], dimnames = dimnames(cA$labels))
  
  mzlabels <- cbind(
    rbind(cbind(cMan$labels,cZ$labels),
          cbind(cZ$labels,cMan$labels),
          t(cACEM$labels)),
    rbind(cACEM$labels,
          cbind(cA$labels,cZ$labels,cZ$labels,cA$labels,cZ$labels,cZ$labels),
          cbind(cZ$labels,cC$labels,cZ$labels,cZ$labels,cC$labels,cZ$labels),
          cbind(cZ$labels,cZ$labels,cE$labels,cZ$labels,cZ$labels,cZ$labels),
          cbind(cA$labels,cZ$labels,cZ$labels,cA$labels,cZ$labels,cZ$labels),
          cbind(cZ$labels,cC$labels,cZ$labels,cZ$labels,cC$labels,cZ$labels),
          cbind(cZ$labels,cZ$labels,cZ$labels,cZ$labels,cZ$labels,cE$labels)))
  dimnames(mzlabels) <- list(dimname, dimname)
  mzfree <- cbind(
    rbind(cbind(cMan$free,cZ$free),
          cbind(cZ$free,cMan$free),
          t(cACEM$free)),
    rbind(cACEM$free,
          cbind(cA$free,cZ$free,cZ$free,cA$free,cZ$free,cZ$free),
          cbind(cZ$free,cC$free,cZ$free,cZ$free,cC$free,cZ$free),
          cbind(cZ$free,cZ$free,cE$free,cZ$free,cZ$free,cZ$free),
          cbind(cA$free,cZ$free,cZ$free,cA$free,cZ$free,cZ$free),
          cbind(cZ$free,cC$free,cZ$free,cZ$free,cC$free,cZ$free),
          cbind(cZ$free,cZ$free,cZ$free,cZ$free,cZ$free,cE$free)))
  dimnames(mzfree) <- list(dimname, dimname)
  
  dzvalues <- cbind(
    rbind(cbind(cMan$values,cZ$values),
          cbind(cZ$values,cMan$values),
          t(cACEM$values)),
    rbind(cACEM$values,
          cbind(cA$values,cZ$values,cZ$values,0.5*cA$values,cZ$values,cZ$values),
          cbind(cZ$values,cC$values,cZ$values,cZ$values,cC$values,cZ$values),
          cbind(cZ$values,cZ$values,cE$values,cZ$values,cZ$values,cZ$values),
          cbind(0.5*cA$values,cZ$values,cZ$values,cA$values,cZ$values,cZ$values),
          cbind(cZ$values,cC$values,cZ$values,cZ$values,cC$values,cZ$values),
          cbind(cZ$values,cZ$values,cZ$values,cZ$values,cZ$values,cE$values)))
  dimnames(dzvalues) <- list(dimname, dimname)
  dzlabels <- cbind(
    rbind(cbind(cMan$labels,cZ$labels),
          cbind(cZ$labels,cMan$labels),
          t(cACEM$labels)),
    rbind(cACEM$labels,
          cbind(cA$labels,cZ$labels,cZ$labels,Alab,cZ$labels,cZ$labels),
          cbind(cZ$labels,cC$labels,cZ$labels,cZ$labels,cC$labels,cZ$labels),
          cbind(cZ$labels,cZ$labels,cE$labels,cZ$labels,cZ$labels,cZ$labels),
          cbind(Alab,cZ$labels,cZ$labels,cA$labels,cZ$labels,cZ$labels),
          cbind(cZ$labels,cC$labels,cZ$labels,cZ$labels,cC$labels,cZ$labels),
          cbind(cZ$labels,cZ$labels,cZ$labels,cZ$labels,cZ$labels,cE$labels)))
  dimnames(dzlabels) <- list(dimname, dimname)
  dzfree <- cbind(
    rbind(cbind(cMan$free,cZ$free),
          cbind(cZ$free,cMan$free),
          t(cACEM$free)),
    rbind(cACEM$free,
          cbind(cA$free,cZ$free,cZ$free,cA$free,cZ$free,cZ$free),
          cbind(cZ$free,cC$free,cZ$free,cZ$free,cC$free,cZ$free),
          cbind(cZ$free,cZ$free,cE$free,cZ$free,cZ$free,cZ$free),
          cbind(cA$free,cZ$free,cZ$free,cA$free,cZ$free,cZ$free),
          cbind(cZ$free,cC$free,cZ$free,cZ$free,cC$free,cZ$free),
          cbind(cZ$free,cZ$free,cZ$free,cZ$free,cZ$free,cE$free)))
  dimnames(dzfree) <- list(dimname, dimname)
  matrixS <- list(MZvalues = as.data.frame(mzvalues), MZlabels = as.data.frame(mzlabels), MZfree = as.data.frame(mzfree),
                  DZvalues = as.data.frame(dzvalues), DZlabels = as.data.frame(dzlabels), DZfree = as.data.frame(dzfree))
  
  
  
}

# 7. calculate expected covariance matrix

makeExpCov <- function(A, S, m, l) {
  library(Ryacas)
  inputA <- A[["labels"]]
  inputSMZ <- S[["MZlabels"]]
  inputSDZ <- S[["DZlabels"]]
  Fil1 <- diag(m)
  Fil0 <- matrix(0,m,l)
  Fil <- cbind(Fil1,Fil0)
  Fil <- ysym(Fil)
  Iden <- diag((m+l))
  Iden <- ysym(Iden)
  
  # A
  for (j in colnames(A[["labels"]])) {
    for (i in rownames(A[["labels"]])) {
      if (A[["free"]][i,j] == FALSE) {
        inputA[i,j] <- A[["values"]][i,j]
      }
      
    }
  }
  inputA <- as.matrix(inputA)
  dimnames(inputA) <- NULL
  inputA <- ysym(inputA)
  # S (MZ)
  for (j in colnames(S[["MZlabels"]])) {
    for (i in rownames(S[["MZlabels"]])) {
      if (S[["MZfree"]][i,j] == FALSE) {
        inputSMZ[i,j] <- S[["MZvalues"]][i,j]
      }
      
    }
  }
  inputSMZ <- as.matrix(inputSMZ)
  dimnames(inputSMZ) <- NULL
  inputSMZ <- ysym(inputSMZ)
  # S (DZ)
  for (j in colnames(S[["DZlabels"]])) {
    for (i in rownames(S[["DZlabels"]])) {
      if (S[["DZfree"]][i,j] == FALSE) {
        inputSDZ[i,j] <- S[["DZvalues"]][i,j]
      }
      
    }
  }
  inputSDZ <- as.matrix(inputSDZ)
  dimnames(inputSDZ) <- NULL
  inputSDZ <- ysym(inputSDZ)
  
  # compute (I-A)^-1
  ImA <- Iden - inputA
  eCMZ <- Fil%*%solve(Iden-inputA)%*%inputSMZ%*%t(solve(Iden-inputA))%*%t(Fil)
  eCDZ <- Fil%*%solve(Iden-inputA)%*%inputSDZ%*%t(solve(Iden-inputA))%*%t(Fil)
  result <- list(MZ = eCMZ, DZ = eCDZ)
  return(result)
}


# 7. Output
reduceMyModel <- function(fit, nodrop) {
  #fit <- fitACE
  #
  outputlist <- list()
  outputlist[["Baseline"]] <- summary(fit)
  output <- summary(fit)[["parameters"]]
  output$z <- summary(fit)[["parameters"]]$Estimate/summary(fit)[["parameters"]]$Std.Error
  output$p <- 2*pnorm(-abs(output$z))
  output <- output[,c("name","p")]
  output <- output[order(-output$p),]
  continue <- TRUE
  while (continue == TRUE) {
    drop <- output[1,c("name","p")]
    pval <- drop[,"p"]
    drop <- drop[,"name"]
    if (pval > 0.05) {
      output <- output[!(output$name == drop),]
      if (!(drop %in% nodrop)) {
        cat("Now I drop the following parameter: ", drop, "\n It's p-value is: ", round(pval,2), "\n\n")
        name = paste0("Drop '", drop,"'")
        fitold <- fit  
        model1   <- mxModel(fitold, name = name)
        model1   <- omxSetParameters(model1, 
                                     labels= drop, 
                                     free=FALSE, values=0)
        fit    <- mxTryHard(model1)
        print(nrow(summary(fit)[["parameters"]]))
        comp1 <- mxCompare(fitold,fit)
        if (comp1@results[["p"]][2] >= 0.05) {
          outputlist[[name]] <- fit
        }
        else {
          cat("Model fitting stopped. Dropping the parameter", drop, "resulted in a significant deterioration of model fit. If there are other non-significant parameter, I will try to drop them.\n\n")
          fit <- fitold 
        }
      }
      else {
        continue <- TRUE
      }
    }
    else {
      cat("Model fitting stopped. There are only significant (p < 0.05) parameters in the model.\n\n")
      break
    }
  }
  return(outputlist)
}

makeOutput <- function(sumACE, rounded = FALSE, n = 2) {
  # Generate Output
  sig.level <- function(coef,std.error){
    ifelse(abs(coef/std.error)>= qnorm(p = 0.001/2, lower.tail = FALSE),"$^{***}$",
           ifelse(abs(coef/std.error)>= qnorm(p = 0.01/2, lower.tail = FALSE),"$^{**}$",
                  ifelse(abs(coef/std.error) >= qnorm(p = 0.05/2, lower.tail = FALSE), "$^{*}$",
                         ifelse(abs(coef/std.error) >= qnorm(p = 0.1/2, lower.tail = FALSE), "$^{+}$","")
                  )
           )
    )
  }
  
  parlist <- sumACE[["parameters"]]
  zval <- parlist$Estimate/parlist$Std.Error
  pval <- 2*pnorm(-abs(zval))
  star <- sig.level(parlist$Estimate, parlist$Std.Error)
  Par <- parlist[1]
  Est <- parlist[5]
  SE <- parlist[6]
  if (rounded == TRUE) {
    zval <- round(zval, n)
    pval <- round(pval, n)
    Est <- round(Est, n)
    SE <- round(SE, n)
  }
  parlist <- cbind(Par, Est, SE, zval, pval, star)
  
  colnames(parlist) <- c("Par","Est","SE","z","p","")
  return(parlist)
}

# make Thresholds

makeThresholds <- function(nt, vars, data) {
  nlevels <- function(variable) { 
    nolevels <- sort(unique(data[,variable]))
  }
  levelslist <- lapply(ordinalwide, nlevels)
  data[,ordinalwide] <- mxFactor(data[,ordinalwide], levels = levelslist)
  ordinallength <- unlist(lapply(levelslist,length))
  nTh       <- ordinallength-1 # No of thresholds
  ntvo      <- length(ordinalwide) # Total No of ordinal vars
  freeThresholds <- function(nt) { # assumes mxMatrix(byrow = FALSE)
    values <- NULL
    if(nt == 1) {
      values <- c(TRUE,rep(FALSE,(max(nTh)-1))) 
    }
    else if (nt > 1 & nt < max(nTh)) {
      values <- c(rep(FALSE,2),rep(TRUE,(nt-2)),rep(FALSE,(max(nTh)-nt)))        
    }
    else if (nt > 1 & nt == max(nTh)) {
      values <- c(rep(FALSE,2),rep(TRUE,(nt-2))) 
    }
  }
  valThresholds <- function(nt) { # assumes mxMatrix(byrow = FALSE)
    values <- NULL
    if(nt == 1) {
      values <- c(.5,rep(0,(max(nTh)-1))) # unused thresholds set to zero (check: post #20 by T. Bates in: https://openmx.ssri.psu.edu/node/4538)
    }
    else if (nt > 1 & nt < max(nTh)) {
      values <- c(0,rep(1,(nt-1)),rep(0,(max(nTh)-nt)))           
    }
    else if (nt > 1 & nt == max(nTh)) {
      values <- c(0,rep(1,(nt-1))) 
    }
  }
  
  lbThresholds <- function(nt) { # assumes mxMatrix(byrow = FALSE)
    lb <- NULL
    if(nt == 1) {
      lb <- rep(NA,max(nTh))
    }
    else if (nt > 1 & nt < max(nTh)) {
      lb <- c(-3,rep(.0001,(nt-1)),rep(NA,(max(nTh)-nt)))         
    }
    else if (nt > 1 & nt == max(nTh)) {
      lb <- c(-3,rep(.0001,(nt-1)))  
    }
  }
  
  labelThresholds <- function(vars) { # assumes mxMatrix(byrow = FALSE)
    labels <- NULL
    noTh <- (length(levels(data[,vars])))-1
    varsnew <- unlist(sapply(strsplit(vars, split=sep, fixed=TRUE), function(x) (x[1])))
    if (noTh == 1) { # binary
      labels <- c(paste0("th",varsnew,1),rep(NA,(max(nTh)-1)))
    }
    else if (noTh > 1 & noTh < max(nTh)) { # ordinal
      labels <- paste0("th",varsnew,(1:noTh),rep(NA,(max(nTh)-noTh)))
    }
    else if (noTh > 1 & noTh == max(nTh)) { # ordinal
      labels <- paste0("th",varsnew,(1:noTh))
    }
  }
  frTh <- unlist(lapply(nTh, freeThresholds))
  svTh <- unlist(lapply(nTh, valThresholds))
  lbTh <- unlist(lapply(nTh, lbThresholds))
  labTh <- unlist(lapply(vars, labelThresholds))
  print(svTh)
  # Thresholds definieren
  thinG     <- mxMatrix(type="Full", nrow=max(nTh), ncol=ntvo, free=frTh, byrow = FALSE, values=svTh, lbound=lbTh, labels=labTh, name="thinG") # matrix of threshold increments
  inc       <- mxMatrix(type="Lower", nrow=max(nTh), ncol=max(nTh), free=FALSE, values=1, name="inc") # matrix of lower 1
  threG     <- mxAlgebra(expression= inc %*% thinG, name="threG")
  return(list(ThMat = threG, thinMat = thinG, incMat = inc, numberTh = nTh, free = frTh, starting = svTh, lbound = lbTh, label = labTh))
}


exph <- function(fit, posdv, tryHard = FALSE) {
if (is.null(fit)) {
  stop("This function needs a fit object to work properly")
}
  if (is.null(posdv)) {
    stop("I need to know the position of your outcome to calculate the heritabilities")
  }
  cat(paste0("Hello, you ordered the heritability estimates for fit object '", deparse(substitute(fit)),"'.\n\n"))
  cat(paste0("You gave me the following information: \n\n Your outcome is in position ",deparse(substitute(posdv)),
             "\n\n So there are ",posdv-1, " predictors which might account for the heritability of you outcome\n\n"))
  cat("You will get two output object stored in a list: \n1) A fit object with the heritabilities as model parameters. \n2) A data frame with the heritabilities.\n\n")
  
  indmat <- mxMatrix(type = "Full", nrow = 1, ncol = posdv, values = c(rep(0,posdv-1),1), name = "IndMat")
  Aeffmat <- mxAlgebra(expression = IndMat%*%(top.a^2), name = "AeffMat")
  Yvar <- mxAlgebra(expression = IndMat%*%top.ACE%*%t(IndMat), name = "VarY")
  Avar <- mxAlgebra(expression = IndMat%*%top.A%*%t(IndMat), name = "VarA")
  h2 <- mxAlgebra(expression = VarA/VarY, name = "heritability")
  h2source <- mxAlgebra(expression = AeffMat/VarY, name = "heritabilitysource")
  model <- mxModel(fit,indmat,Aeffmat,Yvar,Avar,h2,h2source)
  if (tryHard == FALSE) {
    fit <- mxRun(model)
  } else {
    fit <- mxTryHard(model)  
  }
  h2est <- fit@algebras[["heritability"]]$result
  h2sourceest <- fit@algebras[["heritabilitysource"]]$result
  
  SEh2 <- mxSE("heritability", fit)
  SEh2source <- mxSE("heritabilitysource", fit)
  
  h2df <- data.frame(c(h2est, h2sourceest),c(SEh2, SEh2source))
  colnames(h2df) <- c("Estimate","Standard Error")
  rownames(h2df) <- c("Overall h2", paste("Due to genes associated with Var",seq(1:posdv)))
  
  result <- list(fit, h2df)
  names(result) <- c("Fit Object", "Decomposed h2")
  return(result)
}