Showing posts with label fixed effects. Show all posts
Showing posts with label fixed effects. Show all posts

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

Monday, January 28, 2013

HLM comparison with OLS - 2 levels, random coefficient on constant

do file

* xtmixed is capable of estimating a variance for the random effects on multiple levels.

* random effects are normalized to have mean 0.

* Our initial model is y_ij = bij + xij*B1 + uij

* With bij = B0 + v0j

* We can represent our model as:

* y_ij = B0 + v0j + xij*B1 + uij = B0 + xij*B1 + v0j + uij

* With v0j + uij an unobserved error term.

* Let's see an example.

clear

set obs 60
* We will have 60 different schools (level 2 indexes)

gen school=_n

gen v0j = rnormal()*2
* Generate some school effect with standard deviation = 2

expand 200
* There are 200 individuals in each school.

gen x1 = rnormal()
* Each individual has a continuous predictor (say academic performance).

gen u = rnormal()*5
* each individual has an unobserved error.

gen y = 1 + x1*2 + v0j + u
* Each individual's outcomes are a function of individual ability plus the school effect plus individual error.

xtmixed y x1 || school:
* Now let's estimate both the effect of x1 on y as well as the effect of school variation on predicting outcomes.

* The standard deviation of the estimate of school effect is close to 2 which is the true so xtmixed is working well.

* In principal we can attempt to estimate the same thing using OLS with dummies.

qui tab school, gen(sch_id)
* Generate a set of dummy indicator variables for each school.

* School 1 is left out of the estimation as the reference school.
reg y x sch_id2-sch_id60

predict u_hat, resid

* We would like to calculate the standard deviations of the school estimates in order to compare with the xtmixed standard deviation estimates.

* Since we left school 1 out of the estimation it is considered a 0 effect.  All other school effect estimates are relative to school 0.
gen sch_coef_est = 0 if school==1

* The following code will save the school estimates to the coefficient variable that can then be used to find the standard deviation of the estimated school effects.
forv i=2/60 {
  qui replace sch_coef_est = _b[sch_id`i'] if school==`i'
}

sum sch_coef_est u

* We can see that the standard deviation of sch_coef_est is similar to that estimated by xtmixed as well as the standard deviation of the residuals.

* So both methods seem to be effective at estimating the variance of the school level effect.

corr sch_coef_est v0j

* We can also see that the OLS dummy variable method has produced individual school effects that highly correlated with the true school effects.

* In summary both methods seem to work well though we would probably favor the xtmixed (hierarchical linear model) estimator if we did not care about the actual individual school estimates because it provides an estimator that maintains more degrees of freedom.

Wednesday, October 24, 2012

Hierarchical linear modelling


* Explorations of in heirachical linear modeling

* Hierarchical modeling is a a type of model that assigns to different levels different portions of the unexplained error term.

* What does this mean?

* Imagine that you are interested in predicting student success.

* However, you believe that there is a invidual effect size for each student, each teacher, each school, and each district.

* You would like to know what portion of the unexplained variance can be attributed to each level.

* Let's see how this works

clear

* Imagine that there are 5 districts you have data on
set obs 5

gen dist_id=_n

* Each district has an effect size sd=2
gen dist_fe = rnormal()*2

* Each district has 5 schools
expand 5

gen school_id=_n

* Each school has an effect size sd=3
gen school_fe = rnormal()*3

* Within each school you have 8 teachers
expand 8

gen teach_id = _n

* Each teacher has an individual effect size sd = 4
gen teacher_fe = rnormal()*4

* Each teacher has twenty students.

expand 20

gen student_id = _n

* Each of these students has an effect size equal to sd = 5

gen student_fe = rnormal()*5

* Let's imagine for a second that each student stays with the same teacher for two semesters and we have information about student progress at the end of each semester.

expand 2

sort student_id

gen semester = mod(_n-1,2) + 1

* There is some transient error (shocks) that affect student performance during each semester sd = 6

gen u = rnormal()*6

* Finally let's imagine that there is some treatment such as extra tutoring that is randomly assigned to 20% of our students on a semester basis.

gen tutoring = rbinomial(1,.2)

* Let's now generate our test results

gen perform = 2*tutoring + dist_fe + school_fe + teacher_fe + student_fe + u

* For our first cut at our analysis of this data let's compare the mixed effect HLM model with a standard fixed effect model and random effects model.

* The standard fixed effect model estimates an effect for every individual.  In this case let us focus on students.

xtset student_id

xtreg perform tutoring, fe
* We can see that our estimate of the effect tutoring is close and that our variance in estimated student effects (sigma_e) is very close to 6.

* We might also try to estimate tutoring with a random effects model at the student level.
xtreg perform tutoring, re

* We expect that the re estimator outperforms the fe estimator in estimating the returns to tutoring because in this case we know that the assignment to tutoring is trully orthoganol to the student fixed effects.

* Now let's compare our previous estimates to the single level mixed effects
xtmixed perform tutoring || student_id:

* We can see that the xtmixed command produces an estimate on tutoring very nearly identical to that of the random effects model.

* As with the previous two estimates, we are generally satisfied with the estimates of the standard deviation of the student effect.

* Now let's start adding layers.

xtmixed perform tutoring || student_id: || teach_id: || school_id: || dist_id:

* Our estimate of tutoring has not improved.

* Likewise our ability to estimate the variance at each level does not seem to have much ability to distinguish effect size between the district level with a (sd of 2), the school (sd of 3), the teacher (sd of 4), and the student (sd of 5).

* Perhaps this is a result of the sample size being too small.

* At this point I will return to the beginning of the simulation and increase the number of districts to 50 with the resulting number of observations at 80,000.  However, the xtmixed command is very slow and computationally intensive so I will paste in the results directly.

/* The First etimation
Performing EM optimization:

Performing gradient-based optimization:

Iteration 0:   log restricted-likelihood =  -283974.7
Iteration 1:   log restricted-likelihood =  -283974.7

Computing standard errors:

Mixed-effects REML regression                   Number of obs      =     80000
Group variable: student_id                      Number of groups   =     40000

                                                Obs per group: min =         2
                                                               avg =       2.0
                                                               max =         2


                                                Wald chi2(1)       =    873.44
Log restricted-likelihood =  -283974.7          Prob > chi2        =    0.0000

------------------------------------------------------------------------------
     perform |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
    tutoring |   1.975427   .0668413    29.55   0.000     1.844421    2.106434
       _cons |  -.4480866   .0439606   -10.19   0.000    -.5342478   -.3619255
------------------------------------------------------------------------------

------------------------------------------------------------------------------
  Random-effects Parameters  |   Estimate   Std. Err.     [95% Conf. Interval]
-----------------------------+------------------------------------------------
student_id: Identity         |
                   sd(_cons) |   7.231256   .0354216      7.162163    7.301015
-----------------------------+------------------------------------------------
                sd(Residual) |   5.984532   .0211588      5.943205    6.026147
------------------------------------------------------------------------------
LR test vs. linear regression: chibar2(01) = 17369.05 Prob >= chibar2 = 0.0000

* The second estimation (this took a very long time)

Performing EM optimization:

Performing gradient-based optimization:

Iteration 0:   log restricted-likelihood = -284364.73
Iteration 1:   log restricted-likelihood = -283987.22  (not concave)
Iteration 2:   log restricted-likelihood = -283981.69  (backed up)
Iteration 3:   log restricted-likelihood =  -283974.7  (not concave)
Iteration 4:   log restricted-likelihood =  -283974.7  (backed up)
Iteration 5:   log restricted-likelihood =  -283974.7  (not concave)
Iteration 6:   log restricted-likelihood =  -283974.7  (backed up)

Computing standard errors:

Mixed-effects REML regression                   Number of obs      =     80000

-----------------------------------------------------------
                |   No. of       Observations per Group
 Group Variable |   Groups    Minimum    Average    Maximum
----------------+------------------------------------------
     student_id |    40000          2        2.0          2
       teach_id |    40000          2        2.0          2
      school_id |    40000          2        2.0          2
        dist_id |    40000          2        2.0          2
-----------------------------------------------------------

                                                Wald chi2(1)       =    873.44
Log restricted-likelihood =  -283974.7          Prob > chi2        =    0.0000

------------------------------------------------------------------------------
     perform |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
    tutoring |   1.975427   .0668413    29.55   0.000     1.844421    2.106434
       _cons |  -.4480866   .0439606   -10.19   0.000    -.5342478   -.3619253
------------------------------------------------------------------------------

------------------------------------------------------------------------------
  Random-effects Parameters  |   Estimate   Std. Err.     [95% Conf. Interval]
-----------------------------+------------------------------------------------
student_id: Identity         |
                   sd(_cons) |   3.621455   11.10185      .0089014     1473.35
-----------------------------+------------------------------------------------
teach_id: Identity           |
                   sd(_cons) |   3.608619   9.417657      .0216722    600.8684
-----------------------------+------------------------------------------------
school_id: Identity          |
                   sd(_cons) |   3.681086   8.431311       .041338    327.7954
-----------------------------+------------------------------------------------
dist_id: Identity            |
                   sd(_cons) |   3.550187   8.323714       .035854    351.5322
-----------------------------+------------------------------------------------
                sd(Residual) |   5.984525   .0211586      5.943198    6.026139
------------------------------------------------------------------------------
LR test vs. linear regression:       chi2(4) = 17369.05   Prob > chi2 = 0.0000

Note: LR test is conservative and provided only for reference.              */

* So?  As a result we can see that our estimators do not improve substantially as a result of more data.  They may even be worse.

Monday, September 3, 2012

Robust Hausman Test


* The Huasman test is a commonly used to indicate an ideal choice between fixed effect and random effect estiamtors (in a panel data context).  This robust estimator was first proposed by Arellano (1993) {http://ideas.repec.org/a/eee/econom/v59y1993i1-2p87-97.html}.

* If I understand this properly, the RE estimator is a GLS estimator that should only be used when the individualized effect of each person (referred to as their fixed effect) is uncorrelated with the explanatory variables and uncorrelated with the outcome variables.

* This exogeneity of individual heterogeneity is often better understood in the situations when it fails rather than when the assumption is upheld.

* Imagine that motivation is relatively constant for individuals.

* If we have multiple years of GPA, which we are trying to predict and number of hours spent studying, then accross individuals it might be difficult to estimate GPA as a function of hours worked if we ignore the unobserved factor motivation because motivation may cause individuals to both study more hours and do better in general regardless of hours spent studying.

* Let's see a simple simulation of this:

clear
set obs 10000

gen id=_n

gen motivation = runiform()
  label var motivation "Unobserved student motivation"

expand 3
* We have three years of data per student

* The more motivated students are the more they study
gen hours_study = runiform()*2+motivation

gen attendance = runiform()

gen u = rnormal()*5*hours_study
* This is creating some heterogeneity in the error proportional to hours of study.

gen GPA = motivation + hours_study + attendance/2 + u

* Now we have this data that we are concerned might not be suitable for RE but we would like to if we could since RE is more efficient that FE when the assumptions are met.

xtset id

* Stata has a built in command to do the traditional Hausman test:
xtreg GPA hours_study attendance, fe
est store fe
xtreg GPA hours_study attendance, re
est store re
hausman fe re

* Alternatively using the Chamberlain-Munlack Device, we can do a similar estimation:
foreach v in hours_study attendance {
  bysort id: egen mean_`v' = mean(`v')
}

xtreg GPA hours_study mean_hours_study attendance mean_attendance, re
test mean_hours_study mean_attendance
* This test result is not exactly the same. I think it is due to the tests being asympotically equivalent while in finite samples, not equivalent.

* I think this second form of the test is more informative.  We are adding the mean values of each of our explanatory variables (by individual) and seeing if those mean values have additional explanatory power outside of that of their levels.

* This was somewhat disarming for me.  I thought, well what about the unexplained variation uncorrelated with the mean explanatory variables?

* Well, since a FE model can only control for fixed unexplained variation then controlling for that unexplained variation through use of means is surprisingly comprehensive.

* If the means of explanatory variables by individuals is uncorrelated with the error then using a fixed effect approach is not going to improve the estimation outcomes.

* The additional benefit of this form of the Hausman test is that it is extremely easy to make this estimator robust.

xtreg GPA hours_study mean_hours_study attendance mean_attendance, re vce(cluster id)
test mean_hours_study mean_attendance

* Since the mean variables are jointly significant, this suggests to us that we must assume there is unobserved heterogeneity that is correlated with the explanatory variable and the outcome variable and is therefore problematic to effective RE estimation, therefore FE is preferred.

* Note also, this same kind of logic can be applied to a decisions between FE and pooled OLS since it can be shown that RE is a weighted estimator between FE and Pooled OLS.