R/chmethod.R

In sgaure/lfe: Linear Group Fixed Effects

# $Id: chmethod.R 1693 2015-04-07 09:36:29Z sgaure $
findfe <- function(dd,Rhs,se=FALSE) {
  # find references
  refnames <- attr(dd,'refnames')
  nm <- c(colnames(dd),refnames)
  refcnt <- attr(dd,'refcnt')
  # add in the reference in dg

  dg <- c(diag(dd),refcnt)
  ok <- 1:ncol(dd)

  if(se) sev <- double(length(nm))
  alphacoef <- double(length(nm))

  # the super-nodal algorithm
  # is default and far better, but it consumes more memory

  trysolve <- try(solve(dd,Rhs))
  if(inherits(trysolve,'try-error')) {
    if(grepl('problem too large',geterrmessage())) {
      message(paste('Never mind, trying *slower* non-supernodal algorithm, nnz=',nnzero(dd)))
      message(paste(date(),'This may be an opportunity for a nice cup of tea. Or two.'))
      gc()
      ch <- Cholesky(dd,super=FALSE,perm=TRUE)
      trysolve <- solve(ch,Rhs)
      rm(ch)
      gc()
    } else {
      stop(geterrmessage())
    }
  } 


  alphacoef[ok] <- as.vector(trysolve)
  if(se) {
    # is there a faster way to find the diagonal of the inverse?
    sev[ok] <- sqrt(diag(solve(dd)))*attr(se,'sefactor')
    alpha <- data.frame(effect=alphacoef,se=sev,obs=dg)
  } else {
    alpha <- data.frame(effect=alphacoef,obs=dg)
  }
  rownames(alpha) <- nm
  alpha
}


#makedummies <- function(factors) {
#  nm <- c()
#  dummies <- Matrix(0,0,length(factors[[1]]))
#  for(i in 1:length(factors)) {
#    f <- factors[[i]]
#    dummies <- rBind(dummies,as(f,'sparseMatrix'))
#    nm <- c(nm,paste(names(factors)[[i]],levels(f),sep='.'))
#  }
#  rownames(dummies) <- nm
#  dummies
#}

makedd.full <- function(factors) {
#  dm <- makedummies(factors)
  dm <- t(makeDmatrix(factors))
  nm <- rownames(dm)
  dd <- tcrossprod(dm)
  rownames(dd) <- colnames(dd) <- nm
  attr(dd,'dummies') <- dm
  attr(dd,'nm') <- nm
  dd
}

makeddlist <- function(factors) {
  if(length(factors) > 2) {
    if(is.null(attr(factors,'references'))) {
      # find references by fiddling with Cholesky
      message('*** More than two groups, finding refs by Cholesky pivots, interpret at own risk')
      # first the full matrix, find small pivots
      dd <- makedd.full(factors)
      orignm <- attr(dd,'nm')

      # add small amount to diagonal
      eps <- sqrt(.Machine$double.eps)
      Ch <- try(Cholesky(dd,super=TRUE,perm=TRUE,Imult=eps))
      if(inherits(Ch,'try-error') && grepl('problem too large',geterrmessage())) {
        Ch <- Cholesky(dd,super=FALSE,perm=TRUE,Imult=eps)
      }
      # strangely enough, coercing to sparseMatrix doesn't take care of
      # the permutation, we apply it manually.  Let's hope it's never fixed.
      rm(dd); gc()
      pivot <- Ch@perm
      ch <- as(Ch,'sparseMatrix')
      rm(Ch); gc()
      dg <- diag(ch)[order(pivot)]**2
      rm(ch); gc()
      refs <- (dg < eps**(1/3))
      refnames <- orignm[refs]
      message(paste('***',length(refnames),'references found'))
    } else {
      refnames <- attr(factors,'references')
      orignm <- unlist(lapply(names(factors),
                              function(n) paste(n,levels(factors[[n]]),sep='.')))
    }
    # there may be references in more than one factor
    # remove all of them
    # create factor list with named levels
    nf <- lapply(names(factors),function(n) {
      f <- factors[[n]]
      levels(f) <- paste(n,levels(f),sep='.')
      f
    })
    # remove reference levels, and remove the prefix
    # find the levels
    lev <- lapply(nf,function(f) which(levels(f) %in% refnames))
    nnf <- mapply(function(f,l) factor(f,exclude=levels(f)[l]),factors,lev,SIMPLIFY=FALSE)
    dd <- makedd.full(nnf)
    attr(dd,'keep') <- 1:length(nnf[[1]])
    attr(dd,'refnames') <- refnames
#    attr(dd,'refcnt') <- rep(1,length(refnames))
# find the number of occurences
    cntlst <- unlist(lapply(refnames,function(n) lapply(nf,function(f) sum(f == n))))
    attr(dd,'refcnt') <- cntlst[cntlst > 0]
    attr(dd,'comp') <- 1
    res <- list(dd)
    attr(res,'nm') <- orignm
  } else {
    # 2 or fewer factors, find references by component
    cf <- compfactor(factors)
    nml <- lapply(factors,function(f) levels(f))
    nm <- unlist(lapply(names(nml),function(n) paste(n,nml[[n]],sep='.')))
    res <- list()
    li <- 1
# this loop suffers from too much copying and stuff
# when there are many components (e.g. like 10000)
    remfact <- factors
    fullidx <- 1:length(factors[[1]])
    for(l in levels(cf)) {
      # find those in this level
      keep <- which(cf == l)
#      cat(date(),'comp fact',li,'size',length(keep),'\n')
      fcomp <- lapply(remfact,function(f) factor(f[keep]))
      remfact <- lapply(remfact,function(f) factor(f[-keep]))
      cf <- factor(cf[-keep])
      # then the reference level
      maxrefs <- lapply(fcomp,function(f) {tf <- table(f); m <- which.max(tf); tf[m]})
      # in which factor
      rfac <- which.max(unlist(maxrefs))
      # which level
      reflevel <- names(maxrefs[[rfac]])
      # drop that level from the factor
      fcomp[[rfac]] <- factor(fcomp[[rfac]],exclude=reflevel)
      refname <- paste(names(remfact)[[rfac]],reflevel,sep='.')
      # remove those without levels
      len <- unlist(lapply(fcomp,nlevels))
      fcomp <- fcomp[len > 0]
      dd <- makedd.full(fcomp)
      # the keep attribute should be relative to
      # the full factor, not to remfact
      attr(dd,'keep') <- fullidx[keep]
      fullidx <- fullidx[-keep]
      attr(dd,'refnames') <- refname
      attr(dd,'refcnt') <- max(unlist(maxrefs))
      attr(dd,'comp') <- li
      res[[li]] <- dd
      li <- li+1
#      res <- c(res, list(dd))
    }
    attr(res,'nm') <- nm
  }
  res
}