Open Data - Stack Exchange

As a fan of Stack Exchange I am very excited about the new exchange being made available. This exchange is called 'Open Data' and allows users to post questions on where to find publicly available data on a topic and other users will respond with one or more answers.

Open Data - Stack Exchange

Tobit fitted values not "fitting" data

* I was recently asked by a reader why it might be that the predicted values from 
* a tobit regression might have a constant significantly below zero and many fitted
* values unrealistically below zero.

* Is that a problem?

* Let's do a simple simulation of what the Tobit might be doing:

clear
set obs 10000

gen x = rnormal()

gen y_true = -2 + x*2 + rnormal()*4

gen y_observed = y_true
replace y_observed = 0 if y_true < 0

hist y_ob

tobit y_ob x, ll(0)

predict y_hat

graph twoway  (scatter y_hat y_true) (lfitci y_hat y_true)

* Not the best fit but, okay.

* So the take away. It is not inconsistent with the Tobit model at all that
* the constant and many fitted values may be significantly below zero.

* In a way, that is typically the result when it is most important to use a Tobit
* because you observe few values which are positive indicating that the underlying
* function is typically having many values fit below zero which have been censored.

* When you have the case in which the constant and many fitted values above zero,
* you have the least gains in terms of reducing bias from using a Tobit estimation
* method.

Formatted By Econometrics by Simulation

dapply - tapply for dataframes

# Every so often I run into a situation when I find I need a function that seems to me that it should be very straightforward but because it seems so much so, it also becomes extremely difficult to find online.
 
# In this case, I am interested in a function that applies a function across and index of values and returns the results of that function as an equal length as the vector fed into the function.
 
# In the search for this specific function, I found the base function tapply(vector, index, function) which does very close to what I want.
 
# However, tapply returns a vector of length equal to the number of factors in the index value.
 
# From the tapply function there might be a very easy way to get what I want.
 
# The reason I want this function is for data analysis purposes.  Often times we might want to see how well a variable that changes on the individual level (such as student education) might affect a variable that changes on the district level (per pupil spending).
 
# However, nothing occured to me.  So I programmed up an alternative function which has a few more bells and whistles.  (Check out the comments below for efficient coding methods).

 
# So hear it is.
 
dapply <- function(v1, index, fun, data=NULL, each=F) {
  # If a data.frame or list is defined then grab v1 and index from them.
  if (is.null(data)==F) {
    # The deparse and substitute commands together
    print(deparse(substitute(v1)))
    v1 <- data[[deparse(substitute(v1))]]
    index <- data[[deparse(substitute(index))]]
  }
  # Allow index to be a pattern in which case it is replicated accross the range of v1.
  if (length(index)<length(v1)) {
    if (length(index)*2>length(v1)) print(paste("Warning: Input vector is less than the index in length but not greater than twice."))
    # Calculate the number of replications required to duplicate the length of v1 with the index.
    ndup <- length(v1)/length(index)
    # If the number is not a whole number warn user.
    if (!is.integer(ndup)) print(gettextf("Input vector is of length %i can not be evenly divided by index of length %i.  Pattern truncated at end.", length(v1), length(index)))
    # Repeat the index as neccessary
    if (!each) index <- rep(index, times=ceiling(ndup))[1:length(v1)]
    if (each)  index <- rep(index, each =ceiling(ndup))[1:length(v1)]
  }
 
  # Calculate the vector which will be returned.
  vout <- tapply(v1, index, fun)
 
  # Expand the vout to cover multiple returns 
  returner <- NULL
  for (i in 1:length(index)) returner <- c(returner, vout[index[i]==names(vout)])
  returner
}
 
# For the vector 1:100 we will find the sum of the each block of twenty numbers.
dapply(1:100, ceiling(1:100/20),sum)
 
# We could have instead specified the index from 1:5 and by specifying each=T it tells it to multiply each element by what is neccessary to equal the length of the input vector.
dapply(1:100, 1:5, sum, each=T)
 
# Without each we will instead have the sum of all of the every fifth number.
dapply(1:100, 1:5, sum)
 
# But more in line with the uses I am looking for let's first construct a sample data frame:
district <- ceiling(1:1000/100)
stud.achievement <- rnorm(1000)+district*.1
sampleframe <- data.frame(id=1:1000, stud.achievement, district)
 
# Now let's see if we can find average student achievment by district and save it as a new variable.
sampleframe$Ave.ach <- dapply(stud.achievement, district, mean, data=sampleframe)
 
# Note, this function works very similarly to the egen command in Stata.

Syntax Highlighted by Pretty R at inside-R.org

# Gethub repository.

https://github.com/EconometricsBySimulation/2013-05-28-dapply

Quandl Package - 5,000,000 free datasets at the tip of your fingers!

# Yes, you read that correctly and no Quandl (http://www.quandl.com/) did not pay me anything.

# Quandl is a new database management tool which seeks to become the place to find datasets.  They boast of having over 5x10^6 data sets available though after examining them, I have decided that they are not entirely what everybody might think of as data sets.  That is, each unique indicator is considered an independent data set.  This helps them to seem to have a ginormous quantity of data sets.

# That said, they are not wrong in calling each indicator its own data set since much of their data, like financial data or government data is collected by disjoint teams.  The scope of their ambition is fantastic yet it is doable and frankly someone needed to do it.

# Currently, data seekers can access the Inter-University Consortium for Political and Social Research (IPCSR).  This great resource is composed mostly of cross section and panel data sets which are great for much analysis but IPCSR resricts access to data to member universities.  In addition, the kind of data that Quandl is indexing is a lot of data that would not show up on IPCSR database.  In addition, Quandl is integrating an automated structure that will be self-updating.

# For an example of how Quandl is a good step ahead of the game take a look at this search quiery:

http://www.quandl.com/search/lansing,%20michigan

# In this search, I searched out Lansing, Michigan where I live and returned results of data for the last decade or earlier up to today from sources such as the Federal Reserve and the US Energy Information Administration.

http://www.icpsr.umich.edu/icpsrweb/ICPSR/studies?q=Lansing%2C+Michigan&permit%5B0%5D=AVAILABLE

# In constrast when queirying ICPSR, I found a few databases listed but they were historical databases that spanned back generally between 30 and 70 years.  That said both sources could provide valuable information depending upon what I am interested in modeling.

# Quandl is very clever for a number of reasons.  One of these reasons is that they have simultaneously released 8 software packages that can be used in a number of statistical packages such as R, Stata, and Excel.

# In order to demonstrate the use of Quandl I will grab a few data sets from the Lansing quiery drawn from the Federal Reserve.

install.packages("Quandl")
library(Quandl)

# Employment numbers (thousands of people") for Lansing, Michigan
NonFarm = Quandl("FRED/LANS626NAN")
CivLaborForce = Quandl("FRED/LANS626LFN")
PerCapitaIncome = Quandl("FRED/LANS626PCPI")

# Now let's combine the data so that we can related data values.
Labor = merge(NonFarm, CivLaborForce, by="Date")
Combined = merge(Labor, PerCapitaIncome, by="Date")
colnames(Combined) = c("Date", "NonFarm", "CivLaborForce", "PerCapitaIncome")
  # Notice that though our data had many more data points, the default option of merge only keeps data that exists in both data sets.  In this case, it is per capital income that has the least number of data points.

# Let's see if we can predict income as a function of employment:
summary(lm(PerCapitaIncome~NonFarm+CivLaborForce, data=Combined))

# Our naive prediction as a result of this is that as the Civilian Labor Force increases, wages rise.  This is of course a naive example ignoring completely issues of causation and endogeneity not to mention probable random walks and other challenging features of this kind of data.

# The overall take away though, should be "cool", I think.  Maybe this data bank does not provide information currently on many issues of interest to those looking for data.  But it does make things easier and self-updating, which are great features.

Use Expand to Manually Reshape Data

do file

* Imagine that you have a wide data set that you would like to convert to a long data set.

* Your data set is structured in the following manner.

clear
set obs 1000

gen id = _n

forv i=1/4 {
  gen var1_`i' = rnormal()
  gen var2_`i' = rnormal()
  gen var3_`i' = rnormal()
  gen var`=`i'+3' = rbinomial(1,.5)
}

order _all, alphabetic

* You should now have data in which there are three variables 1,2, and 3 which have four different records and four variables var4-var7 which are time invariant.

* We can use expand to reshape the data.

expand 4

* Now we have four instances of each of the original observations.

* We want now to create variables var1 var2 var3 and var4 which represent each different values for each panel period.

* First let's sort and label our different panel periods.

* Clump all of the same id's together
bysort id: gen year=_n

* Now we just copy our variables so that they only occur in the appropriate time period.

* I will do this manually though it would be very easy to do by macros.

gen var1 = var1_1 if year == 1
replace var1 = var1_2 if year == 2
replace var1 = var1_3 if year == 3
replace var1 = var1_4 if year == 4

gen var2 = var2_1 if year == 1
replace var2 = var2_2 if year == 2
replace var2 = var2_3 if year == 3
replace var2 = var2_4 if year == 4

gen var3 = var3_1 if year == 1
replace var3 = var3_2 if year == 2
replace var3 = var3_3 if year == 3
replace var3 = var3_4 if year == 4

* Now we just need to drop the extra variables.

drop var?_?

************************************************
* This can also be accomplished with the reshape command:

clear
set obs 1000

gen id = _n

forv i=1/4 {
  gen var1_`i' = rnormal()
  gen var2_`i' = rnormal()
  gen var3_`i' = rnormal()
  gen var`=`i'+3' = rbinomial(1,.5)
}

order _all, alphabetic

reshape long var1_ var2_ var3_, i(id)

* The tricky thing to remember with reshape is that it requires exact syntax on your variables to be converted.

* That is, reshape is looking for a number at the end of each variable name.  This number it will turn into the j variable.