Nothing
# $Id: felm.old.R 1655 2015-03-18 18:51:06Z sgaure $
felm.old <- function(formula, fl, data) {
mf <- match.call(expand.dots = FALSE)
if (missing(fl)) {
# we should rather parse the formula tree
# find the terms involving G
trm <- terms(formula, special = "G")
feidx <- attr(trm, "specials")$G - 1
festr <- paste(labels(trm)[feidx], collapse = "+")
if (festr == "") stop("No factors specified")
# remove the G-terms from formula
formula <- update(formula, paste(". ~ . -(", festr, ")"))
mf[["formula"]] <- formula
# then make a list of them, and find their names
felist <- parse(text = paste("list(", gsub("+", ",", festr, fixed = TRUE), ")", sep = ""))
nm <- eval(felist, list(G = function(t) as.character(substitute(t))))
# collapse them in case there's an interaction with a funny name
nm <- lapply(nm, paste, collapse = ".")
# replace G with as.factor, eval with this, and the parent frame, or with data
# allow interaction factors with '*'
iact <- function(a, b) interaction(a, b, drop = TRUE)
if (missing(data)) {
fl <- eval(felist, list(G = as.factor, "*" = iact))
} else {
G <- as.factor
fl <- local({
"*" <- iact
eval(felist, data, environment())
})
}
gc()
names(fl) <- nm
} else {
# warning('The fl-argument is obsolete')
}
if (!is.list(fl)) stop("need at least one factor")
fl <- lapply(fl, as.factor)
if (is.null(names(fl))) names(fl) <- paste("fe", 1:length(fl), sep = "")
# mf <- match.call(expand.dots = FALSE)
m <- match(c("formula", "data"), names(mf), 0L)
mf <- mf[c(1L, m)]
mf$drop.unused.levels <- TRUE
mf[[1L]] <- as.name("model.frame")
mf <- eval(mf, parent.frame())
mt <- attr(mf, "terms")
y <- model.response(mf, "numeric")
# try a sparse model matrix to save memory when removing intercept
# though, demeanlist must be full. Ah, no, not much to save because
# it won't be sparse after centering
# we should rather let demeanlist remove the intercept, this
# will save memory by not copying. But we need to remove it below in x %*% beta
# (or should we extend beta with a zero at the right place, it's only
# a vector, eh, is it, do we not allow matrix lhs? No.)
x <- model.matrix(mt, mf)
rm(mf)
icpt <- 0
icpt <- which(attr(x, "assign") == 0)
if (length(icpt) == 0) icpt <- 0
ncov <- ncol(x) - (icpt > 0)
if (ncov == 0) {
# No covariates
fr <- demeanlist(y, fl)
z <- list(r.residuals = y, fe = fl, p = 0, cfactor = compfactor(fl), residuals = fr, call = match.call())
class(z) <- "felm"
return(z)
}
# here we need to demean things
dm <- demeanlist(list(y = y, x = x), fl, icpt)
yz <- dm[[1]]
xz <- dm[[2]]
rm(dm)
gc()
badconv <- attr(xz, "badconv") + attr(yz, "badconv")
dim(xz) <- c(nrow(x), ncov)
attributes(yz) <- attributes(y)
# here we just do an lm.fit, however lm.fit is quite slow since
# it doesn't use blas (in particular it can't use e.g. threaded blas in acml)
# so we have rolled our own.
# we really don't return an 'lm' object or other similar stuff, so
# we should consider using more elementary operations which map to blas-3
# eg. solve(crossprod(xz),t(xz) %*% yz)
# Or, even invert by solve(crossprod(xz)) since we need
# the diagonal for standard errors. We could use the cholesky inversion
# chol2inv(chol(crossprod(xz)))
cp <- crossprod(xz)
ch <- cholx(cp)
# ch <- chol(cp)
# beta <- drop(inv %*% (t(xz) %*% yz))
# remove multicollinearities
badvars <- attr(ch, "badvars")
b <- crossprod(xz, yz)
if (is.null(badvars)) {
beta <- as.vector(backsolve(ch, backsolve(ch, b, transpose = TRUE)))
inv <- chol2inv(ch)
} else {
beta <- rep(NaN, nrow(cp))
beta[-badvars] <- backsolve(ch, backsolve(ch, b[-badvars], transpose = TRUE))
inv <- matrix(NaN, nrow(cp), ncol(cp))
inv[-badvars, -badvars] <- chol2inv(ch)
}
rm(b)
if (icpt > 0) names(beta) <- colnames(x)[-icpt] else names(beta) <- colnames(x)
# cat(date(),'projected system finished\n')
z <- list(coefficients = beta, badconv = badconv)
N <- nrow(xz)
p <- ncol(xz) - length(badvars)
# how well would we fit with all the dummies?
# the residuals of the centered model equals the residuals
# of the full model, thus we may compute the fitted values
# resulting from the full model.
zfit <- xz %*% ifelse(is.na(beta), 0, beta)
rm(xz)
zresid <- yz - zfit
rm(yz)
z$fitted.values <- y - zresid
z$residuals <- zresid
# insert a zero at the intercept position
if (length(fl) > 0) {
if (icpt > 0) ibeta <- append(beta, 0, after = icpt - 1) else ibeta <- beta
pred <- x %*% ifelse(is.na(ibeta), 0, ibeta)
z$r.residuals <- y - pred
} else {
z$r.residuals <- zresid
}
# z$xb <- pred
rm(x)
rm(y)
gc()
z$cfactor <- compfactor(fl)
numrefs <- nlevels(z$cfactor) + max(length(fl) - 2, 0)
numdum <- sum(unlist(lapply(fl, nlevels))) - numrefs
z$numrefs <- numrefs
# if(length(fl) <= 2) {
# numdum <- sum(unlist(lapply(fl,nlevels))) - nlevels(z$cfactor)
# } else {
# numdum <- sum(unlist(lapply(fl,nlevels))) - length(fl) + 1
# }
z$df <- N - p - numdum
vcvfactor <- sum(z$residuals**2) / z$df
z$vcv <- inv * vcvfactor
z$se <- sqrt(diag(z$vcv))
z$sefactor <- sqrt(vcvfactor)
z$tval <- z$coefficients / z$se
z$pval <- 2 * pt(abs(z$tval), z$df, lower.tail = FALSE)
z$terms <- mt
z$fe <- fl
z$N <- N
z$p <- p + numdum
z$xp <- p
z$call <- match.call()
class(z) <- "felm"
return(z)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.