Avoiding Loops in R: An Example with Principal Minors

Yesterday, I found myself wanting to compute a large subset of the second order principal minors of a matrix (diagonal-preserving minors; the ones for which the rows and columns kept are the same). Don't judge me for wanting to do this, and bear with me, there is a lesson about R programming here (mostly for me, read the comments for testing).

If you don't know what a principal minor is or forgot what it was (because most of us don't spend much time with principal minors), here is a short refresher on how to get a second-order principal minor (other orders are possible and there's a lot of related theory, but let's try to focus on this one task for now).

Step 1: Deletion. Take any matrix. Delete all but two rows and all but two columns. For example, suppose the matrix is 5x5 and let's keep rows 1 and 3 and columns 2 and 4 (lots of other combinations are possible).

In the back of your mind, think how easy this is in R.

Step 2: Calculation. Once we delete the irrelevant parts of the big matrix to get our 2x2 submatrix, obtain the principal minor by taking the determinant of the resulting 2x2 matrix. For our example, we compute 3x2 - 3x1 = 3.

Both steps in this calculation are easy to implement in R. For the first step, merely use R's fantastic syntax for extracting rows and columns of matrices. If the matrix is called mat, we efficiently perform the first step by:

little_mat = mat[c(1,3), c(2,4)]

That's easy enough. With little_mat in hand, the second step is even easier.

pm = det(little_mat)

The real trick to this programming problem is picking out all of principal minors I want to compute (diagonal-preserving ones that keep the same-numbered rows and columns; little_mat won't do because c(1,3) is not c(2,4)), and doing it efficiently. I also don't want to rewrite my code as my matrix changes size.

So, what do we do? A really bad An idea that will eventually compute the right answer is to write this using for() loops.




Barring a coding error, this code should eventually work, but it will be slow (or fast depending on how you write it). If you want to run through this code many times (like I did yesterday), efficiency is at a premium and you should immediately think of faster techniques -- like using the apply() function (or preallocating your storage vector, which actually saves time).

Before moving to apply(), we need a quick way to enumerate all of the ways to keep two out of N rows. If we were counting, that's N choose 2, also known as a combination. Fortunately, this step has already been implemented in the gregmisc library with a function called combinations(). Conveniently, combinations(N, 2) will return a matrix with all of the ways to select two numbers less than or equal to N as its rows. There are other options, but this is the default, which is perfect for our application.

For example, here is some output for N=5.


With this way of computing the rows/ columns we want to keep, all we need to do is put our simple code for computing a second order principal minor into a function. Let's call it minor2(). Then, use apply() to apply this function to each row (MARGIN=1 is for rows) of the combinations matrix. Here's the complete set of code that will work if you have a matrix called mat for which you want the diagonal-second order principal minors.



Because we avoid doing loops, this method is much faster. For a big project, this added efficiency is worth it. More generally, this technique of using apply() to avoid repeatedly looping is a good lesson to apply to other programming problems.

Update: Thanks to Berend for going through the trouble of system.time() testing the code in this post (and spiceman for a good comment as well). Read on into the comments for more, but for a preview of what we learned: (1) The copying via c() is what slows the loop down, (2) apply is slightly slower than writing the loop yourself, (3) combn(), the combinations() alternative, is really slow, and (4) one should always stand behind system.time() methods before talking about the efficiency of a method.

I much appreciate the feedback. One of the reasons I have this blog is to learn more about R (and those comments were definitely informative).

Recreational R: Simulating a Card Trick

In this post, I simulate an interesting card trick, which was described by a friend of mine named Xan. Here's Xan's description of the card trick:

• I put a deck of cards down face up on the table. Meanwhile you think of a secret number between 1 and 10 -- for exposition let's say you pick 4.
• One by one, I discard cards from the top of the deck. When we get to the 4th card -- and let's call it your special card -- you look at the number. For exposition, let's say your special card happens to be a 7. Then 7 secretly becomes your new number. Note that I don't know your special card.
• I keep flipping cards, and 7 cards later, you have a new special card and number yet again. Note that I still don't know your special card.
• This process continues -- me flipping cards at a constant rate, you secretly updating your special card and counting up to the next one -- until I decide to stop.

But I don't just stop on any card. I stop on your current special card. Which I'm not supposed to know.

The following script file is a crude first approximation to simulating this card trick for the first 150 cards. This crude strategy is simple: (1) Simulate a sequence of 150 cards, (2) Pick a starting value -- might as well do all of them at once, and (3) Implement Xan's trick on all of them (that's what the for loop does).



Here's the result of running this script three times and plotting 10 secretnum sequences: one for each starting value. Each run starts with 10 distinct secret numbers (with a default chance of 0.1 of the trick "working"), and in my three runs, every starting value of the secret number eventually converges on the same sequence of secret numbers. More on the intuition for this here.



There are a couple of things wrong with this simulation that may have stacked the deck in favor of it working. I fixed these before writing this post, but I think it is useful to see the quick and dirty solution before adding frills and realism (which is what I do now).

1. I merely simulated my 150 cards using a uniform distribution on 1 to 13 (then forcing the face cards 11, 12 and 13 to be 10s). This is equivalent to sampling with replacement, which may have different properties than sampling without replacement.
2. There are only 52 cards in a deck. How would this work with only one shuffle, but passing through the deck (almost) three times until 150 cards? This is a different problem. Maybe you wouldn't expect much improvement after one pass through the deck (because you've seen that sequence of cards before).

A clever use of the rep() command can create a "fresh" deck of cards without suit distinction (which doesn't matter for the trick). Now, how to use R to shuffle a deck of cards? The sample() command, of course. With this technique in mind, I wrote the following function to compute the probability that the trick works (that is, the probability that two independent guesses end up at the same secret number) given a number of decks and a number of passes through the entire stack of cards (for a given shuffle).



Once you define this function, we can save the probability of the trick working for 1000 different shuffles of the cards. This script does it for one deck and one pass through the deck.



Then, we can take the average probability, which feels like the probability the trick works (unconditional on the particular shuffle or guess of a secret number). In the simulations I have run, this number turns out to be around 0.74 for one deck and one pass through the deck and about 0.86 using one deck and passing through twice. There's not much improvement after the second pass through of the deck.

In case you're interested, here are the results for up to 3 decks and up to 5 passes through the set of cards. It is always better to have more decks, but two passes through two decks is about enough to virtually guarantee that the trick works.



Update: In case you are wondering, I also coded it up with a 0.01 probability of forgetting to update on any given flip of the card. Here's the table of results in this case.



It's no surprise: For forgetful people, the trick is less successful and there isn't much point in reaching for a second deck.

An Application of boot() to IV regression

Bootstrapping standard errors can be a useful technique when obtaining closed form for the standard error formula is difficult or intractable. In this post, I give an example of how to use R to create a bootstrap sampling distribution in the context of IV regression. Specifically, I use boot() to automatically augment a function of mine to resample the indices of my data set with replacement (see the code in the function below).

In my application, I present a function that uses the boot library to report bootstrap standard errors for an instrumental variables regression. My ivboot() function builds on the iv() command I wrote for the tonymisc library.*

Now, onto the ivboot() example. Start by loading the tonymisc library and some data from the library, as well as the boot library.



We can run the following R script to define the function ivboot().



When applied to data and an IV regression model, the ivboot() function creates an IV regression object that -- when tonymisc is loaded -- is compatible with mtable() output. The only difference between an ivboot()-created object and an iv() object is that the ivboot() object has standard errors that are based on the bootstrap distribution of the coefficient estimates (Statistics other than the second-stage coefficient estimates are not bootstrapped).

Here is some code to compare the bootstrap output to the analytical standard error output:



On this small sample (N=100) of simulated data, the ivboot() command took less than a minute to run on my computer (timing of may vary depending on your computer). For much larger data sets, this will be slower. If you have a larger problem or lower standards (or higher standards, your choice), you can use the boots option to ivboot() to specify the number of bootstrap samples. Currently, I have set the default to 500, but you could specify boots = 200 if you want the command to run faster (boots = 10 will make it run even faster, but I don't recommend that!).

Here is the mtable() output, which can easily be ported into LaTeX using the toLatex() command.


*This standard output from an mtable() extension to my iv() command provides quite a bit of information in a convenient format. Another nice feature of iv() is that iv()-created objects have first-stage summary information readily stored in the object for extraction and analysis.

Defaults, Lists and Classes: A Functional Post

In this post, I demonstrate a couple of useful tricks to writing functions in R. The context is a function I wrote called called samp() that allows for an easy demonstration of sampling distribution properties.



Defaults

By default, this function draws K = 500 samples of size N from a normal distribution with a mean of mu = 0 and a standard deviation of std = 1. As you can see, the syntax for specifying a default for an argument to the function is to set it equal to its default value inside the function statement. Arguments can be given no default by naming them as an argument to the function without specifying their default value.

Sensibly chosen defaults can save a lot of time for people who use your function. For example, although this samp() function has four arguments, I don't need to type samp(100, 0, 1, 500) every time I wish to take 500 samples of size N=100 from a Normal(0,1) random variable. If sampling from a Normal(0,1) is especially common, setting the defaults this way can save a lot of keystrokes.

On the other hand, using defaults instead of hard coding the Normal(0,1) choice has the distinct advantage that users get flexibility from your function giving them the option to specify arguments other than the default. The right use of defaults strikes the perfect balance between ease of use and flexibility.

Returning a List... with Names

If your function warrants the special treatment, you have probably performed several useful calculations that you would like to return. Maybe you have conducted several related hypothesis tests that you would like to reference later (report and/or use in another calculation). A great way to do this is to store the results of your function into a list.

For example, in the samp() function, I wanted to store vectors of the sample means, sample standard deviations and standard errors (500, one for each replication). For good measure, I wrote the function to store the means and standard deviations of these vectors as well as the true parameter values that went into the function.

Although it is often a good idea to store the results in a list, it is a better idea to name the elements of the list for easy extraction. If I ran the command

mysamp = samp(100)

The object mysamp would contain a list of the results from the function samp(). Suppose I want to extract the vector of means (for example, to compute a histogram). Without running the names command, I would have to remember that the vector of means was in the first position in the list and use the command

hist(mysamp[1])

to plot a histogram of the means vector from my sampling object. As it is easy to forget the precise order in which the elements of a function are stored, this syntax can lead to too much thinking. It is better to name the first list element something like means. Naming the list elements allows the user to use the $ extractor on objects created with the function. That is, the syntax becomes

hist(mysamp$means)

which is much easier to remember and easier to read. This latter fact is often unappreciated, but if you get in the habit of naming elements in a list, you will be a better collaborative coder. Even if you're not into collaboration, when you return to that project after two weeks of doing something else, it is much easier to remember where you left off.

A Benefit of Keeping It Classy

Another coding practice that often goes unappreciated is the use of classes to make your life easier (and your console cleaner). At this point, if you type mysamp into your console window, R will bombard you with more than a screen's worth of output. Most of this is output that I don't want to see printed out. I might want to save it for later, but I don't want to save it.

This is where the class(result) = "samp" command comes in handy. By "classing" my object this way and writing a short a print method for this class (named print.samp()), I can cut down on the amount of output I have to see when I inspect my object. Here's code to define a print method.



Now, try typing mysamp into R. You will see much less output (only the output from print.samp), but the object mysamp still has all of the information you want to keep around. If you don't believe me, just type mysamp$means to be sure.

Just as we could "class and print," we can also "class and plot." This can be handy if you want there to be a standard meaning for the command plot(mysamp). In our working example, let's make this command mean plot three histograms, one for means, one for standard deviations and one for standard errors. Here's some code that imposes this meaning on objects with class "samp"



In case you were wondering what this plotting out looks like, here is the picture that is produced from the command plot(mysamp).