Thursday, 19 August 2010

New work in progress up: "Group Reciprocity"

We've got a first draft of the paper on Group Reciprocity up. I blogged about the results from the first experimental session before. The end results aren't quite so strong but are still interesting.

Experimental economists have worked a lot on reciprocity: when you harm (or help) me, I harm (help) you back. Group reciprocity adds a twist: you harm me, so I want to harm somebody else in your group. Group could be almost anything: gender, ethnicity, football team supporters ....

Group reciprocity might be an important cause of things like war, ethnic tension and racial discrimination, so it would be useful if we could study it in the lab. The paper is a first attempt to do that. I think there's still a lot more to be done.

Atlanta Race Riot 1906


Friday, 13 August 2010

Calculating the standard error of marginal effects using the delta method, in R

Here's a little piece of code I used to calculate marginal effects on probabilities in a non-linear (bivariate probit) model, using the Delta Method.

Warning: I am not a proper statistician. Use at your own risk!

# to create the model:
library(VGAM)
mbvp.1 <- vglm(cbind(decision_matrices.1-1, decision_matrices.6-1) ~ defect_other, binom2.rho, data=ds, trace=TRUE)

# needed to compute the gradient of my function numerically
library(numDeriv)
# needed for a multivariate normal distribution
library(mvtnorm)
# bv is a vector of betas
# this calculates the change in predicted probabilities
f <- function (bv) {
    rho <- rhobit(bv[3], inverse=T)
    p <- array(dim=c(2,2,2))
    for (c1 in 0:1) {
        for (c2 in 0:1) {
            for (d_o in 0:1) {
                sgn1 <- 2*c1 - 1
                sgn2 <- 2*c2 - 1
                ############################################
                # the probability of making choice (c1,c2) 
                # when the independent variable
                # defect_other=d_o:
                ############################################
                p[c1+1,c2+1,d_o+1] <- pmvnorm(upper=c(sgn1*(bv[1]+bv[4]*d_o),sgn2*(bv[2]+bv[5]*d_o)), corr=matrix(c(1,sgn1*sgn2*rho,sgn1*sgn2*rho,1), nrow=2))
            }
        }
    }
    # I use the global firstdec so that the function conforms to what grad() expects
    # the total probability of choosing 1 for either the first
    # or the second decision:
    if (firstdec) p[2,2,2]+p[2,1,2]-p[2,2,1]-p[2,1,1] else p[2,2,2]+p[1,2,2]-p[2,2,1]-p[1,2,1]
}


firstdec <- TRUE
pred[1] <- f(coef(mbvp.1))
gradient <- grad(f, coef(mbvp.1))
# ===== this is the delta method ====
predse[1] <- sqrt(gradient %*% vcov(mbvp.1) %*% gradient)

firstdec <- FALSE
pred[2] <- f(coef(mbvp.1))
gradient <- grad(f, coef(mbvp.1))
predse[2] <- sqrt(gradient %*% vcov(mbvp.1) %*% gradient)

margprobs <- cbind(pred, se=predse, Z=pred/predse, p=2*pnorm(-abs(pred/predse)))
print(margprobs, digits=3)