#' Interface between the Zelig Model twosls and
#' the Pre-existing Model-fitting Method
#' @param formula a formula
#' @param ... additonal parameters
#' @param data a data.frame
#' @return a list specifying '.function'
#' @export
zelig2twosls <- function (formula, ..., data) {
loadDependencies("systemfit")
# Helper function to perform set-difference
"%w/o%" <- function(x, y)
x[!x %in% y]
formula<-parse.formula(formula, "twosls")
tt<-terms(formula)
ins<-names(tt) %w/o% names(attr(tt,"depVars"))
if(length(ins)!=0)
if(length(ins)==1)
inst <- formula[[ins]]
else
inst <- formula[ins]
else
stop("twosls model requires instrument!!\n")
class(formula) <- c("multiple", "list")
# Return
list(
.function = "callsystemfit",
formula = formula[names(attr(tt,"depVars"))],
method = "2SLS",
inst = inst,
data = data,
...
)
}
#' @S3method param twosls
param.twosls <- function(obj, num=1000, ...) {
# Produce a vector of all terms
big.coef <- coef(obj)
# Produce a pretty sparse matrix containing 3 vcov matrices.
#
# Note that this matrix will give a value of zero to any invalid row-column
# combination.
# In particular, any terms that do not belong to the same equation will have
# a zero value.
big.vcov <- vcov(obj)
# This is a complete list of the terms. This is largely ignored, aside from
# the fact that we need a list of the formulae. In general, terms.multiple
# produced a pretty unwieldy list of items.
all.terms <- terms(obj)
# This list stores the results
simulations.list <- list()
# Iterate through the set of terms, and simulate each list separately.
for (key in names(all.terms)) {
# Extract the terms for an individual model.
eq.terms <- terms(all.terms[[key]])
# Extract the labels for the terms
eq.term.labels <- attr(eq.terms, "term.labels")
# Add the labeled for the intercept column, if it should exist
if (attr(eq.terms, "intercept"))
eq.term.labels <- c("(Intercept)", eq.term.labels)
# Format the title, this should look like:
# <list-item-name>_<term-label>
#
# So for the list: list(mu1 = y ~ x + sin(x))
# We get:
# "mu1_(Intercept)" "mu1_x" "mu1_sin(x)"
entries <- paste(key, eq.term.labels, sep = "_")
# Extract the mean-value of this term (from the lumped-toegether vector)
eq.coef <- big.coef[entries]
# Extract the vcov matrix of this term (from the lumped-together matrix)
eq.vcov <- big.vcov[entries, entries]
# Simulate the parameters
eq.simulations <- mvrnorm(num, eq.coef, eq.vcov)
# Name the columns
colnames(eq.simulations) <- eq.term.labels
# Add to the list
simulations.list[[key]] <- eq.simulations
}
# Return.
list(
coef = simulations.list,
linkinv = NULL
)
}
#' @S3method qi twosls
qi.twosls <- function(obj, x=NULL, x1=NULL, y=NULL, num=1000, param=NULL) {
# Compute the expected value of multistage LS methods
compute.ev <- function (obj, x, param) {
#
if (is.null(x) || is.na(x)) {
return(NA)
}
# If 'x' has too many rows, there will currently be errors. This is an issue
# in Zelig-core
if (nrow(x$matrix) > 1) {
warning("This package does not currently support pooled results.")
x <- x[1, ]
}
# Begin regular function
terms <- terms(obj)
# :q
coef.list <- coef(param)
# Hold Results
eta <- list()
#
for (key in names(coef.list)) {
#
coef <- coef.list[[key]]
# print(colnames(coef))
small.x <- as.matrix(x$matrix[, colnames(coef)])
#
eta[[key]] <- coef %*% (small.x)
}
# Convert list into a matrix
eta <- Reduce(function (x, y) cbind(x, y), eta)
colnames(eta) <- names(terms)
eta
}
ev1 <- compute.ev(obj, x, param)
ev2 <- compute.ev(obj, x1, param)
fd <- ev2 - ev1
# Name each column after the associated equation
# Return the results
list(
"Expected Values: E(Y|X)" = ev1,
"Expected Values: E(Y|X1)" = ev2,
"First Differences: E(Y|X1)-E(Y|X)" = ev2 - ev1
)
}
#' @S3method describe twosls
describe.twosls <- function (...) {
category <- "continuous"
description <- "Two Stage Least Squares"
authors <- c("Ferdinand Alimadhi", "Ying Lu", "Elena Villalon")
year <- 2007
package <-list(
name = "systemfit",
version = "0.8"
)
parameters <- list()
parameters$mu <-list(
equations=c(2,Inf),
tagsAllowed=TRUE,
depVar=TRUE,
expVar=TRUE
)
parameters$inst<-list(
equations=c(1,1),
tagsAllowed=FALSE,
depVar=FALSE,
expVar=TRUE
)
list(category = category, authors = authors, year = year, description = description, package = package, parameters = parameters)
}
#' @S3method plot sim.twosls
plot.sim.twosls <- function (x, ...) {
# Define locak function to plot a set of quantities of interest
plotSet <- function (title) {
for (col in colnames(qis[[title]])) {
q <- qis[[title]][, col]
plot(density(q), main = paste(col, title, sep=": "))
}
}
# Code begins here
cols <- c( rep(), rep(), rep() )
qis <- as.list.qi(x$qi)
qis <- Filter(function (y) any(!is.na(y)), qis)
qis <- Filter(is.matrix, qis)
max.cols <- max(unlist(Map(ncol, qis)))
layout.matrix <- matrix(0, length(qis), max.cols)
rownames(layout.matrix) <- names(qis)
count <- 1
for (title in names(qis)) {
for (k in 1:ncol(qis[[title]])) {
layout.matrix[title, k] <- count
count <- count + 1
}
}
layout(layout.matrix)
for (key in names(qis)) {
plotSet(key)
}
}
callsystemfit<-function(formula,data,method,inst=NULL,...){
# Call systemfit..
out <- systemfit(
data = data,
formula = formula,
method = method,
inst = inst,
...
)
# Assign class to formula, so that it is correctly parsed
class(formula) <- c("multiple", "list")
# Set the terms explicitly
attr(out,"terms") <- terms(formula)
# Set the class explicitly
class(out) <- c("multiple", class(out))
# Fin. Return the modified object
return(out)
}
as.list.qi <- function (x, names = "") {
class(x) <- "list"
indices <- attr(x, ".index")
attr(x, ".index") <- NULL
rename.keys(x, indices, names(indices))
}
rename.keys <- function (x, keys, to, warn = TRUE) {
all.names <- names(x)
indices <- match(keys, all.names)
if (any(is.na(indices)))
stop("Keys contains values that are not in `x`")
all.names[indices] <- to
names(x) <- all.names
x
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.