Using refset in your R analyses

refset is an R package that creates subsets which refer to an original dataset. Here is a short example of how to use refset to simplify your data analysis.

I have a dataset from a laboratory experiment. In total, 532 subjects participated. Each subject was in a group of 4 and made decisions in 20 rounds. There were several treatments.

I typically create two scripts: R-create-data.R makes a data frame from my raw files, including any computed variables, and saves it to a .RData file. R-analysis.R loads that file and runs any statistics or graphs. That way, I can rerun my analyses without recreating the data.

At the end of R-create-data.R, I have made a data frame called subj. A row of subj records one decision by one subject, so it contains 532 * 20 = 10640 rows. The script then create subsets of subj to represent, e.g. the data from particular treatments:

subjhist <- subj[subj$treatment1=="H",]
subjPX <- subj[subj$treatment1=="PX",]
subjn <- subj[subj$treatment1=="N",]
subjst <- subj[subj$treatment1=="S",]
save(subj, subjhist, subjPX, subjn, subjst,    file="data.RData")

In theory, the script has created all the variables I need, and now I can just run regressions using whichever subset is appropriate. In practice, I will look at the data and think of a new analysis I should do, which requires a new variable.* Hmm... did my subjects play differently in earlier than later rounds?

When that happens, I need to create a new variable in all my datasets:

subj$lateround <- subj$Period > 10
subjPX$lateround <- subjPX$Period > 10 subjhist$lateround <- subjhist$Period > 10
## et cetera...

Eventually, I will get round to including lateround in R-create-data.R, rerunning that file, and removing the code above from my analysis script. But in the meantime, doing this is a hassle, and creates repetitive, verbose, unreadable code.

Let's replace those subsets by refsets.

library(refset)
subjhist %r% subj[subj$treatment1=="H",]
subjPX %r% subj[subj$treatment1=="PX",]
subjn %r% subj[subj$treatment1=="N",]
subjst %r% subj[subj$treatment1=="S",]

Each of these variables now refers to the main data frame subj, rather than being a copy of it.

We'll put this in R-analysis.R, so that it happens after our data is loaded. Otherwise, saving and loading breaks the connection between refset and data frame, and simply creates subsets as before. The data will look just the same, and I can run my existing analyses unchanged.

Now, when I create a new variable, I just do it in the main data frame:

subj$lateround <- subj$Period > 10

The variable is then visible in all of my refsets:

head(subjPX$lateround, 20)
## [1] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE  TRUE  TRUE  TRUE
## [14]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
## same for subjhist, subjn...

Problem solved: I can create new variables once, not several times, and immediately use them in any analysis.

* Yes, I know... one should have prepared all my analyses ex ante and simply run them at the press of a button. Otherwise one is just doing post hoc exploration and curve fitting. I applaud the ideal and strive to live up to it, but the real world isn't always like that -- and in any case, sometimes post hoc exploration is worth doing.

Honore-style fixed effects estimators for censored data in R

Below is some R code to implement fixed effects Tobit, a la Honore (1992) "Trimmed LAD and Least Squares Estimation of Truncated and Censored Regression Models with Fixed Effects", and also fixed effects double-censored Tobit a la Alan, Honore and Petersen (2008).

# FE tobit a la Honore (1992). See Honore (2002) Port Econ J for a simple guide

# y1 partner y2 other

honore <- function (b, x1, x2, y1, y2) {

dxb <- (x1 - x2) %*% b

sum(

(pmax(y1, dxb) - pmax(y2, - dxb) - dxb)^2 +

2*(y1 < dxb)*(dxb-y1)*y2 +

2*(y1 < -dxb)* (-dxb-y2)*y1

)

}

fetobit <- function (x1, x2, y1, y2) {

# x <- model.matrix(form, dataset)

initvals <- lm.fit(x1 - x2, y1 - y2)$coefficients

res <- optim(initvals, fn=honore, x1=x1, x2=x2, y1=y1, y2=y2,

method="BFGS", control=list(reltol=1e-8, maxit=1000))

if (res$convergence != 0) warning("Didn't converge")

res$par

}

# estimate the standard error, asymptotically.

fetstderr <- function(x1, x2, y1, y2, bhat) {

N <- nrow(x1)

x <- x1 - x2

xb <- x %*% bhat

G <- matrix(0, nrow=length(bhat), ncol=length(bhat))

V <- G

# this should be done faster...

for (i in 1:N) {

xtx <- outer(x[i,], x[i,])

G <- G + ( - y2[i] < xb[i] & xb[i] < y1[i] ) * xtx

V <- V + ((y2[i]^2) * (y1[i] <= xb[i]) +

(y1[i]^2) * (xb[i] <= - y2[i]) +

(y1[i] - y2[i] - xb[i])^2 * (- y2[i] < xb[i] & xb[i] < y1[i])

) * xtx

}

G <- G / N

V <- V / N

Ginv <- solve(G)

sqrt(diag(Ginv %*% V %*% Ginv / N)) # right?

}

# two-sided censoring

# y's must be censored at 0 and 1

alan <- function (b, x1, x2, y1, y2) {

d <- pmax(-1, pmin(as.matrix(x1-x2) %*% as.matrix(b), 1))

# sum of R(y1, y2, max(-1, min(,(x1-x2) %*% b, 1)))

r1 <- r2 <- rep(NA, length(d))

r1 <- ifelse(d <= y1 -1, d + d^2/2 +(y1-1)^2/2,

ifelse(y1-1 < d & d <= 0, d*y1,

ifelse(0 < d & d <= y1, d*y1 - d^2/2,

ifelse(y1 < d, 2,<="" div="" y1^2="">

NA))))

r2 <- ifelse(d <= -y2, -y2^2/2,

ifelse(-y2 < d & d <= 0, d^2/2 + d*y2,

ifelse(0 < d & d <= 1-y2, d*y2,

ifelse(1-y2 < d, (y2-1)^2="" -="" 2,<="" 2="" d="" d^2="" div="">

NA))))

- sum(r1 - r2) # Alan et al say minimize, but I say maximize, ie minimize the minus sum

}

fe2tobit <- function (x1, x2, y1, y2) {

initvals <- lm.fit(x1 - x2, y1 - y2)$coefficients

res <- optim(initvals, fn=alan, x1=x1, x2=x2, y1=y1, y2=y2,

method="BFGS", control=list(reltol=1e-8, maxit=1000))

if (res$convergence != 0) warning("Didn't converge")

res$par

}

fet2stderr <- function (x1, x2, y1, y2, bhat) {

N <- nrow(x1)

x <- x1 - x2

xb <- x %*% bhat

V <- G <- matrix(0, nrow=length(bhat), ncol=length(bhat))

u1 <- ifelse(xb > 0, pmax(y1-xb,0), pmin(y1,1+xb))

u2 <- ifelse(xb > 0, pmin(y2, 1-xb), pmax(y2+xb, 0))

u <- u1 - u2

gwt <- ifelse(xb < -1 | xb > 1, 0,

(-1 < xb & xb < y1-1) - (0 < xb & xb < y1) -

(-y1 < xb & xb < 0) + (1-y2 < xb & xb < 1))

for (i in 1:N) {

xtx <- outer(x[i,], x[i,])

G <- G + gwt[i] * xtx

if (-1 < xb[i] & xb[i] < 1) {

v <- u[i]*x[i,]/2

V <- V + outer(v, v)

}

}

G <- G / N

V <- V / N

Ginv <- solve(G)

sqrt(diag(Ginv %*% V %*% Ginv / N))

}

Points to note:

• fetobit expects data to be censored at 0; fe2tobit expects censoring at 0 and 1. There should be no NAs. 
• Warning! I have done only a very quick test on these. They seem to work but I am not a proper econometrician. Use at your own risk. 
• The standard errors may be inaccurate for low N. Consider bootstrap estimation. 
• Warning! Google searches for "Trimmed LAD" may lead to unexpected results.