Estimating Variance as a Function of Treatment Rank Class

Imagine that we have a treatment that we give to five different groups of individuals.  Each individual has a variable response which as a unique mean and variance based on the treatment.  We do not know how the means will change but we believe the variance of responses will expand depending upon what level of treatment the individual gets.  We would like to expressly model both the differences in means and that of the variances.

This code was formulated in response to a question posted on CrossValidated.

We want to solve
$$    \max_{\bf {\hat\beta,\hat\gamma}} (\sum_{i=1}(ln(D(x_i, \hat\mu, \hat\gamma_0+\hat\gamma_1 rank))) $$

# Specify how many individuals are in each of our groups
nobs.group <- 500
 
# Simulate our data
grp1 <- data.frame(values=rnorm(nobs.group,5,1), grp=1)
grp2 <- data.frame(values=rnorm(nobs.group,3,2), grp=2)
grp3 <- data.frame(values=rnorm(nobs.group,6,3), grp=3)
grp4 <- data.frame(values=rnorm(nobs.group,5,4), grp=4)
grp5 <- data.frame(values=rnorm(nobs.group,1,5), grp=5)
 
# Group our data into a single object
mydata <- rbind(grp1,grp2,grp3,grp4,grp5)
 
# Speficy the function to maximize (minimize)
lnp <- function(gamma, x, rank) 
  # I include a negative here because the default option with optim is minimize 
  -sum(log(dnorm(x,gamma[1]*(rank==1)+ 
                   gamma[2]*(rank==2)+
                   gamma[3]*(rank==3)+
                   gamma[4]*(rank==4)+
                   gamma[5]*(rank==5), 
                 gamma[6]+gamma[7]*rank)))
 
ans <- optim(c(
  # Specify initial values for parameters to be estimated
  beta1=1,beta2=1,beta3=1,beta4=1, beta5=1, 
  gamma1=1,gamma2=1), 
  # Specify the function to minimize (maximize)
  lnp, 
  # Input dependent variable as x and the explanatory variable as rank
  x=mydata$values, rank=mydata$grp, 
  # Be sure to inlcude the hessian in the return for 
  # calculating standard errors
  hessian=T)
 
# The standard erros can be estimated using the hessian
stand.error <- sqrt(diag(solve(ans$hessian)))
 
# This will create a nice table of results
cbind(par.est=ans$par, 
      stand.error,
      tstat=ans$par/stand.error,
      pvalue=1-pt(ans$par/stand.error, nrow(mydata)-length(ans$par)),
      CILower=ans$par+stand.error*qt(.05,nrow(mydata)-length(ans$par)),
      CIUpper=ans$par+stand.error*qt(.95,nrow(mydata)-length(ans$par))      
      )

Created by Pretty R at inside-R.org

          par.est stand.error      tstat    pvalue    CILower     CIUpper
beta1   4.9894067  0.04038367 123.550112 0.0000000  4.9229567  5.05585658
beta2   2.0009055  0.10955198  18.264440 0.0000000  1.8206415  2.18116942
beta3   3.7640531  0.19407912  19.394427 0.0000000  3.4447027  4.08340355
beta4   6.4818562  0.23879420  27.144111 0.0000000  6.0889287  6.87478375
beta5  -0.5547626  0.29730735  -1.865957 0.9689176 -1.0439715 -0.06555377
gamma1 -0.4849449  0.05684407  -8.531142 1.0000000 -0.5784798 -0.39140993
gamma2  1.3867000  0.04520519  30.675682 0.0000000  1.3123164  1.46108352

Question and Answer: Generating Binary and Discrete Response Data

I was recently contacted by a reader with two very specific questions and I thought that this would be a good topic to publicity respond to. He would like to simulate his data:
I have firm level data and the model is discrete choice with the main explanatory variable also a binary choice:  First question is how can I calibrate the data generation model? 

Answer: 

This is a fundamental question for any kind of econometric model.  How you calibrate your data implies the inherent structure of your data which in term implies what method you should use to attempt to recover your parameters.  Now some data generating processes exist out there which do not yet have econometric solutions to.  Yet there are many that do.

In general you can calibrate your data by i. modifying the parameters, ii. the distribution of explanatory variables, or iii. the distribution of the errors.

In a binary response case the most common models are probit/logit in which case in order to simulate data you would generate your underlying model and overlay the appropriate CDF over it which gives you probabilities of a success.  Finally you would make a random draw based on those probabilities for each outcome being simulated.

I have numerous example code demonstrating this:
Stata: (Reverse Engineering a Probit) (Probit vs Logit)
 

R:
Nobs <- 10^4
X <- cbind(cons=1, X1=rnorm(Nobs),X2=rnorm(Nobs),X3=rnorm(Nobs))
B <- c(B0=-.2, B1=-.1,B2=0,B3=-.2)
P <- pnorm(X%*%B)
SData <- as.data.frame(cbind(Y=rbinom(Nobs,1,P), X))
summary(glm(Y ~ X1 + X2 + X3, family = binomial(link = "probit"), data = SData))

Discrete Data
As for discrete data, it is less clear what the optimal choice is. I prefer the logistic regression which is basically an extension of the Logit model with a few interesting caveats.

Stata: (Simulating Multinomial Logit)
R: (here is an article dealing specifically with using R to create discrete response data http://works.bepress.com/joseph_hilbe/3/)

Nobs <- 10^4
X <- cbind(cons=1, X1=rnorm(Nobs),X2=rnorm(Nobs),X3=rnorm(Nobs))
# Coefficients, each input vector (c) is associated with a different outcome
B <- cbind(0, c(B0=-.2, B1=-.1,B2=0,B3=-.2), c(B0=.3, B1=0,B2=.6,B3=.4))
# Everything is relative to option 1 which is the default
num <- exp(X%*%B) # Numerator
den <- apply(num,1,sum) # Denominator
P <- num * 1/cbind(den,den,den) # Probability
CP <- cbind(P[,1],P[,1]+P[,2]) # Cumulative probabilities
U <- runif(Nobs) # Draw from the uniform draw
Y <- rep(0,Nobs) ; Y[U>CP[,1]]<-1; Y[U>CP[,2]]<-2 # Calculate outcome

SData <- as.data.frame(cbind(Y=Y, X)) # Combine Datarequire("nnet")
summary(Mlogit <- multinom(Y ~ X1 + X2 + X3, data = SData))

Power Analysis by Simulation: R, RCT, Malaria Example

I have received a number of requests for demonstration code on how to perform a power analysis using simulation in R.  I have already demonstrated howto do this in Stata but lacked the easy to use Stata command “simulate” that I preferred.  However, in a recent post I have written up a command very similar to simulate in R called SimpleSim.

This command takes is capable of taking a vector of parameters and feeding them into a function which returns a vector of results.  The results are turned into a table which can be used to show the rates of type 1 and type 2 error given a particular sample size and underlying effect.

Example:
Let's imagine that we are interested in testing the effect that bednet usage has on the rate of infection of malaria.  We are concerned that the decision to purchase and use bednets might endogenous as a result of those more likely to come down with malaria being more prone to seeing bednets as a valuable investment.

So, in response we are going to design an randomized control trial (RCT) in which we either supply fully paid bednet coupons or half discounted coupons each to different thirds of the population leaving the remaining third as a control.  We are concerned that if supply the bednets at no cost then the recipients will not value them so that is why we give coupons that only reduce the cost of the bednets by 1/2.

After distributing the coupons, we would like to use the random distribution of the coupons as an instrument for use of a bednet.  Ultimately we would like to do two levels of analysis.

1. Using instrumental variables (IV) to see how effective bednets are at preventing Malaria in an active population. This is actually an interesting question because often there is imperfect compliance because few individuals spend all of the hours at dusk and dawn under a bednet.  In this case I should use an estimator that takes the form of the first stage being a binary regressor followed by a second stage binary regressor.  However, since linear models seem to be pretty good at approximating average partial effects I am just going to use a IV estimator.
2. Using ordinary least squares (OLS) we would like to see how effective each of the coupon programs are at increasing use of bednets within the homes.

require("AER") # We will use the ivreg function later

 
# We define a function that both simulates a population and estimates 
# results and returns and returns the results.
MalariaIV <- function(
  nsim = 100,
  npop = 1000,    # Define the sampling population
  treat0 = 1/3,   # Define the proportion which is control
  treat1 = 1/3,   # Define the proportion which is treated with free bednets
  treat2 = 1/3,   # Define the proportion which is treated with 50% cost
  t0comp = .05,   # Bednet useage among control group
  t1comp = .85,   # Bednet useage if recieving free net
  t2comp = .5 ,   # Bednet useage if recieving 50% cost
  malariaRT = .85,# Rate of getting malaria without bednet
  netRT = .45,    # Rate of getting malaria with bednet
  Tdetect = 1,    # Likihood of detecting malaria if present
  Fdetect = 0,    # Likihood of detecting malaria if not present
  alpha = .05     # The alpha level that a p-value must be below
) {
  # Define how the function works
  # First, define a vector to store results
  pvalues <- NULL
 
  #
  for (i in 1:nsim) { # Repeat the simulation a number of times
    # Generate the population by technology useage
    simdata <- data.frame(
      control=c(rep(1,npop*treat0),rep(0,npop*treat1), rep(0,npop*treat2)),
      treat1 =c(rep(0,npop*treat0),rep(1,npop*treat1), rep(0,npop*treat2)),
      treat2 =c(rep(0,npop*treat0),rep(0,npop*treat1), rep(1,npop*treat2)))
    # Calculate the actual population generated (should be 999 in this case)
    npeople <- nrow(simdata)
    # Generate the rate of bednet usage.
    simdata$bednet <-
      simdata$control*rbinom(npeople,1,t0comp) + # Bednet usage control
      simdata$treat1 *rbinom(npeople,1,t1comp) + # Free
      simdata$treat2 *rbinom(npeople,1,t2comp)   # 50% cost
    # Now let's generate the rate of malaria
    simdata$malaria <-
      (simdata$bednet==0)*rbinom(npeople,1,malariaRT)+
      (simdata$bednet==1)*rbinom(npeople,1,netRT)
    # Finally generate the rate of malaria detection as a function
    # of true parasite levels.
    simdata$mdetect <- 
      (simdata$malaria==0)*rbinom(npeople,1,Fdetect)+ # False detection
      (simdata$malaria==1)*rbinom(npeople,1,Tdetect)  # True detection
 
    # Time to do our simple estimation of the effect of treatment on bednet use
    lmcoef <- summary(lm(bednet~treat1+treat2,data=simdata))$coefficients
    # Now let's try the 2SLS to estimate the effect of bednets on contraction
    # of malaria.
 
    ivregest <- ivreg(mdetect~bednet | treat1+treat2, data=simdata)
    ivregcoef <- summary(ivregest)$coefficients
 
    # Save the rejection of the null rates
    pvalues <- rbind(pvalues,c(treat1=lmcoef[2,4]<alpha,
                               treat2=lmcoef[3,4]<alpha,
                               iv=ivregcoef[2,4]<alpha))
 
  }
  # Calculate the mean rejection rate for each coefficient.
  apply(pvalues,2,mean)
}
 
MalariaIV() 
# Running it once we can see that we get a single set of results
# where we easily reject the null.  However, we want to know
# what happens when bednets are not so effective or malaria is harder
# to detect.  We can modify the parameters fed into the model to test
# these questions manually or we could use the SimpleSim function from
# a previous post.
 
SimpleSim(fun=MalariaIV, 
          npop=c(100,1000),
          t0comp=c(.05,.25),
          t1comp=c(.5,.75),
          t2comp=c(.25,.5),
          malariaRT=c(.85,.5),
          netRT=c(.75,.45),
          Tdetect=c(1,.7),
          Fdetect=c(0,.3),
          alpha=.05,
          nsim=10)
# This could take a little while to run since there are 256 combinations
# to try and each of them will be run 10 times.
 
# The above command gives back lots of data but it is not always very easy
# to understand in a matrix form.  It is often easier to just vary one
# paramter at a time.
 
sample.size <- SimpleSim(fun=MalariaIV, 
          npop=c(500,1000*1:10),
          t0comp=c(.25),
          t1comp=c(.75),
          t2comp=c(.375),
          malariaRT=c(.5),
          netRT=c(.35),
          Tdetect=c(.8),
          Fdetect=c(.2),
          alpha=.05,
          nsim=200)
 
# Looking at just sample size
require(ggplot2)
 
# Save the results to single long data format to be useable by ggplot2
results <- with(sample.size, data.frame(reject = as.numeric(
  c(iv,treat2,treat1)),
  id=rep(c("iv","treat2","treat1"), each=length(npop)),
  npop))
 
p <- ggplot(results, aes(npop, reject))
# Normally a 80% detection rate is the minimum rejection rate needed 
# to justify a study. 
p + geom_point(aes(colour =id)) + 
  geom_line(aes(group=id)) + 
  geom_hline(yintercept = .8) 
 

# Thus we want to ensure our study has at least 5000 participants.
 
# I think there might be some additional considerations for the
# number of participants required to ensure that there is no
# false rejection of the null.  However, I am no expert on Power 
# Analysis so this is what I know.

Formatted by Pretty R at inside-R.org

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.

Non-Parametric Regression Discontinuity

* I recently went to an interesting seminar today by Matias Cattaneo from the University of Michigan.

* He was presenting some of his work on non-parametric regression discontinuity design which I found interesting.

* What he was working on and the conclusions of the paper was interesting but even more interesting was a release by him and coauthors of a Stata package that implements RD design for easy

* net install rdrobust, from(http://www-personal.umich.edu/~cattaneo/rdrobust) replace

* Regression discontinuity is a technique which allows identification of a localized effect around a natural or structured policy discontinuity.

* For instance, if you wondering what the effect federal grants have on college attendance, then you may be concerned that just looking at those students who are eligible for federal grants in contrast with those who are not eligible will be problematic because students who are eligible (low income) may different than those who are not eligible for the grant (not low income).

* The RD argument is that if individuals do not, as a response to the grant being available, move their reported income level to become eligible for the grant than those who are near the cut off for the grant and those not near the cut off will be fundamentally very similar.

* This may occur if for instance the income cut-off for the grant is unknown.

* So even if students are systematically under-reporting their income, they are not doing it aware of the actual cut off, so the students sufficiently close, above and below the cut off are arguably the "same" or drawn from the same pool except that one group received the program and another group did not.

* The previous post deals some with assuming a linear structure of the underlying characteristics.

* http://www.econometricsbysimulation.com/2012/08/sharp-regression-discontinuity-example.html

* However, the more interesting case (potentially) may be when we assume a nonlinear response to the income in our dependent variable.

* But before going there let's think about what this method boils down to.

* Like all identification methods in statistics or econometrics when we do not have experimental data, identification of an effect is driven by some exogeneity argument.

* That is, x causes y and is unrelated to u (the error).  In the case when u may be correlated with the error the use an exogenous variable to force the movement in the variable of interest may be sufficient to identify a causal effect.

* In this case, clearly it is not enough to simply see what the average y response (GPA, attendance, graduation rates, whatever) is to a change in grant level because those who receive the grants are systematically different from those who do not.

* However, because the cut off for receiving the grant is unknown, around the cut off the two samples who receive the grant and who do not can arguably be considered the same.

* Thus, we could say that the unknown position of the cut off is the random exogenous variable which near the cut off forces some students into the group that receives the grant and some students into the group that does not.

* Let's imagine some non-parametric relationship between income and performance:

clear

set obs 10000

gen income = 3^((runiform()-.75)*4)
  label var income "Reported Income"

  sum income
gen perf0 = ln(income) + sin((income-r(min))/r(max)*4*_pi)/3 + 3
  label var perf0 "Performance Index - Base"

scatter perf0 income

* Looks pretty non-parametric

* Let's add in some random noise
gen perf1 = perf0 + rnormal()*.5
  label var perf1 "Performance Index - with noise"

scatter  perf1 income

* Using the user written command rcspline, we can see the local average performance as a function of income.

* ssc install rcspline

rcspline perf1 income,  nknots(7) showknots title(Cubic Spline)
* I specify "7" knots which are the maximum allowed in the rcspline command.

* The spline seems to fit the generated data well.

* Now let's add a discontinuity at .5.

gen grant = income&lt;.5
sum grant

* So about 50% of our sample is eligible for the grant.

* Now let's add the grant effect.

* First let's generate an income variable that is centered at the grant cut point.
gen income_center = income-.5

gen perf2 = perf1 + .5*grant - .1*income_center*grant
  * Thus the grant is more effective for students with lower income.
  label var perf2 "Observed Performance"

**** Simulation done: Estimation Start ****

rcspline perf2 income,  knots(.15 .25 .35 .37 .4 .45 .5 .55 .6 .65 .75 .85 1.1 1.25 1.5) title(Cubic Spline)
* This is obviously not the ideal plot and I have had some difficulty finding a command which will generate the plot that I would like.

* However, we can see that there does appear to be "something" going on.

reg perf2 income grant
* We can see that our itial estimate of the effect of the grant is entirely wrong.

* It appears so far that the effect of the grant on performance is actually hindering performance (which we know is false).

* Now, let's try our new command rdrobust

rdrobust perf2 income_center
* The default cut point is at 0.  Thus using income_centered works.

* Though this estimate is negative and thus seems the reverse of what we would expect, it is actually working quite well.

* That is because regression discontinuity is trying to identify the effect of the discontinuity on the outcome variable with the default assumption that at the discontinuity the forcing variable is becoming 1.

* In this case however, the discontinuity is really driving the grant to be equal to zero.

* Thus we must inverse the sign on the rd estimator in order to identify the true effect in this case.

* Alternatively, we could switch the sign of income.

gen nincome_center = income_center*(-1)

rdrobust perf2 nincome_center

* rdrobust is a newly designed command that has some extra bells and whistles that other regression discontinuity commands have as well as some oddities.

* I would suggest also looking to the more official stata command rd (ssc install rd)
rd perf2 nincome_center

* This command is nice because it estimates many bandwidths through the mbw option.

* The default mbw is "100 50 200" which means, use the 100 MSE (mean squared error) minimizing bandwidth, half of it and twice it.

* We can plot our estimates of the treatment effect using a range of bandwidths.

gen effect_est = .
  label var effect_est "Estimated Effect"

gen band_scale = .
  label var band_scale "Bandwidth as a Scale Factor of Bandwidth that Minimizes MSE"


forv i = 1/16 {
  rd perf2 nincome_center, mbw(100 `=`i'*25')
    if `i' ~= 4 replace effect_est = _b[lwald`=`i'*25'] if _n==`i'
    if `i' == 4 replace effect_est = _b[lwald] if _n==`i'
    replace band_scale = `=`i'*25'     if _n==`i'  
}
gen true_effect = .5
  label var true_effect "True effect"

two (scatter effect_est band_scale) (line true_effect band_scale)

* We can see around the 100% MSE bandwidth estimates are fairly steady though they dip a tiny bit.

Path Analysis

Stata do file

* Path analysis is an interesting statistical method that can be used to indentify complex relationships beween variables and an outcome variable.

* As with all statistical methods the modelling framework is essential to derive reasonable results.

* Conviently, I am only interested in simulating data so as usual my data will perfectly conform to the model's specifications.

* Imagine the following model.

* All of the boxes are observable variables.  The arrows indicate the causal direction of the effects.

* There are two exogenous variables: A and D.   These variables are not influenced by any other variables in the model.

* All other variables are endogenous.

* Each of the variables represents a direct effect of a one unit change in one variable on that of the other variable.

* This framework is convient because it allows us to indentify a "total effect" which is a combined result of both the direct and indirect effects of variables on the outcome variable.

* The variable of primary interest in explaining is H.

* The variable G has only a direct effect on H (pHG).

* While the variable C only has an indirect effect on H (pFC*pHF).

* The reason the indirect effect is a product is because C has a pFC effect on F, and F has a pHF effect on H, thus a change in H as a result of a change in C is how much F changes as a result of C and how much that change effects H.

* Variables can have both and indirect and direct effect.

* B for instance has the direct effect: pHB
* Indirect effects: pCB*pFC*pHF + pEB*pHE
* Total effect: pHB + pCB*pFC*pHF + pEB*pHE

* The key feature about this particular example is that all of the arrows are one directional.

* Making a great deal of inference possible that otherwise would not be possible.

* Usually we cannot say that when trying to explain H with explanatory variables A through G that A causes B and B causes H.

* However, if we do the work to indentify reasonable pathways then this type of analysis could be quite interesting.

* Let's generate out data.

clear

* Let's imagine 6000 youth in our sample.

set obs 6000

* Let's first specify our effects

* pEA = .3
* pEB = .13

* pHA = .2
* pHB = .2
* pHE = .3
* pHG = 1.1
* pHF = .2

* pBA = .5

* pCB = .2
* pCD = .1

* pGD = .2

* pFC = .76
* pFB = .4

* For B we can calculate our true effects:

* B Direct: pHB = .2

* Indirect effects: pCB*pFC*pHF + pEB*pHE
* Indirect effects:.2*.76*.2 + .13*.3 = .0694

* Total effect: .0694+.2 = .2694

gen A =                rnormal()
gen B = A*.5 +         rnormal()
gen D =                rnormal()
gen C = B*.2 + D*.1 +  rnormal()
gen E = A*.3 + B*.13 + rnormal()
gen F = B*.4 + C*.76 + rnormal()
gen G = D*.2 +         rnormal()
gen H = E*.3 + A*.2 + B*.2 + F*.2 + G*1.1 + rnormal()

* Simualtion Done

* In order to generate our different effects we simply run OLS for each endogenous variable.

reg A B
  local pBA = _b[B]

reg C B D
  local pCB = _b[B]
  local pCD = _b[D]

reg C B D
  local pCB = _b[B]
  local pCD = _b[D]

reg G D
  local pGD = _b[D]

reg F C B
  local pFB = _b[B]
  local pFC = _b[C]

reg E A B
  local pEA = _b[A]
  local pEB = _b[B]

reg H A B E F G
  local pHA = _b[A]
  local pHB = _b[B]
  local pHE = _b[E]
  local pHF = _b[F]
  local pHG = _b[G]

* In order to estimate the indirect effect say of B on H.
* We just plug our estimates into the equation.

* B direct effect: pHB
* Indirect effects: pCB*pFC*pHF + pEB*pHE
* Total effect: pHB + pCB*pFC*pHF + pEB*pHE

di "B's estimated indirect effect = `pCB'*`pFC'*`pHF' + `pEB'*`pHE'"
di "B's estimated indirect effect = " `pCB'*`pFC'*`pHF' + `pEB'*`pHE'

* Which turns out to be close to our true value.

di "B's total estimated effect on H is " `pHB' + `pCB'*`pFC'*`pHF' + `pEB'*`pHE'

* It is possible to use the user written command pathreg to make things easier.

* Install it by typing the following command. findit pathreg
pathreg (H E B F G) (G D) (C B D) (B A) (E A B) (F B C)

* This command does not currently calculate out all of the indirect and direct effects.

* I am not sure the best way to calculate the standard errors of the different effect estimates.

* My guess is that since this is just a series of fast OLS regressions the easiest thing to do would be to boostrap the entire process.

* This would require slightly more code but definitely easy to do from this point.

R-squared

Original Code

* R-Squared

* R-squared and pseudo r-squared is a useful statistics produced by most regression type estimation routines.

* R-squared (R2) is a measure of how much of the variance in y is explained by the model.

* Thus a model with only an intercept has an R2 of 0.

set seed 101
clear
set obs 10000

gen y1=rnormal()

reg y1

* While in the opposite extreme a model which does not have any unexplained variance has an r2 of 1.

reg y1 y1, noconstant

* Technically this regression should not work but Stata does the math and produces the results.

* Let's see how well R2 approximates explainable variation.

gen x=rnormal()
gen u=rnormal()

gen y2 = (1)*x + u
* The variance of the model is equal to 1 (from 1^2 * var(x))
* The variance of the unexplained error is equal to 1 (var(u))
* Thus our true explained variance should be equal to var(x)/(1^2 * var(x)+var(u)) = 1/2

reg y2 x
* Thus we can see our R2 estimate of the explained variance is very close to the true which is .5

* I made enphasis on noting the coefficient on the x.

* That coefficient significantly scales explainable variation.

* Thus:
gen y3 = 2*x + u

* Should have a much larger R2 because model variance = 2^2*varx = 4
* Var(u) = 1
* R2 = 4/(4+1)=.8

reg y3 x

* If we were to add multiple xs the calculation is similar though if there is correlation between the xs then that will factor into the model.

gen x1 = rnormal()
gen x2 = rnormal()

gen y4 = x1 + x2 + u

* R2 = var(x1) + var(x2) / (var(x1) + var(x2) + var(u) = 2/3
reg y4 x1 x2

* If there is correlation between the xs then that can substantially throw off the calculations.

* In the extreme cases corr(x1, x2)=1 then we are back to the same scenario as y3

* y = x1 + x2 + u (if x1~N(0,1) and x2~N(0,1)) then x1=x2

* y = 2*x1 + u

* Thus R2 = .8

* In the other extreme corr(x1,x2)=-1

* Then, given that they are both N~(0,1), x2=-x1

* y = x1 + x2 + u = x1 - x1 + u = u

* Which is the same as y1

* R2 = 0

* R-squared can also be thought of as the square of the correlation between the predicted values and the observed.

reg y4 x1 x2
predict y4hat

corr y4hat y4
* Thus we can see that there is an 81% correlation between yhat and y observed.

* A high correlation would indicate that our model have done well at predicting observable characteristics.

di r(rho)^2

* A brief note on adjusted R2.

* R2 is known to always be larger the more variables are in your model.


gen z1 = rnormal()
gen z2 = rnormal()
gen z3 = rnormal()

reg y4 x? z?

* Thus: the R2 moved from   R-squared     =  0.6610 to
*                           R-squared     =  0.6611

* This factor being known researchers have developed the Adj-R2 which slightly penalizes the R2 for including more variables.

* Thus Adj R-squared =  0.6609

* This might be appropriate given known facts, however it is trivial and almost always worth ignoring.

* I generally don't pay attention to the AR2 and I don't know anybody else who does either.

* A .0001 difference in R2 is so unimportant as to be completely ignorable without significant loss of content.

Estimating Random Coefficients on X

Original Code

* A frequent assumption in economics is that the coefficient that is being estimated is constant throughout the population.
* That is, the effect of the exogenous variable is the same for all individuals.
* This assumption is implicit in the notation:

* Y = XB + u
* Because B is assumed to be constant for all individuals.
* However, let's for one second imagine that we are not estimating the constant B, but rather a random variable b.
* Thus: Y = Xb + u

* This immediately presents an obvious problem.

* What do we want to know from this formulation?
* One problem is that for however many observations we have, we have the same number of bs, thus we cannot hope to estimate the individual b values.

* Where most people have taken this model is to say we are interested in two primary things:
* 1. What is the average effect of X on Y?  Ie. E(dY/dx|X) which is really what we are after assuming constant coefficient.
* 2. And how much variance exists in b?

* To think about this the following formulation is useful:
* b = B + v
* Where B is the average cofficient and v is the random component.

* Now to insert this into the estimation equation:

* Y = Xb + u = X(B + v) + u = XB + Xv + u = XB +  e
* Now, we can easily show that given E(v|x)=E(e|X)=0 OLS is unbiased.
* Bhat = (X'X)^-1 X'Y = (X'X)^-1 X'(XB + e) = (X'X)^-1 X'XB + (X'X)^-1 X'e
* Bhat = B + (X'X)^-1 X'e = B + (X'X)^-1 X'(Xv + u) = B + v + (X'X)^-1 X'u
* E(Bhat|x) = E(B + v + (X'X)^-1 X'u|x) = B + E(v|X) +  (X'X)^-1 X'E(u|X) = B

* However, OLS is no longer the most efficient estimator of B any longer (because the homoskedasticity assumption is violated).

clear
set obs 10000
set seed 121

gen x = runiform()*5

* I want v and u be drawn from a multivariate normal which allows for them to be correlated.
matrix C = (9, 2.5 \ 2.5, 23)

drawnorm v u, cov(C)

gen y = 4*x + v*x + u

* We can see OLS works "fine"
reg y x

predict uhat, resid
scatter uhat x, sort title(Heteroskedastic Normally Distributed Error)
* This is a form of heteroskedasticity.

* So, we would like to do both 1 and 2 above.
* We would like to estimate not only the average effect of x on y (what we can accomplish with OLS),
* but also the variance of the error terms which is u and which is v.


* In order to do that we will first define the NormalReg.
* That is the regression MLE form assuming the errors are normally distributed.

cap program drop myNormalReg
program define myNormalReg

  * The first argument of any maximum likelihood program is the name of the temporary log likelihood variable created by the stata when the ml procedure is called.
  * Each additional argument is an linear "equation" that Stata maximizes by choosing estimators for.
  args lnlk xb sigma2
  * Thus we have two equations that we are asking Stata to solve.

  * The following is the log likelihood value that we would like Stata to maximize.
  * In this case it is the log of the normal density asking Stata to choose the optimal mean and standard deviation.
  qui replace `lnlk' = -ln(sqrt(`sigma2') * sqrt(2*_pi)) - ($ML_y-`xb')^2/(2*`sigma2')
  * I have opted to program the pdf of the normal in the equation rather than the build in normal pdf command.

  * Notice that unlike my previous post on modelling heteroskedasticity (http://www.econometricsbysimulation.com/2012/11/modeling-heteroskedasticity.html)
  * I am using sigma squared rather than sigma as the primarily form of the variance to be modelled.
  * This choice is not trivial.
end

* The first thing you are probably wondering is "why am I using a model that only has one error term when I really want to learn about two error terms?"

* The answer is that, as seen above, the sum of the error v and u can be expressed in a single form e.

* Thus the variance of e equals: var(e) = var(u) + Var(v)*x^2 + 2*cov(v,u)*x

* This conviently is a linear form that can easily be included in the above equation.

* First we need to generate x2
gen x2 = x^2
gen varu = 1

ml model lf myNormalReg (y = x) (sigma2: varu x2 x, noconstant)
ml maximize

matrix list C

* Looking back at matrix C we can see that our coefficient on varu corresponds well with the variance of u.
* Our coefficient on x2 approximates our variance of v.
* But our coefficient on x is way too large for our covariance terms.
* This is expected, the coefficient on x is supposed to be twice as large as the cov(v,u) from the formulation of the var(e) term.

local varu = [sigma2]_b[varu]
local varv = [sigma2]_b[x2]
local cov_uv = [sigma2]_b[x]/2

matrix Chat = (`varv', `cov_uv' \ `cov_uv', `varu')

matrix list Chat
matrix list C

* In order for the MLE estimator to be efficient in this case we are assuming that both v and u are normally distributed (which they are because we generated them).
* We may also be concerned that the generated covariance maxtrix Chat is not Positive Semi-Definite, a requirement for v and u to be jointly multivariate normally distributed.

* We can check this by attempting to take the inverse:
matrix INVA = inv(Chat)
matrix list INVA

* In this case there was not problem inverting Chat.  But this need not be the case.

* Let's try simulating some outcomes for a matrix a hairs breathe away from not being PSD.

cap program drop PDcheck
program define PDcheck, rclass
  clear
  set obs 1000

  gen x = runiform()*5
  matrix C = (2, 4 \ 4, 11)

  drawnorm v u, cov(C)

  gen y = 4*x + v*x + u
  gen x2 = x^2
  gen varu = 1

  ml model lf `1' (y = x) (sigma2: varu x2 x, noconstant)
  ml maximize

  return scalar beta_coef = [eq1]_b[x]

  local varu = [sigma2]_b[varu]
    return scalar varu = `varu'
  local varv = [sigma2]_b[x2]
    return scalar varv = `varv'
  local cov_uv = [sigma2]_b[x]/2
    return scalar cov_uv = `cov_uv'

  local PD = 0
  if (`varu'*`varv' - `cov_uv'^2 >0)&(`varu'>0)&(`varv'>0) local PD = 1
  return scalar PD=`PD'

end

PDcheck myNormalReg
return list

simulate PD=r(PD) beta_coef=r(beta_coef) varu=r(varu) ///
         varv=r(varv) cov_uv=r(cov_uv), rep(50): PDcheck myNormalReg
sum
* We can see that in this formulation only about half the time is the covariance matrix PD.

* This could present a problem.
* There are various ways of trying to correct for this.
* All of them require more advanced programming than I currently understand using the MLE syntax.
* I hope/need to learn how to make such corrections in the near future.

* Note, Stata is already doing some clever behind the scenes manipulations to make sure only feasible values of the parameters are chosen.

* Stata must ensure that sigma >= 0 for the psd (because it is taking the square root).



cap program drop myNormalReg2
program define myNormalReg2

  * The first argument of any maximum likelihood program is the name of the temporary log likelihood variable created by the stata when the ml procedure is called.
  * Each additional argument is an linear "equation" that Stata maximizes by choosing estimators for.
  args lnlk xb var_u var_v cov_uv

  tempvar sigma2 psd_check
  gen `sigma2' = (exp(`var_u') + exp(`var_v') + 2*`cov_uv')

  gen `psd_check' = (`var_u'*`var_v' - `cov_uv'^2)

  qui replace `lnlk' = -ln(sqrt(`sigma2') * sqrt(2*_pi)) - ($ML_y-`xb')^2/(2*`sigma2') * ///
                        (1^sqrt(`psd_check')) * (1^sqrt(`var_u'))
  * I have opted to program the pdf of the normal in the equation rather than the build in normal pdf command.

  * Notice that unlike my previous post on modelling heteroskedasticity (http://www.econometricsbysimulation.com/2012/11/modeling-heteroskedasticity.html)
  * I am using sigma squared rather than sigma as the primarily form of the variance to be modelled.
  * This choice is not trivial.
end

  ml model lf myNormalReg2 (y = x) (var_u: ) (var_v: x2, noconstant) (cov_uv: x, noconstant)
  ml maximize
  *  Note that because I was able to specify simga directly I have already multiplied by the constant 2 causing the cov estimate to be scaled properly.

  * This new specification is guaranteed not to create covariance matrices that are not positive definite.
  * However, it frequently fails to converge which means that it is really not a very good algorithm.
  * I will continue to experiment with this problem in the near future.

  * This next step is using a Cholesky Decomposition to pre-specify the variance matrix (as suggested by my advisor Jeff Wooldridge).

  * I will post more on that as I make steps forward.

  * Also, if anybody knows any built in commands in Stata that can do this, please post them.

  * I have used xtmixed previously to identify random coefficients.

  * However, I think the command imposes orthogonality on the random coefficients.

Simulating Spatial Data

# Spatial data tags are an increasingly recorded for data that is being generating as a result of widescale implementation of GPS technology.

# In this post I will present a simulation in which the population is distributed around a single town center.

# The important characteristic of the population is that their characteristics are autocorrelated.  This can either be in response values or in error terms.

# Let's first specify an imaginary population.

population = 1000

# Specify the location in which the population is centered

center.X = 0
center.Y = 0

# Dispersion, the average distance of an individual from the center

dispersion.km = 40

# Now let's generate our population position

# First off.  We would like our population to be distributed in a circle around the center.  So we need to pick an angle in radians.

angle <- pi*runif(population)*2

# This is the distance from the center that each member of the population will be.  If this were set to a constant then the following commands would end up drawing a circle.
distance <- rnorm(population)*dispersion.km

mydata <- data.frame(x=center.X+distance*cos(angle), y=center.Y+distance*sin(angle))

smoothScatter(mydata, nrpoints=0, main="Population Density")

# Let's imagine that there is an unobserved variable called "soil quality" which varies in both the x and y.

mydata$soil.quality <- sin(mydata$x*pi/50+mydata$y*pi/100)+1

# Farm size is random. But on average smaller as the plots get closer to the city center.

mydata$size <- runif(population)/4+abs(distance)/100

# In order to produce some cool graphs we will need to install a new package:
install.packages("ggplot2")
library(ggplot2)

qplot(x, y, data=mydata, size=size, colour = soil.quality, main="Farms tend to be smaller near town")

# Rainfall also is spatially correlated.
mydata$rainfall <- sin(mydata$x*pi/60+mydata$y*pi/160) + sin(mydata$y*pi/60)+1

qplot(x, y, data=mydata, size=size, colour = rainfall, main="Rainfall is also distributed spatially")

# Now, let's imagine some technology usage, say fertilizer.

mydata$fert.use = mydata$rainfall-mydata$soil.quality+mydata$size+rnorm(population) + 4

qplot(x, y, data=mydata, size=size, colour = fert.use, main="Fertilizer use as a result should also be spatially distributed")

# Now let's see if we can't test if we can if fertilizer use is spatially correlated.

# The trick is figuring out what that means.

# I will define it as this, spatial correlation is the test to see if the use of fertilizer by one person is correlated with the use of fertilizer by another person.

# So I need to figure out a way of finding out what the closest person fertilizer usage is.

neighbor.fert.use <- neighbor.x <- neighbor.y <- 0

for (i in 1:population) {
  # We will look at each person i and find the person which is closest.
  # First let us constuct a variable that measures distance.
  # This is that standard Euclidean distance formula.
  distance.from.i = ((mydata$x[i]-mydata$x)^2 + (mydata$y[i]-mydata$y)^2)^.5

  # The following set of nested statements can be somewhat confusing.  Read from the inside statement first.
  # rank() will create a vector of length population that ranks distance from i.
  # Outside of that is a logical operator that will create a vector of length population which is all False except from rank == 2 which is true which means that mydata$fert will be drawn from that single value of rank==2.
  # Finally we assign it to the ith place in the neighbor.fert.use vector.
  neighbor.fert.use[i]<- mydata$fert.use[rank(distance.from.i, ties.method="random")==2]
  neighbor.x[i]<- mydata$x[rank(distance.from.i, ties.method="random")==2]
  neighbor.y[i]<- mydata$y[rank(distance.from.i, ties.method="random")==2]
}

plot(mydata$x,mydata$y, main="Arrows indicate closest farm")
for (i in 1:population) arrows(x0=mydata$x[i], x1=neighbor.x[i],y0=mydata$y[i], y1=neighbor.y[i],length = 0.075)

# From the plot we can see that the farm matching algorithm specified above appears to be working well.  We can see from that plot that every farm has a farm that is closest to it however there is some farms that are not the closest farm to any other farm (which makes sense).

cor(neighbor.fert.use, mydata$fert.use)
# From the correlation between fertilizer uses between each farm and it's closest neighbor I get a correlation of around .4
# This indicates that fertilizer use is spatially correlated.  However, as a result of how we set up the model this correlation is rather simple to handle.  It is not because a neighbor uses fertilizer that a farmer will use fertilizer but rather the effect of unobservables which is driving the use of fertilizer.  Mainly rainfall variation and soil quality variation.

# Imagine that we observe rainfall but not soil quality.  Let's see how well we can predict fertilizer usage.

# first let's remember how fertilizer use is calculated:
# mydata$fert.use = mydata$rainfall-mydata$soil.quality+mydata$size+rnorm(population) + 4

summary(lm(fert.use~rainfall+size,data=mydata))
# We can see that the coefficient on rainfall is too small.  This is because within the construction of this sample rainfall is correlated with soil quality.
cor(mydata$rainfall,mydata$soil.quality)

# Let's try to include technology choice of closest farm neihbor as a proxy for the spatial correlation of soil quality.
summary(lm(fert.use~rainfall+size+neighbor.fert.use,data=mydata))

# We can see that the coefficient on rainfall is even smaller.  This indicates that controlling for the neighbors choice is not helping.  Why is that?

# Sure, controlling for the neighbors fertilizer use is controlling for some of the soil quality variable.  However, what else is in the neighbor's decision = rainfail+size.  Since size is correlated spatially and rainfall is correlated spatially controlling for the neighbor's choice in effect controls for some of the effects of the explanatory variables.  Thus, both rain and size variables suffer as a result of controlling for the nearest neighbor's technology choice.