Showing posts with label Chamberlain Munlack. Show all posts
Showing posts with label Chamberlain Munlack. Show all posts

Thursday, December 5, 2013

Incidental Parameters Problem with Binary Response Data and Unobserved Individual Effects

It is a well known problem that in some models as the number of observations becomes large, econometric estimators fail to converge on consistent estimators.  The leading case of this is when estimating a binary response model with panel data with potential "fixed effects" correlated with the explanatory variables appearing in the population.

One method that is typically implemented by researchers is to observe inviduals or organizations over multiple periods of time.  This is called panel data.  Within the context of panel data it is often assumed unobserved effects: genetics, motivation, business structure, and whatever other unobservables that might be correlated with the explanatory variables are unchanging over time. 

If we do then it is relatively easy to remove the effect of unobservables from our analysis.  In my previous post I demonstrate 3 distinct but equivalent methods for accomplishing this task when our structural model is linear.

However, when our model is a binary response variables (graduate from college or not, get married or not, take the job or not, ect.) it is usually no longer logically consistent to stick with a linear model.

In addition, not all of our remidies which worked for the linear model provide consistent estimators for non-linear models.  Let us see this in action.

First we will start with generating the data as we did in the December 4th post.

# Let's say: x = x.base + fe
 
set.seed(2)
 
nperson <- 500 # Number of persons
nobs <- 5      # Number of observations per person
 
 
fe.sd <- 1 # Spefify the standard deviation of the fixed effed
x.sd  <- 1 # Specify the base standard deviation of x
 
# First generate our data using the time constant effects
constantdata <- data.frame(id=1:nperson, fe=rnorm(nperson))
 
# We expand our data by nobs
fulldata <- constantdata[rep(1:nperson, each=nobs),]
 
# Add a time index, first define a group apply function
# that applies by group index.
gapply <- function(x, group, fun) {
  returner <- numeric(length(group))
  for (i in unique(group)) 
    returner[i==group] <- get("fun")(x[i==group])
  returner
}
 
# Using the generalized apply function coded above
fulldata$t <- gapply(rep(1,length(fulldata$id)), 
                     group=fulldata$id, 
                     fun=cumsum)
 
# Now we are ready to caculate the time variant xs
fulldata$x <- fulldata$fe + rnorm(nobs*nperson)
 
# This is were our simulation diverges from a linear model.
 
# Let us define the linear component as:
lc <- .5*fulldata$x + .5*fulldata$fe
 
# Now we will define the logit probability as
fulldata$pl <- exp(lc)/(1+exp(lc))
 
# Now we define the probit probability as
fulldata$pp <- pnorm(lc)
 
cor(fulldata$pl,fulldata$pp)
 
plot(fulldata$pl,fulldata$pp, main="Probit and Logit Probabilities", 
     xlab="Logit", ylab="Probit")
 
plot(sort(lc), sort(fulldata$pl), main="Probit and Logit Probabilities",
     ylab="P(y=1|x,a)", xlab="xb+a",
     type="l", lwd=2, lty=3)
  lines(sort(lc), sort(fulldata$pp), lwd=2, lty=2)
legend(-3.3,.8,c("logit","probit"), lty=c(3,2), lwd=2) 


# We can see the probit tends to have a shorter range in which the 
# action is happening. 
  
# We should really think of the probit and the logit as now
# two different sets of data.  Now let's generate out outcomes.
fulldata$yl <- fulldata$pl>runif(nobs*nperson)
fulldata$yp <- fulldata$pp>runif(nobs*nperson)
 
# Let's try to estimate our parameters with the logit.
glm(yl ~ x, data = fulldata, family = "binomial")
 
# We can see our estimator is upwards biased as we expect.
 
# Now we will try to inlcude fixed effects as if we were not
# aware of the incidental parameters problem.
glm(yl ~ x+factor(id), data = fulldata, family = "binomial")
 
# We can see, that including a matrix of dummy variables
# seems to actually make our estimator worse.
 
# Instead let's try the remaining fix that is available to us
# from the previous post listing 3 fixes.
 
# We will include an average level of the explanatory variable
# for each individual.  This is referred to as the 
# Chamberlain-Munlak device.
fulldata$xmean <- ave(fulldata$x, group=fulldata$id)
 
glm(yl ~ x+xmean, data = fulldata, family = "binomial")
# We can see that including an average effect significantly
# reduces the inconsistency in the estimator.
 
# Now, let's see what happens if we do the same things in the
# probit model.
glm(yp ~ x, data = fulldata, family = binomial(link = "probit"))
# Probit experiences similar upward bias to that of the logit.
 
glm(yp ~ x+factor(id), data = fulldata, 
    family = binomial(link = "probit"))
 
glm(yp ~ x+xmean, data = fulldata, family = binomial(link = "probit"))
# Interestingly, including the Chamberlain-Munlak device in the probit
# though theoretically inconsistent does seem produce estimates
# comparably good as including the device with the logit at least
# in the sample sizes simulated here.
Created by Pretty R at inside-R.org

Wednesday, December 4, 2013

Unobserved Effects With Panel Data

It is common for researchers to be concerned about unobserved effects being correlated with observed explanatory variables.

For instance, if we were curious about the effect of meditation on emotional stability we may be concerned that there might be some unobserved factor such as personal genetics that might  predict both likelihood to meditate and emotional stability.

In order to remove this potentially biasing effect we could think about taking measurements over multiple periods for the same individual inquiring about frequency of meditation and emotional stability.

If we observe that within the same individual, removing the time constant effects which (presumably) genetics is a component of that there is still a relationship between meditation and emotional stability, then we may feel on firmer ground as to our hypothesis that mediation may lead to more emotional stability.

In order to accomplish the goal of estimating this relationship we may experiment with a "fixed effects" model defined as:

$$y_{it}=x_{it}\beta + a_i+u_{it}$$

In this typical linear model with panel data, there is no problem including an arbitrary number of dummy variables.  Let's see this in action.

nperson <- 300 # Number of persons
nobs <- 3      # Number of observations per person
 
# In order for unobserved person effects to be a problem they must be
# correlated with the explanatory variable.
 
# Let's say: x = x.base + fe
 
fe.sd <- 1 # Spefify the standard deviation of the fixed effed
x.sd  <- 1 # Specify the base standard deviation of x
 
beta <- 2
 
# First generate our data using the time constant effects
constantdata <- data.frame(id=1:nperson, fe=rnorm(nperson))
 
# We expand our data by nobs
fulldata <- constantdata[rep(1:nperson, each=nobs),]
 
# Add a time index, first define a group apply function
# that applies by group index.
gapply <- function(x, group, fun) {
  returner <- numeric(length(group))
  for (i in unique(group)) 
    returner[i==group] <- get("fun")(x[i==group])
  returner
}
 
# Using the generalized apply function coded above
fulldata$t <- gapply(rep(1,length(fulldata$id)), 
                         group=fulldata$id, 
                         fun=cumsum)
 
# Or a more simplified function
indexer <- function(group) {
  returner <- numeric(length(group))
  for (i in unique(group)) 
    returner[i==group] <- 1:sum(i==group)
  returner
}
 
# Is a special case of gapply
fulldata$t <- indexer(fulldata$id)
 
# Now we are ready to caculate the time variant xs
fulldata$x <- fulldata$fe + rnorm(nobs*nperson)
 
# And our unobservable error
fulldata$u <- rnorm(nobs*nperson)
 
# Finally we are ready to simulate our y variables
fulldata$y <- .5*fulldata$x + .5*fulldata$fe + fulldata$u
 
# First lets see how our standard linear model performs:
summary(lm(y~x, data=fulldata))
 
# Adding a dummy variable removes the bias
summary(lm(y~x+factor(id), data=fulldata))
 
# The same result can be taken by removing the mean from
# both the explanatory variables and the dependent variables.
# Why is that?


Think of the problem as:
$$y_{it}=x_{it}\beta + a_i+u_{it}$$
So $$y_{it}-mean_t(y_i)=(x_{it}-mean_t(x_{i}))\beta + a_i-mean_t(a_i)+u_{it}-mean(u_i)$$

Because the unobservable effect is constant over time it drops out.  And so long as their was a term controlling for the average unobservable effect (the dummy variables) then the average per person unobserved error must by definition be equal to zero.

thus: $$y_{it}-mean_t(y_i)=(x_{it}-mean_t(x_{i}))\beta + u_{it}$$

fulldata$ydemean <- fulldata$y-ave(fulldata$y, group=fulldata$id)
fulldata$xdemean <- fulldata$x-ave(fulldata$x, group=fulldata$id)
 
summary(lm(ydemean~xdemean-1, data=fulldata))
 
# We can also accomplish this by adding the Chamberlain device to the 
# that regression is the total or mean of the explanatory variables at
# the level of each individual.
fulldata$xmean <- ave(fulldata$x, group=fulldata$id)
 
fulldata$xsum <- gapply(fulldata$x, group=fulldata$id, fun=sum)
 
# This is a little trickier to figure out how it accomplishes the task
# of differencing out the unobserved effect.  
 
# This is how I think of it. The unobserved individual effect must be 
# correlatedwith the explanatory variable in aggrogate to be a problem. 
# However, that correlation can only be on the individual level since
# by definition the "fixed effect" is constant on the individual level.
# Thus by creating a new variable which is the average or total for
# each individual, we are allocating to that variable any variation
# which correlates with the explanatory variable. 
 
summary(lm(y~x+xmean, data=fulldata))
summary(lm(y~x+xsum, data=fulldata))
Created by Pretty R at inside-R.org

Saturday, February 2, 2013

Chamberlain Mundlak Device and the Cluster Robust Hausman Test


* Unobserved variation can be divided in a useful manner.

* y_it = X_it*B + c_i + u_it

* c_i is fixed individual effect or random individual effect if c_i is uncorrelated with Xi (the time averages of X_it).

* In order for c_it to bias B is if there is some correlation between X_it and c_i.

* Therefore if we were to create a new variable which "controls" any time constant variation in X_it then the remaining v_i must be uncorrelated with X_it.

* Thus the Chamberlain Mundlak Device was born.

* Let's see it in action!

clear

set obs 100

gen id = _n

gen A1 = rnormal()
gen A2 = rnormal()

* Let's say we have 2 observations per individual
expand 2

gen x = rnormal()+A1
gen u = rnormal()*3

gen y = -5 +2*x + A1 + A2 + u

* Now in the above model there is both a portion of the unobserved variance correlated with the average x (A1) and a random portion uncorrelated with the average x (A2) the individual level.

* The fixed effect varies with the x variable while the random one does not.

* The standard approach in this case would be to use the Hausman test to differentiate between fixed effect and random effect models.

xtset id
* Let's first set id as the panel data identifier.

xtreg y x, fe
estimates store fe
* We store the estimates for use in the Hausman test

xtreg y x, re

hausman fe, sigmamore
* We strongly reject the null which we should expect so in classical econometric reasoning we choose to use the fixed effect estimator.

* An alternative method of estimating the fe estimator is by constructing the Chamberlain-Mundlak device.

* This device exploits the knowledge that the only portion of the time constant variation in X that can be correlated with u must be correlated only with the time average X for each individual.

bysort id: egen x_bar = mean(x)

reg y x x_bar

* Amazingly we can see that the new estimator is the same as the fe estimator above.

* Notice however the degrees of freedom.

* In the fe esimator we have used up half of our degrees of freedom.

* Yet our x estimate is the same size and our standard errors are very similar?

* If we double our sample size should not our standard errors decrease substantially?

* The answer is no. Why?

* I am going to run the fixed effect estimator manually.

reg y x i.id

* Look at our SSE or the R2.  In the fixed effect model the R2 is much larger.

* This is because in terms of the random effect (A2), the fixed effect model controls for both the portion of the unobserved individual level variance which is correlated with the average x for each student as well as that portion uncorrelated with the average x.

* The Chamberlain-Mundlak (CM) device however only controls for the portion of the variance correlated with the average Xs.  Thus there is much more unexplained variance which ends up reducing the power of our test which is approximately accounted for when we adjust the sample size.

* The CM can be additionally useful because it provides an alternative form of the Hausman test.

reg y x x_bar

* The significance of the generated regressor x_bar indicates the exogeneity of the unobserved individual effects.

* The test can be easily adjusted to be robust to cluster effects by specifying cluster in the regression.

reg y x x_bar, cluster(id)

Thursday, May 17, 2012

Unobserved fixed effects model


* Often times we are concerned that there are some unobserved
* factors which are correlated with our explanatory variables x
* as well as with our error term u.

* For example, we might be concerned that intelligence is
* correlated with years of schooling as well as future
* expected earnings.  However, fortunately, intelligence
* is thought of as a time constant factor.

* Therefore, if we remove time constant factors we might
* be able to approximate the returns to education.
* (This is assuming the returns to years of education is
* constant.  If it is a function of intelligence then
* we are going to need to think about being more clever
* about this.)

* Stata code
clear

* Imagine we have 200 individuals that we track
set obs 200
set seed 101

gen c = rnormal()
  label var c "Time Constant Heterogeniety (individual specific)"

gen id = _n
  label var id "Individual specific ID"

* create 5 observations for each initial observation
expand 2

bysort id: gen year=_n
  label var year "Year of observation"

tab year

gen x = rnormal()+c
  label var x "explanatory variable X (with time constant and time varying components)"

gen u = 3*rnormal()+3*c
  label var u "Error term (correlated with unobservables c)"

gen y = x + u
  label var y "Outcome variable"

reg y x
* We can see that OLS is biased

xtset id year
* Tells stata to use id as a panel data individual identifier

xtreg y x, fe
* However, the fixed effect estimator is unbiased because it
* successfully eliminates the correlation between the time
* constant correlation between the x and the error u.

* Note: an identical command is:
reg y x i.id

* Or:
areg y x, absorb(id)

* An alternative approach is the Chamberlain Munlack device.
* If we fear that the constant part of x might be correlated
* with u then we can easily control for that by including
* it in the regression:
bysort id: egen x_mean = mean(x)

reg y x x_mean

* When there is only two time periods difference in difference
* the same but in time periods more than two it tends to be
* different.  Though it is also effective at removing time
* constant effects.
gen y_diff = y-l.y
gen x_diff = x-l.x

reg y_diff x_diff