Matrix operations in R

# There are many ways of inputing matrices into R

# Cbind will bind vectors or other matrices together by adding on the columns of one to the other.
A = cbind(c(2,3),c(3,2),c(1,3))

# by the way, the c() function binds a set of elements seperated by commas into a vector.

A

# Rbind does the same as cbind but uses rows instead
B = cbind(c(-1, 1, 4), c(0, 1, 5), c(3, -1, 1))

B

# We can also create a matrix by using the matrix command.
# Data input in this manner is read from a vector into columns.

C = matrix(c(1,4,3,2),nrow = 2, ncol=2)

C

# An array can accompish the same effect as a matrix.

# The largest difference is that an array can take on more than two dimensions.

D = array(c(10,4,5,2), dim=c(2,2))

D

# One need not populate an matrix/array at creation when specifying matrix or array.

E = matrix( NA, nrow = 4, ncol = 4)

E
# We can see the matrix is filled with empty values

# We can replace individual elements by specifying their positions
E[1,1] = 3
E[2,1] = 0
E[3,1] = 4
E[4,1] = -1

# As well as entire columns or rows
E[,2] = c(0,2,3,2)

# Or subsection of the matrix with another matrix
E[,3:4] = cbind(c(2,3,2,1),c(-1,2,1,0))

E

# Sometimes you do not need to specify every element of the matrix if there are common elements
F = matrix(0, nrow=3, ncol=3)

F[1,1] = 4
F[2,2] = 3
F[3,3] = 6

F
# F is a diagnol matrix which is populated first with 0s then filled with the diagnol values.

# Sometimes we start with a data set and want to convert it to a matrix
G = data.frame(score1=c(5,10,-7),score2=c(-1,17,3),score3=c(0,5,2))

G

# Now let's convert it to a matrix
G = data.matrix(G)

G

# Now let's do some matrix operations with the matrices we have defined:

# a. C x D

C
D

C %*% D

# Which is different from element wise multiplication which is

C * D

# b. A x B

A
B

A %*%B

# c. B x A

# d. E x E' = E times the transpose of E

E
E %*% t(E)

# e. F^-1 or F inverse.  We know the inverse of F exists because it is a diagnol matrix with poisitive diagnol non-zero elements.

solve(F)

require(MASS)
ginv(F)

# Neither commands work for this particular application.

# However, it is easy to see that in a diagnol matrix the inverse is just.

# First we will make Finv the same as F
Finv = F

# This is an interesting bit of R functionality

# On the right the diag command is retrieving the diagnol of F and doing a element wise 1/x operation.

# The diag on the left is specifying target values to be replaced.
diag(Finv) = 1/diag(F)

# This might seem a little odd.  Let's try this kind of thing on Z
Z <- Finv

Finv

diag(Z) <- 23
Z
# In this case because 23 is a single number that can be duplicated throughout the diagnol of Z, there is no problem.


# We can check that this is really the inverse of F

Finv %*% F

# Or equally good

F %*% Finv

# Interestingly element by element multiplation yeilds the same result in this example

F * Finv

# f. Rank(D)

D
qr(D)$rank

# F only has rank 1.  This can be seen by dividing col 1 by col 2.

D[,1]/D[,2]

# Both numbers are two meaning item [1,1] = [1,2]*2 and [2,1]=[2,2]*2

# Interestingly this trick works for rows as well:

D[1,]/D[2,]

# g. B + G

B
G

B + G

# h. B - G

B - G

# Addition and subtraction is element by element in matrix notation

# i. Rank(C)

C
qr(C)$rank

# C does not suffer from linearity problems
C[1,]/C[2,]

# j. -G
G
-G

# k. Trace(E)

# Trace is the sum of diagnol elements

E
sum(diag(E))

# l. Rank(F)

# Because F is a diagnol matrix (with no zero diagnol elements) the rank must be equal to the min of the dimensions, 3.
qr(F)$rank

# kF, where k = -7

-7*F

# n. BF

B
F
B %*% F

# o. FB

F %*% B
# With matrices, BF != FB

# p. determinate of C or |C|
C
det(C)

# q. |D|
D
det(D)
# Which looks kind of funny but that is because R usings search algorithms to find the determinate and they are not exact.

# But calculating the determinate of a 2d matrix is easy:
D[1,1]*D[2,2]-D[1,2]*D[2,1]

# r. | CD |
C; D
H = C %*% D
det( H )
# Which also happens to be zero

H[1,1]*H[2,2]-H[1,2]*H[2,1]

# Sorry for the repetitive examples.  I had homework for a class and thought I might as well turn it into a post that someone might find useful.

Drawing jointly distributed non-normal random variables

* This method only approximates joint non-normal draws (which is really what any method does).

* I was recently told that it was "impossible" to draw joint non-normal distributions.

* But you will see that the approximation looks pretty good.

* It is easy to draw jointly distributed non-normal draws so long as you can start by drawing jointly distributed normal draws.

set more off

clear
set obs 10000

* For instance let's draw four variables.
* 1. a chi2 with 5 degrees of freedom
* 2. a poisson k = 5
* 3. a uniform variable with min = -5 and max = 5
* 4. a random f distribution draw with 5 and 5 degrees for numerator and denominator degrees of freedom.

* First we will specify the correlation matrix.
* The only constraint as far as I know is that the covariance matrix has to be PSD.
* This in practicality limits the possible correlations between variables since cross terms tend to cause vialations more likely in the PSD requirement.

matrix c = (  1, .7,-.3,  .2 \ ///
             .7,  1, .2, -.1 \ ///
    -.3, .2,  1,  .3 \ ///
     .2,-.1, .3,   1 )

* If we do not specify a mean or covariance matrix then the default draws are standard normals which is what we want for simplicity.
drawnorm x1 x2 x3 x4, corr(c)

corr x?
spearman x?

* Now all we need to do is turn our normal draws into uniform draws.
* Note: if x~N(0,1) and THETA is the CDF of the normal then y=CDF(x)~uniform
* So for any new distribution with CDF ALPHA and inverse INVALPHA the variable z=INVALPHA(y) ~ alpha.

gen y1 = normal(x1)
gen y2 = normal(x2)
gen y3 = normal(x3)
gen y4 = normal(x4)

sum
* Looking good.  The next step is that we take the inverse CDF of the distributions of interest.

gen z1 = invchi2(5, y1)
  label var z1 "chi2"
* The inverse poisson distribution seems to be incorrectly defined in Stata so that it uses 1-p rather than p to calucalate the inverse.
gen z2 = invpoisson(5, 1-y2)
  label var z2 "Poisson"
* It is easy to transform a uniform (0,1) to (a,b) by subtracting a and multiplying by (b-a)
gen z3 = y3*10-5
  label var z3 "Uniform"
gen z4 = invF(5, 5, y4)
  label var z4 "F distribution"
corr z?
spearman z?
* We can see that the spearman rank correlation is maintained with the standard pearson correlations are only slightly diminished by the non-linear transformations.

* In general the correlations are slightly drawn towards zero so if possible it might be worth it to exagerate the correlations in the matrix c so that they end up being drawn more closely to the desired levels.

hist z1, saving(chi2, replace) nodraw
hist z2, saving(poisson, replace) nodraw
hist z3, saving(normal, replace) nodraw
hist z4, saving(invF, replace) nodraw

graph combine chi2.gph poisson.gph normal.gph invF.gph

* Much of the content of this post was covered in a previous post under the title: Drawing Rank Correlated Random Variables.  It might be worth looking over the previous post if you have additional questions.

Hot Decking!

* Hot decking is a method commonly used in statistics to imput values where missing data is present.

* Let's see how it works!

* Imagine we have a data set of 200,000 people.

* Some of the questions were not answered by those people but most people did answer the majority of the questions.

set seed 1010

clear
set obs 200000

gen male = rbinomial(1,.51)

gen age = rpoisson(40)

gen education = rpoisson(12)

gen legality = rbinomial(1,.15)

gen parents_assets = rpoisson(4)

gen social_network = rbinomial(4, .1)
* Number of close friends

gen race = ceil(runiform()*4)
* Let's say there are 4 "races" of approximately equal representation

* Note that if we are hot decking across all of these characteristics then there are going to be a number of potential hot decks equivalent to the number of options from these different variables ie 2 for gender 2 for legality 4 for social networks and 4 for race so 2*2*4*4=64 without taking into account age which has a mean of 40 and variance of 40.  Likewise parents_assets with a mean of 4 and variance of 4 thus sd of 2.  Since the poisson distribution does not have an upper limit: age, education, and parents assets could present a problem but in practice there is a low probability of draws larger or smaller than 2 standard deviations as the poisson begins to look very much like a discrete normal distribution as k gets large.

* Now let's have some prediction variable

gen u = 2*rpoisson(40)

gen earnings = male + .01*age + .3*education + .1*social_network + legality + parents_assets + race + u

* The ideal case would be if we could to the OLS

reg earnings male age education social_network legality parents_assets race

/*


      Source |       SS       df       MS              Number of obs =  127112
-------------+------------------------------           F(  7,127104) =  772.07
       Model |  867810.062     7  123972.866           Prob > F      =  0.0000
    Residual |  20409480.7127104   160.57308           R-squared     =  0.0408
-------------+------------------------------           Adj R-squared =  0.0407
       Total |  21277290.8127111   167.39142           Root MSE      =  12.672

------------------------------------------------------------------------------
    earnings |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
        male |   .8483197   .0711003    11.93   0.000     .7089644     .987675
         age |   .0118882   .0056213     2.11   0.034     .0008705     .022906
   education |    .315772   .0102739    30.74   0.000     .2956353    .3359086
social_net~k |   .0225761   .0592116     0.38   0.703    -.0934775    .1386298
    legality |   1.055723   .0992171    10.64   0.000     .8612589    1.250187
parents_as~s |    1.01567   .0177606    57.19   0.000     .9808596     1.05048
        race |   .9607688    .031791    30.22   0.000     .8984591    1.023079
       _cons |   79.86437   .2833975   281.81   0.000     79.30892    80.41982
------------------------------------------------------------------------------

*/

* However, in actuallity some of our data is missing.

foreach v in male age education legality parents_assets race social_network {

* There is a 1/16 chance of the value being missing

gen miss = rbinomial(1,`=1/16')

replace `v' = . if miss==1

drop miss

}

* One way of handling this would be to do our estimation but by dropping the missing data.

reg earnings male age education social_network legality parents_assets race

* We can see that though our individual explanatory variables only represent 1/6 missing of 100,000 because different people have different values missing accross all of our explanatory variables we have a supstantial drop in the number of observations because few people are not missing answering at least one of the questions.

* So we, are going to try to impute our missing values to strengthen our estimation power.

* First let us install a user written command for hot decking in Stata (http://ideas.repec.org/c/boc/bocode/s366901.html):

* We will need to temporary change stata's missing values

sum
recode male age education legality parents_assets race social_network (.=-9999)

* This is a somewhat tricky bit of code that I believe is working correctly but I easily may be mistaken.
foreach v in male age education legality parents_assets race social_network {

  * I want to create a list of variables that is absent of the current looping variable
  local byvars = subinstr("male age education legality parents_assets race social_network", "`v'", "",.)

  * Create a initial group
  qui egen grp = group(`byvars')

  qui sum grp

  di "For variable `v' we have " r(max) " hot decks"

  * Now let's generate a variable that indicates how many potential values to choose from in each hotdeck
  qui gen missing = 1 if `v' == -9999

  * Count the number of missings in each group
  bysort grp: egen missing_count = sum(missing)

  * Count the number of items in each group
  bysort grp: gen all_count = _N

  * The hot deck size is the number of observations in the deck less then number of missing observations
  bysort grp: gen draw_deck = all_count - missing_count

  * Finally we need to figure out were to start counting each group from the greater distribution
  sort grp

  qui gen n = _n
  * This is the best trick I have right now for specifying where the group "starts" in the vertical distribution.

  * Ie the minimum of the _n values is the starting position on the vertical list for that group
  bysort grp: egen pos_min = min(n)

  * Specify which card to replace for each by adding a random draw from available hot decks a draw to the starting place of the group
  qui gen replace_card = floor(draw_deck*runiform()) + pos_min if `v' == -9999

  * Now we have to make sure our data is arranged properly
  sort grp `v'

  * Replace our `v' value with the card from the relevant hot deck
  qui replace `v' = `v'[replace_card] if `v'==-9999

  drop grp missing missing_count all_count draw_deck n pos_min replace_card

}

recode male age education legality parents_assets race social_network (-9999=.)

sum

reg earnings male age education social_network legality parents_assets race

/*


      Source |       SS       df       MS              Number of obs =  172180
-------------+------------------------------           F(  7,172172) =  918.86
       Model |  1033602.53     7  147657.505           Prob > F      =  0.0000
    Residual |    27667352172172  160.696002           R-squared     =  0.0360
-------------+------------------------------           Adj R-squared =  0.0360
       Total |  28700954.6172179  166.692538           Root MSE      =  12.677

------------------------------------------------------------------------------
    earnings |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
        male |   .9096129   .0611052    14.89   0.000      .789848    1.029378
         age |   .0133704   .0050211     2.66   0.008     .0035292    .0232116
   education |   .3026236   .0091769    32.98   0.000     .2846371    .3206102
social_net~k |   .0415584   .0530193     0.78   0.433    -.0623583    .1454752
    legality |     .97673   .0899224    10.86   0.000     .8004842    1.152976
parents_as~s |      .9697   .0158518    61.17   0.000     .9386308    1.000769
        race |   .9613812    .027394    35.09   0.000     .9076896    1.015073
       _cons |   80.19244   .2510821   319.39   0.000     79.70032    80.68456
------------------------------------------------------------------------------

*/

* First thing to notice is that we have increased the number of usable observations to 172K, which is 50% more than that of the first regression.

* If the code is working properly, the reason I believe we have not regained our 200k observations is that many hot decks are only populated by missing values.

* Thus the hot decking algorithm has nothing to work with.

* I have been told that hot decking is different from classical measurement error in that it does not lead to attenuation bias.

* Which is interesting but also problematic.

* Observe the t values and rejection rates in the second estimation compared with the first.

* Uniformly the t-values are getting larger despite the estimates not always getting better.

* This is because some of the regressors are imputed and therefore cannot be trusted in an identical fassion to that of standard exogenous regressors.

Robust Hausman Test Fail?

* Robust Hausman Test Fail?

* The Huasman test is a commonly used to indicate an ideal choice between fixed effect and random effect estiamtors (in a panel data context).

* In this post I will attempt to violate the underlying assumptions in the Hausman test to see how well the test performs under non-experimental situations.

* To execute this post I will use the robust form of the test purposed by Arellano (1993) {http://ideas.repec.org/a/eee/econom/v59y1993i1-2p87-97.html}.

clear
set obs 10000

gen id=_n

expand 5
* We have 5 years of data per id


bysort id: gen year = _n
 
* Exlpanatory variables are serially correlated accross years
gen x1 = abs(rnormal())+year
gen x2 = abs(rnormal())+year

gen u = rnormal()*5

* Let's create a set of variables that are the means of x1 and x2.
bysort id: egen x1_mean = mean(x1)
bysort id: egen x2_mean = mean(x2)

xtset id

gen y1 = x1 + x2 + u
xtreg y1 x1 x2 x1_mean x2_mean, cluster(id) re
test x1_mean x2_mean
* It is not a requirement that the explanatory variable be independent and failure of independence of draws does not cause problems for the hausman test.

* Let's see what happens when y is no longer a linear function of our explanatory variables

gen y2 = x1^.97 + x2^.98 + u
xtreg y2 x1 x2 x1_mean x2_mean, cluster(id) re
test x1_mean x2_mean

* Non-linearities do not seem to have an obvious and problematic effect on the Hausman test (though both FE and RE are now inconsistent generally).

* Perhaps if there is noise in the measurement of x1 and x2, the Hausman test will suffer.

gen x1_true = x1+rnormal()
gen x2_true = x2+rnormal()

bysort id: egen x1t_mean = mean(x1_true)
bysort id: egen x2t_mean = mean(x2_true)

gen y3 = x1_true + x2_true + u
xtreg y3 x1_true x2_true x1t_mean x2t_mean, cluster(id) re
test x1t_mean x2t_mean
* Interestingly the test fails very badly.  As far as I know, under measurement error in the explanatory variables, there is no reason to use a FE estimator above a RE estimator.

xtreg y3 x1_true x2_true , cluster(id) re
xtreg y3 x1_true x2_true , cluster(id) fe

* Finally we would like to know what would happen to the test if the error (u) is correlated inviduals?

sort id
gen pctile = _n/(_N+1)
gen u2 = normal(pctile)*5
sum u2

gen y4 = x1 + x2 + u2
xtreg y4 x1 x2 x1_mean x2_mean, cluster(id) re
test x1_mean x2_mean
* We can see that the Hausman test once again seems to be working.

* So, the take way? Hausman works well even when the model is slightly misspecified or when errors are serially correlated or when there exists measurement error in the explanatory variable.

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.

Three ways to Create Unique Identifiers

* Imagine that you have just recieved a data set and there are three variables household_ID, city_ID, state_ID which are characteristics of the observations.

* There is a the possibility that you have multiple observations per combination representing different responses from different surveys.

* You would like to create a unique identifier for each of the potential outcomes

* Let's start with 25,000 observations
clear
set obs 25000

* Household IDs go from 1 to 50
gen household_ID = ceil(runiform()*50)

* State IDs go from 1 to 50
gen state_ID = ceil(runiform()*50)

* City IDs go from 1 to 20
gen city_ID = ceil(runiform()*20)

di "Thus " 50*50*20 " potential outcomes"

* By chance there will be some mutiple outcomes and some outcomes that are never drawn.

* To count the multiple outcomes we can use the duplicates report command

duplicates report household_ID city_ID state_ID

* Let's imagine that we have a couple of variables of interest

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

* Now let's create our unique identifiers

* The easiest way is by using the group command

egen ID = group(household_ID city_ID state_ID)
* This command will just create a unique identifier of for each combination of household_ID city_ID state_ID in a manner consistent with sorting those variables.

sort household_ID city_ID state_ID
* We can easily see this in the numbering of ID.

sum ID
* I get a max ID of a little less than 20k.  This means that there are a little less than 20k unique household_ID city_ID state_ID combinations.

* We can manually do a similar command as the group function through the following set of nested loops.
local i = 0
gen ID2 = .

 qui levelsof household_ID, clean
 local levelsv = r(levels)
   qui levelsof city_ID, clean
   local levelsvv = r(levels)
    qui levelsof state_ID, clean
    local levelsvvv = r(levels)

foreach v of local levelsv  {
  foreach vv of local levelsvv  {
    foreach vvv of local levelsvvv  {
      qui replace ID2 = `i' if `vvv' == state_ID & `vv' == city_ID  & `v' == household_ID
 local i = `i'+1
}
  }
}
* This method is obviously much slower than the group function and does create uniform numbering when each combination is not at least once represented.

* An alternative way of creating a unique ID would be to create a text variable which is a combination of the other IDs.

* This method, though more complex (and using up more system resources) may be desireable in that it yeilds a unique identifier which is easily identifiable.

gen ID3 = "H" + string(household_ID,"%2.0f" ) + "C" + string(city_ID) + "S" + string(state_ID)

* ID3 is less easy to count unique identifiers with.

* Now the only problem is that some of our data has multiple observations.

* In order to resolve these multiple observations we will decide a rule on how to handle x1 and x2 that reduces our observations to 1 per identifier (thus making it finally unique)

* I am not a huge fan of the collapse command because it is a destructive command but it can save a few steps.
collapse (mean) mean_x1=x1 mean_x2=x2 (median) med_x1=x1 med_x2=x2, by(ID)

sum ID
di "Now max ID (" r(max) ") equals the number of observations (" r(N) ") indicating that we have achieved our goal of creating a unique identifier.

Fun with Macros

* Stata can use macros in a manner that can be very creative.

* One use could be in a manner similar to that of functions in R.
* If we would like to find the value of y at different levels of x where y is defined as y=x^2

local x = 1

* This immediately evaluates the value of y
local y = `x'^2

di "`y' = " `y'

* Thus the following display will be the same as the previous one despite x changing in value
local x = 2
di "`y' = " `y'

* If on the other hand, we specify y as the following:
local y = "\`x'^2"

local x = 1
di "`y' = " `y'
local x = 2
di "`y' = " `y'

* Then thw two different displays are different because y is waiting to evaluate the x until the display command.

* We can do this with multiple locals
local z = "100/(\`y')"

local x = 1
di "`z' = " `z'
local x = 2
di "`z' = " `z'

* Globals work in an identical fassion
global a = "(\`z')^(\`x') - \`y'"

local x = 1
di "$a = " $a

local x = 2
di "$a = " $a

global b = "(\$a)*cos(\$a)"

local x = 1
di "$b = " $b

local x = 2
di "$b = " $b

# R uses functions in a standard manner for programmers
y = function(x) x^2
z = function(x) 100/y(x)
a = function(x) z(x)^x - y(x)
b = function(x) a(x)*cos(a(x))

y(1); y(2)
z(1); z(2)
a(1); a(2)
b(1); b(2)
# Produces the same results as above