Export R Results Tables to Excel - Please don't kick me out of your club

This post is written as a result of finding the following exchange on one of the R mailing lists:

Is-there-a-way-to-export-regression-output-to-an-excel-spreadsheet

[Make sure to check out the many great comments on the bottom of the post.  That is where some better answers to this problem can be found.]

Question: Is there a way to export regression output to an excel spreadsheet?
Translation: I would like to be able to do a very simple thing that almost any statistical programming language can easily do, please suggest a basic command to do that.

Response1: ?lm ?coef ?write.csv ...
Translation: Read the manual and try this bit of incomplete code.

Questioner: I am very very new with R... Is there some simple code I could just paste?
Translation: Really? Isn't there anything you could suggest?

Response2: This is the help-you-learn-R mailing list, not the do-my-work-for-me mailing list...
Translation: Go F%$# yourself, freeloader. We only answer interesting questions.

Me:
I just wanted to say that it is just this type of response that gives R-users a bad reputation. I am an active R user and very happy to contribute to R in whatever way is possible but when I see posts like this, it makes me want to switch to a language in which the users are NICE people.  Okay, I know, I know.  A few bad apples should not spoil the basket, but sometimes things just taste rotten.

First thing, the reason I even stumbled across this post was because I had the same very similar question. Looking at the hits on the bottom of the page I can see that there are over 400 people who have viewed this discussion I am guessing most of them because they were looking for a specific solution rather than being interested in seeing how quickly experienced R users could could offend new users (for which there are numerous other examples).

In all likelihood, a lot of other new R users have come across this same post and been equally confounded but this rude and ridiculous response.

The original user who asked this question asked a very simple question for which any statistical language should have a very simple canned response.  Something along the lines:

lmOut(mylm, file="results", filetype="csv")

Yet the response that was instead produced was one which was overly complex, patronizing, and ultimately needlessly insulting.

I have written a little program to help format regression summary statistics into spreadsheet formats easily read by excel. Sorry if this is redundant. I am sure hundreds of people have programmed similar solutions. But I think it might be useful to many users who are not very familiar with how R constructs results.


Find the code on this gist:
https://gist.github.com/EconometricsBySimulation/6274532

Output looks something like this:

Looking over the original exchange it does not look like this code would even work for a logit for which it was originally needed.  However, I will post it on github and perhaps others will find the concept useful enough to make revisions (unlikely).  I might take another go at making it more general in the future though in all likelihood some user will send me an angry message saying "this has already been done by ...".

Easy Monte Carlo Sampler Command

R script

# In R I am often times unsure about the easiest way to run a flexible Monte Carlo simulation.
# Ideally, I would like to be able to feed a data generating function into a command with an estimator and get back out interesting values like mean estimates and the standard deviation of estimates.
# Stata has a very easly command to use called "Simulate".
# I am sure others have created these types of command previous but I don't know of them yet so I will write my own.

MC_easy = function(dgp, estimator, obs, reps, save_data_index=0) {
  # dgp (Data Generating Proccess) is the name of a function that creates the data that will be passed to the estimator
  # estimator is the name of a function that will process the data and return a vector of results to be analyzed
  # obs is the number of observations to be created in each spawned data set
  # reps is the number of times to loop through the dgp and estimation routine
  # save data index is a vector of data repetitions that you would like saved

  # I typically start looping programming by manually going through the first loop value before constructing the loop for both debugging and demonstrative purposes.

  # This command will create a a starting sample data set.
  start_data = get(dgp)(obs)
  # This command runs the estimate values on that starting data set.
  MC_results = get(estimator)(start_data)

  # Create an empty list of save data values
  save_data = list()
  # If the first data set generated has been specified as a data set to save then save it as the first entry in the save_data list.
  if (sum(1==save_data_index)>0) save_data[[paste("1",1,sep="")]] = start_data

  for (i in 2:(reps-1)) {
    temp_data = get(dgp)(obs)
 
    # If this repetition is in the save_data_index then save this data to the list save_data
    if (sum(i==save_data_index)>0) save_data[[paste("i",i,sep="")]] = temp_data
 
    MC_results=rbind(MC_results, get(estimator)(temp_data))
  }
  # Display the results of the estimations
  MC_results = as.data.frame(MC_results)
  print(rbind(mean=mean(MC_results), sd= sd(MC_results), var= sd(MC_results)^2))

  # If the number of items for which the data was saved is equal to zero then only return the results of the estimation.
  if (sum(save_data_index>0)==0) return(MC_results)

  # If the number of items is greater than zero also return the saved data.
  if (sum(save_data_index>0)>0)  return(list(MC_results,save_data))
}

#### Done with the Monte Carlo easy simulation command.

# Let's try see it in action.


# Let's define some sample data generating program
endog_OLS_data = function(nobs) {
  # Let's imagine that we are wondering if OLS is unbiased when we have a causal relationship between an unobserved variable z and the observed variables x2 and x3 and the unobserved error e.
  # x4 is an observed variable that has no relationship with any of the other variables.
  z = rnorm(nobs)
  x1 = rnorm(nobs)
  x2 = rnorm(nobs) + z
  x3 = rnorm(nobs) - z
  x4 = rnorm(nobs)

  e = rnorm(nobs)*2 + z

  y = 3 + x1 + x2 + x3 + 0*x4 + e
  data.frame(y, x1, x2, x3, x4)
}

sample_data = endog_OLS_data(10)
sample_data
# Everything appears to be working well

# Now let's define a sample estimation
OLS_est = function(temp_data) {
  # Summary grabs important values from the lm command and coef grabs the coefficients from the summary command.
  lm_coef = coef(summary(lm(y~x1+x2+x3+x4, data=temp_data)))

  # I want to reduce lm_coef to a single vector with values for each b, each se, and each t
  lm_b = c(lm_coef[,1])
    names(lm_b) = paste("b",0:(length(lm_b)-1),sep="")
  lm_se = c(lm_coef[,2])
    names(lm_se) = paste("se",0:(length(lm_se)-1),sep="")
  lm_se2 = c(lm_coef[,2])^2
    names(lm_se2) = paste("se2_",0:(length(lm_se2)-1),sep="")

  lm_rej = c(lm_coef[,4])<.1
    names(lm_rej) = paste("rej",0:(length(lm_rej)-1),sep="")
  c(lm_b,lm_se, lm_se2,lm_rej)
}
OLS_est(sample_data)

# From this is it not obvious which estimates if any are biased.
# Of course our sample size is only 10 so this is no surprise.
# However, even at such a small sample size we can identify bias if we repeat the estimation many times.

# Let's try our MC_easy command.

MC_est = MC_easy(dgp="endog_OLS_data", estimator="OLS_est", obs=10, reps=500, save_data_index=1:5)
# We can see that our mean estimate on the coefficient is close to 3 which is correct.
# While b1 is close to 1 (unbiased hopefully) while b2 is too large and b3 too small and b4 is close to zero which is the true value.
# What is not possible to observe with 500 repetitions but capable of being observed with 5000 is that the standard error is a downward biased estimate of the standard deviation.
# This is a well known fact.  I think even there is a proof that there is no unbaised estimator for the standard error for the OLS estimator.
# However, se^2 is an unbiased estimator of the variance of the coefficients.
# Though the variance of se^2 is large, making it neccessary to have a large number of repetitions before this becomes apparent.
# The final term of interest is the rejection rates.  We should expect that for the coefficients b0-b3 that the rejections rates should be above random chance 10% (since there trully exists an effect) which is what they all ended up being though with b1 and b3 the rejection rate is only slightly above 20% which indicates that our power given our sample size given the effect size of x1 and x3 is pretty small.
# The final note is that though we do not have much power we are at least not overpowered when it comes to overrejecting the null when in fact the null is true.  For b4 the mean rejection rate is close to 10% which is the correct power.

# We can recreate the results table by accessing elements of the MC_est list returned from the function.
rbind(mean=mean(MC_est[[1]]),sd=sd(MC_est[[1]]),var=sd(MC_est[[1]])^2)

# We may also be curious about the medians and maxes.
summary(MC_est[[1]])

# Or we may like to check if our data was generated correctly
MC_est[[2]]

ANOVAs and MANOVAs

# Anova's are frequently used in experimental setting when treatments is a categorical value and there are one or more response variables.  Rather than testing the statistical significance of each categorical value with respect to the response variable.  We instead test the joint significance.

# Let's see how this works.

# Imagine we are a pharmaceutical company and would like see is the effect on sleep of three different drug combinations at three different levels each.

# In order to test each level for all of the drug combinations we would have 27 combinations.  If we tested the individual significance of each combination then without adjusting the alpha value we would over-reject far too many times.

#
asbermatimitite=rbinom(1000,3,.5)*5
zugrimiitosoite=rbinom(1000,3,.5)*5
crelitotiserite=rbinom(1000,3,.5)*5

# Let's say each of our subjects on average sleeps a number of hours as a random draw from the poisson distribution.
base.sleep = rpois(1000,8)

sleep.hrs = base.sleep+asbermatimitite/5-zugrimiitosoite/2.5+crealitotiserite*0

# Let's say the most frequent side effect of these drugs is drowsiness.  We would also like to know the effect of the drugs on drowsiness levels.

base.drowsiness = runif(1000)

drowsiness = base.drowsiness + asbermatimitite*.01 + zugrimiitosoite*.02 - crealitotiserite*.01

# From this setup we can see that crelitotiserite does not actually increase sleep but does reduce drowsiness.

boxplot(sleep.hrs ~ asbermatimitite+zugrimiitosoite+crelitotiserite, horizontal = T, main="# of hours slept", ylab="Asbermatimitite.ZugrimiitosoiteCcrelitotiserite (Mg)")

# From the box plot it is clear there is movements in the means as the explanatory variables changes.  Other that that general statement I am no sure how else to read the grpah.

boxplot(drowsiness ~ asbermatimitite+zugrimiitosoite+crelitotiserite, horizontal = T, main="Self reported drowsiness", ylab="Asbermatimitite.Zugrimiitosoite.Crelitotiserite (Mg)")

# Likewise with drowsiness levels.  There is movement but it is hard to tell where and why.

man <- manova(cbind(sleep.hrs,drowsiness)~asbermatimitite+zugrimiitosoite+crelitotiserite)

# This will save the results of the manova in the list called man.

summary(man)

# We can see that all of the explanatory variables are significant.  crelitotiserite less so than the other two.  Unlike regression analysis we do not know what to do with this information now.  All we can say is that given random treatment there are differences in one or more of the dependent variables as a result of the different levels of treatment   I next logical approach would be to attempt to estimate the form of the effects.  If one of the variables was not significant then it might have been reasonably dropped from the analysis.

Stata is to Accounting as R is to Tetris

Both Stata and R handle many of the same data computation needs.

However, researchers must subset data within them very differently.

For Stata subsetting can be extremely easy.

If you want to restrict your data you can simply post after most commands an if statement.
 
  * Stata code
  clear
  set obs 100
  gen y = rnormal()
  gen x1 = rnormal()
  gen x2 = rnormal()
  gen u = rbinomial(1,.5)
  reg y x1 x2 if s==1

Thus the OLS regression of y on x1 and x2 will only occure if s=1.  In R this operation can be a little more tricky. Imagine you have a data set called mydata which has four variables y x1 x2 s.  The easiest way to restrict the data would probably be to create a new subset data set.

  # R code
  # Create your data set
  mydata = data.frame(y=rnorm(100), x1 =rnorm(100), x2 =rnorm(100), s = rbinom(100,1,.5))
  # Create a sub-set of your data by specifying the subset mydata[mydata$s==1,]
  lm(y~x1+x2, data=mydata[mydata$s==1,])

Thus same operation as above.  

Those of you with less experience in R are probably wondering how using brackets accomplished the same task.

In R like in Stata you use brackets to indicate subcripts.

For instance:  The vector "letters" is a built in environmental vector in R that contains all of the letters from a to z.

Thus:

  letters[1] # Displays "a"
  # Vectors can also be subscripted by vectors (with repetition)
  v = c(1,2,3,2,1,5)
  letters[v]
  # Displays "a" "b" "c" "b" "a" "e"
 
  # Vectors can also be subsetting using logical operators.
  vv = rep(c(TRUE,FALSE),13)
  # Creates a vector 26 elements long alternating between TRUE and FALSE
  letters[vv]
  # Will display every other letter starting with a.
 
  # This brings us back to how we subsetted a dataframe.
 
  # Let's make a new data frame called mysample
  mysample = data.frame(a = letters, b = 1:26, c = rnorm(26))
 
  # We can subset the data frame by using two subscripts now
  mysample[4,2] # Displays 4
  mysample[3,]  # Displays an entire row of the data frame.
  mysample[vv,] # Will display every other row
  # Replacing subsets is notated in a similar manner as subsetting
  mysample[vv,1] <- "z" # Will replace every other letter with "z"
 
  # Thus mydata[mydata$s==1,] is telling R to use any row in which variable s of data frame mydata is equal to 1.
 
At this point you are probably thinking that R is overly complicated and that Stata handles data much better.  This is not true.

Stata handles data in a manner similar to that of an accountant.  If you want your accountant to add within rows different values there is no problem.  You can even use subscript to move values from one row to another.  R on the other hand takes data and transforms it and combines it into new forms often much easier than Stata but with more complex notation.

Matrix Operations in Mata

* This post demonstrates a few methods for how to input matrices into Mata and how to do some basic matrix operations.  Mata is a matrix programming language so basic matrix operations are extremely easy.

mata
// Let's first build some matrices in Mata
// They can be built directly
A = ( 2 , -1, 5 \ 3, 0 , -1 \ 3, 3 , -1)
A

B1 = (0 , 3 ,-1)
B2 = (3 , 0 , 0)
B3 = (0 , 2 , 0)

// Or through a combination of vectors
B = (B1 \ B2 \ B3)

B

// We can also start with an empty matrix and fill in values

C = J(3,2,4)
// The J command creates a matrix with 3 rows, 2 colums, and with default values of 4

C
// We can replace individual elements once the matrix is created
C[1,2] = 3
C[3,1] = 7

// Or entire submatrices
C[2,] = (2,6)
C

// Now let's see how various matrix operations perform in Mata:

// a. AB
A*B

// b. BA
B*A

// c. A+B
A+B

// d. A'B'
A'*B'

// e. B'A'
B'*A'

// f. A'B
A'*B

// g. (AB)'
(A*B)'

// h. ABC
A*B*C

// i. C'AC
C'*A*C

// j. CAC'
C*A*C'
// This does not exist because C*A cannot be multiplied as C is 3x2 and A is 3x3

// k. trace(C'AC)
trace(C'*A*C)

// l. trace(CAC')
trace(C*A*C')
// Likewise the trace does not exist

end

* Matrices can also be input into Mata from data sets.

clear
set obs 5
gen y = 3
replace y = 6  if _n == 2
replace y = 10 if _n == 3
replace y = 8  if _n == 4
replace y = 2  if _n == 5

mata
// The command st_data retrieves data from stata variables to be used in Mata
y = st_data(. , "y")
y

// Once input, matrices can easily be manipulated
y*y'

// The square of the norm of y
norm(y)^2

// This happens to be identical to:
sum(y:^2)

// The norm is the equclidean norm which is the square square root of the sum of all of the squares of a vector.
// Thus the square of it is just the sum of the squares.

// As should be clear, manipulating matrices in Mata is extremely easy.
// Thus Stata is able to pack a powerful Matrix programming language inside an effective high level user language.

// The largest frustration that I have had with Mata is the relative quality of the documentation.
// I find the documentation of Mata much harder to use than that of Stata (at least in version 11, perhaps version 12 has better documentation).
end

Mata speed gains over Stata

* The inclusion of Mata as an available alternative programming language for Stata users was a great move by Stata.

* Mata in general runs much quicker than programming on the surface level in Stata.

* In Stata each loop that runs is compiled (interpretted into machine code) as it runs creating a lot of work for the machine.

* In Mata on the other hand, the entire loop is compiled prior to running.

* Let's see how this works.

* Let's say we want to add up the square of the numbers 1 through 100000

* Method 1: Surface loop

timer clear 1
timer on 1
local x2 = 0

forv i = 1/1000000 {
  local x2 = `x2'+`i'^2
}

di `x2'

timer off 1
timer list 1

* On my laptop, this takes about 13.5 seconds

* Method 2: Mata loop
timer clear 1
timer on 1
mata
  x2=0
  // This command can be read as start i at 1,
  // keep looping so long as i is less than 1000000,
  // the third argument looks a little fishy but it is syntax
  // that has been around for a while (at least since C).
  // It would be identical to writing i=i+1, in other words, add 1 to i.
  // Following the for loop we can immediately place a since line command.
  for (i = 1; i <= 1000000; i++) x2=x2+i^2
  // If there is nothing done with the value x2 then mata displays this value.
  // R handles this identically
  x2
end

timer off 1
timer list 1

* In contrast, my computer completed the loop using mata in .27 seconds, many magnitudes of speed faster.

* However this does not mean you need to learn to use mata (since it has its own limitations and syntax) in order to speed up your commands.

* Method 3: Use Stata's data structure to accomplish vector tasks
timer clear 1
timer on 1

clear
set obs 1000000
gen x2 = _n^2

* The sum command will calculate the mean of x2 which is the same as the sum of x2 divided by it's number of observations.
sum x2
* We can reverse that operation easily.
di r(N)*r(mean)

timer off 1
timer list 1
* Using a little knowledge of how Stata stores post command information this method does the same trick in .2 seconds

* Method 4: The speed gains in 3 was as a result of using the vector structure of data columns.  Mata can do very similar things even easier.

timer clear 1
timer on 1
// This command looks a little fishy, but it is easy to understand.
// Order of operations must be taken into account.
// First the 10^6 is evaluated which equals 1000000
// Then the vector 1..10^6 is made which looks like 1 2 3 ... 1000000
// The .. tells mata to make a count vector.
// If I had written :: then mata would have made a column vector instead.
// Once the vector is made then the command :^2 tells stata to do a piece wise squaring of each term in the vector.
// Finally the sum command adds all of the elements of the vector together to generate the result we were looking for.
mata: sum((1..10^6):^2)
timer off 1
timer list 1
* The result is that this command only took .04 seconds to run through efficient coding in Mata.

# As a matter of comparison, this command
system.time(sum((1:10^6)^2))
# took .04 seconds in R

# And the loop:
x=0
system.time(for(i in 1:10^6) x=x+i^2)
# 1.3 seconds

# Thus Mata in this example is significantly faster than Stata and about the same speed as R.

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.