
outputReg_pw <- function(regModels, varNames=NULL, varLabels=NULL, addVIFs=FALSE){
  varList <- as.character()
  numCols <- length(regModels)
  #Find longest set of variable listed
  for(i in 1:numCols){
    m <- regModels[[i]]
    vars <- names(m$coefficients)
    coefNames <- names(m$coefficients)
    mvars <- coefNames[which(grepl(":", coefNames))]
    rVars <- vars[which(substr(vars, 1, 4) == "fact" | substr(vars, 1, 2)=="ID")]  #remove list
    rVars <- rVars[which(!(rVars %in% mvars))]  ##removing from the remove list
    vars <- vars[which(!(vars %in% rVars))]

    varList <- unique(c(varList, vars))
  }
  
 # rowOrder <- varNames
  #  rowLab <- varLabels
  
  diagLabs <-  c("N", "Countries", "R-Squared")
  
  numR <- length(varList) + length(diagLabs)
  regdf <- as.data.frame(matrix(NA, ncol = 1, nrow = numR))
  row.names(regdf) <- c(varList, diagLabs)
  names(regdf)[1] <- "vars"
  regdf$vars <- row.names(regdf)
  
  for(i in 1:numCols){
    m <- regModels[[i]]
    vars <- names(m$coefficients)
    coefNames <- names(m$coefficients)
    mvars <- coefNames[which(grepl(":", coefNames))]
    rVars <- vars[which(substr(vars, 1, 4) == "fact")]  #remove list
    rVars <- rVars[which(!(rVars %in% mvars))]  ##removing from the remove list
    vars <- vars[which(!(vars %in% rVars))]
   
    coefs <- m$coefficients[which(names(m$coefficients)%in%vars)]
    coefs <- round(coefs, 3)
    pvals <- m$pvals[which(names(m$coefficients)%in%vars)]
    pcses <- round(m$pcses[which(names(m$coefficients)%in%vars)],3)
    
    # get the stars
    coefs <- as.character(coefs)
    for(j in 1:length(vars)){
      p <- pvals[j]
      s <- ifelse(p < .001, "***", 
                  ifelse(p < .01, "**", 
                         ifelse(p < .05, "*", "")))
      coefs[j] <- paste(coefs[j], s, " (",  pcses[j], ")", sep="")
    }

    df <- as.data.frame(coefs)
    df$vars <- vars
    row.names(df) <- df$vars
    mName <- paste("model",i,sep="")
    
    names(df)[1] <- "c1"
    df <- df[,c("vars", "c1")]

    ##add diagnostics
      diagDF <- as.data.frame(matrix(NA, ncol = 1, nrow = length(diagLabs)))
      row.names(diagDF) <- diagLabs
      names(diagDF)[1] <- "c1"
      
      pan <- m$panelStructure
      
      Diag_obs <- sum(pan$obs.mat) #"Total obs."  
      Diag_obs_c <- nrow(pan$obs.mat) # "Countries"
      Diag_obs_y <- ncol(pan$obs.mat)  # "Years"
      Diag_obs_av <- round(mean(apply(pan$obs.mat, 2, sum)))  #"Avg obs. per panel"
      Diag_obs_r  <- round(m$r2,4)
      
      diagDF$c1[1] <- Diag_obs  #"Total obs."  
      diagDF$c1[2] <- Diag_obs_c  
      diagDF$c1[3] <- Diag_obs_r
      
      diagDF$vars <- row.names(diagDF)
      diagDF$c1 <- as.character(diagDF$c1)
      df_coefs <- dplyr::bind_rows(df, diagDF)
      regdf <- dplyr::left_join(regdf, df_coefs, by="vars")
    }
  #row.names(regdf) <- c("Intercept", rowLab, diagLabs)
  #regdf <- regdf[,-1]
  nc <- ncol(regdf)
  
    if(length(which(regdf$vars=="MIL:YR"))>0 & length(which(regdf$vars=="YR:MIL"))>0){
      for(k in 2:nc){
        if(is.na(regdf[which(regdf$vars=="MIL:YR"),k])){
        regdf[which(regdf$vars=="MIL:YR"),k] <- regdf[which(regdf$vars=="YR:MIL"),k]
        }}
        regdf <- regdf[-which(regdf$vars=="YR:MIL"),]
        }
    
if(length(which(regdf$vars=="MIL:incomeLow Income"))>0 & length(which(regdf$vars=="incomeLow Income:MIL"))>0){
  for(k in 2:nc){
    if(is.na(regdf[which(regdf$vars=="MIL:incomeLow Income"),k])){
      regdf[which(regdf$vars=="MIL:incomeLow Income"),k] <- regdf[which(regdf$vars=="incomeLow Income:MIL"),k]
    }}
  regdf <- regdf[-which(regdf$vars=="incomeLow Income:MIL"),]
  }

  if(addVIFs==TRUE){
    source("myFunctions/calc_vifs.R")
    
    NC <- ncol(regdf)
    for(h in 2:NC){
      m <- regModels[[h-1]]
      vfResult <- calc_vifs(m)
      
      for(v in 1:nrow(regdf)){
        vrN <- regdf$vars[v]
        vifN <- which(names(vfResult)==vrN)
        if(length(vifN)>0){
          vf <- round(vfResult[[vifN]],2)
          regdf[v,h] <- paste(regdf[v,h], " [", vf, "]", sep="")
        }
      }}}
    
  return(regdf)
}