Two things worth knowing before we start, because they will save you time all semester. The three stacked lines in the bottom left corner open a table of contents, so you can jump straight to a section rather than arrowing through every slide. And pressing the letter o gives you a panel view of the whole deck at once, which is the fastest way to find the slide you half remember. Both matter more in this chapter than most, because we move back and forth between the loop material and the parallel material, and you will want to compare them.
Click on the three horizontally stacked lines at the bottom left corner of the slide, then you will see table of contents, and you can jump to the section you want
Hit letter “o” on your keyboard and you will have a panel view of all the slides
The code boxes on these slides are live. The faintly blue area is an editor you can type in, and Run Code executes everything in it. If you only want part of it, highlight that part and press Command Enter on a Mac, or Control Enter on Windows. Two buttons in the top right corner matter: the stacked-paper icon copies the code so you can paste it into RStudio on your own machine, and the reload button puts the original back if you have edited yourself into a corner. Do use these. This is a chapter where you learn far more by running things than by reading them. One more, at the right end of the toolbar: the eye icon hides that cell’s output and a second click brings it back. Some of these results run long, and once you have read one it is just pushing the rest of the slide out of view.
Four things by the end of this chapter. First, writing your own functions, which is the foundation for everything else here. Second, repeating work with a for loop and with lapply. Third, and this one surprises people, learning when not to loop at all, because a great deal of what looks like a loop problem is really a vectorization problem and runs far faster written that way. Fourth, parallelizing genuinely slow work with future underscore lapply. Notice the order. Parallelization comes last deliberately. It is the most exciting tool and the least often the right one, and you should reach for vectorization first almost every time.
lapply() to complete repetitive jobsfuture_lapply() function from the future.apply packageThe honest answer to when you should write a function is: when you notice yourself about to do the same thing a second time with different inputs. The mean example on screen makes the point cheaply. You could write sum of x divided by length of x every single time you need an average, and it would work, but nobody does that, because someone wrote mean for you. Your own functions are the same idea at a smaller scale. And the payoff grows with the length of the task. For a one-line calculation a function barely pays for itself. For a fifteen-line procedure you run on forty datasets, it is the difference between maintainable and unmanageable.
It is beneficial to write your own function when you expect to repeat the same action with different inputs to the action.
Example: mean()
Calculating the average of a variable is such a common task
sum(x)/length(x) every time you get a meanmean() functionA function is more useful when the task is longer and more complicated.
Here is the skeleton, and it really is this simple. You give the function a name, you list its arguments inside function, and then in the body you do something with those arguments and return a result. Two things to notice. The return call is how the value gets out. Without it you may still get something back, but being explicit makes your intent obvious to whoever reads it next, including you in six months. And the note at the bottom matters more than it looks: the argument does not have to be called x. That is just a placeholder. Name your arguments after what they actually are.
Read the first line as an assignment. The left arrow stores the function object under the name function name, function of x declares one input, and the curly braces mark the body that runs each time you call it. The numbered lines inside are placeholders for the calculation and the returned result, not literal R code, which is why this demonstration chunk is marked not to evaluate. In the second block, function name of x is the call: R matches the value you supply to the argument x, runs the body, and gives the returned value back to the calling code.
Here is the general structure of a user-defined function:
Then, you can use the function like this:
Note: the argument does not have to be named x
Now the skeleton with something in it. This function is called square it, it takes one argument named a, it squares a and stores that in squared a, and then returns squared a. Map it against the structure on the previous tab and you will see the two numbered comments are exactly the two steps. Then run the second block, which calls square it with ten and stores the result in ten squared. The thing to notice is that you never told R what a is anywhere outside the function. The ten you passed in becomes a for the duration of the call, and then it is gone.
The caret two is the squaring operation, so the call produces one hundred. The first live block is set to autorun, which defines square it as soon as the slide loads and makes it available to the second block. The second block deliberately separates assignment from display: the left arrow saves the returned value as ten squared, and writing that object name on the final line asks R to show it. That is the normal pattern when you want to keep a function’s result for later work instead of merely printing it once.
The following function takes numeric numbers, square them, and return the squared values:
Try the function:
This tab is about where things live. Inside this function we create two objects, squared a and original a, and we return original a. Run it, then look at your environment tab in RStudio. Neither of them is there. Objects created inside a function exist only while the function is running, and then they vanish. That is a feature, not a limitation. It means a function cannot accidentally clobber a variable you are using elsewhere, and you can reuse obvious names like x or result inside functions without worrying. And read the last line carefully: we return original a, but only its contents come back, not the name.
Trace the shown call with eight. Inside the function, a caret two makes squared a equal sixty-four, and square root of squared a makes original a equal eight for this input. The first live block autoruns so the function exists before you try the call below it. The call displays the returned eight, but it does not export either internal name to the global environment. This is why checking the environment after the call is part of the demonstration: a visible result is not evidence that the function’s local working objects were created globally.
Any objects that are defined within a function are not registered (created) to the global environment.
For example, squared_a and original_a are defined only internally, and are not registered on the global environment.
You can confirm this by looking at the environment tab in RStudio after running the following:
Even though we are returning original_a, only its content is returned.
Now the direction that catches people out. Look at the function. It uses something called multiplier, but multiplier is not one of its arguments, and it is never created inside the function either. So run it, and it fails, because R has nowhere to find that name. Then define multiplier in the global environment and run exactly the same function again, and now it works. That is R’s scoping rule: when a name is not found locally, R goes looking outward. This is convenient and it is also dangerous, because a function that silently depends on a global variable will work on your machine and break on someone else’s. Pass things in as arguments.
The first call still computes ten squared as one hundred, but it stops at multiplier times squared a with an object-not-found error because multiplier is missing. The next block creates a global multiplier of ten and immediately repeats the call, so z becomes ten times one hundred, or one thousand, and that is the value returned. The contrast between those two calls is the evidence for the scope rule. If multiplier were an explicit argument, the function’s requirements would be visible at the call site and the result would no longer depend on whatever happens to be in your workspace.
When R sees objects that are not provided explicitly as arguments to the function, then R looks for them in the global environment:
Here, multiplier is not provided as an argument to the function, and it is not defined inside it either.
Try this:
Now, define multiplier on the global environment, and then run it again:
A small addition with a large payoff. Writing a equals five in the argument list gives that argument a default, so calling square it with no arguments at all still works and uses five. Defaults are how you make a function easy to call in the common case and still flexible in the unusual one. The pattern you will see in real packages is a function with two or three arguments you almost always supply and a dozen you almost never touch, all defaulted to something sensible. Your own functions should do the same. Default the things that usually stay the same, and require the things that always change.
In the live definition, autorun makes this revised version of square it replace the earlier one when the slide loads. The empty parentheses in square it of nothing mean that you supply no value, so R substitutes the default five and returns five squared, which is twenty-five. If you called square it with ten, the supplied ten would override the default. A default is therefore a fallback, not a fixed value.
You can set default values for function arguments by argument = value like below:
Try this:
Adding arguments is just adding them, separated by commas. This function takes a and b and returns a squared plus b squared. Then look at the two calls, because the difference between them is the lesson. In the first, four and three are matched to a and b by position: first value to first argument. In the second we name them explicitly, b equals three and a equals four, and now the order does not matter at all. Positional is fine for one or two obvious arguments. The moment there are several, or the moment you are looping, name them. You will see exactly why in the exercises.
Both arguments also have defaults, five for a and two for b, but the first call supplies four and three, so both defaults are overridden and the result is four squared plus three squared, or twenty-five. Autorun defines the function and performs that first call on load. The second live block then proves that named matching takes precedence over position: although b appears first in the call, R still sends three to b and four to a, so the answer remains twenty-five.
It is easy to create a function with multiple arguments. You just simply add more arguments within function() like below:
As you are likely to have noticed, the order of input arguments is assumed to be the same as the order of the arguments of the function. Above,
a = 4b = 3You can mess with the order of input arguments if you want if you name the input arguments as follows:
Three exercises, and they build on each other deliberately. The first is a straight translation: you are given the Fahrenheit to Celsius formula, so wrap it in a function. The second asks for a yield response function, taking a nitrogen rate and returning predicted corn yield from that log relationship. The third is the interesting one, because it asks you to calculate revenue using the function you wrote in exercise two rather than rewriting the yield formula. That is the habit worth building. Functions are meant to be built out of other functions, and if you find yourself copying the yield equation into a second place, you have missed the point.
Use the inner tabs to check each answer after you have worked in the empty live box. In Exercise 1, f to c subtracts thirty-two, multiplies by five, divides by nine, and returns the Celsius value. In Exercise 2, n to yield takes N, applies R’s natural logarithm to N, adds the baseline one hundred twenty, and returns yield. In Exercise 3, calc rev takes both corn price and N, calls n to yield with N for the intermediate yield, multiplies that yield by price, and returns revenue. The answer chunks are folded and marked not to evaluate, so opening an answer reveals reference code without silently defining those functions for you. That keeps the exercise honest while still letting you compare your structure line by line.
Define a function that takes temperature in Fahrenheit, convert it into Celsius, and return it.
Here is the formula: temp_C <- (temp_F - 32) * 5 / 9
Work here
Answer
After running a randomized nitrogen trial, you found the following relationship between corn yield (bu/acre) and nitrogen rate (lb/acre):
\text{corn yield} = 120 + 25 \times log(\text{nitrogen rate})
Write a function that takes a nitrogen rate as an argument, calculate the estimated yield for the nitrogen rate, and then return the estimated yield.
Work here
Answer
You would like to calculate the revenue of corn production (defined below) as a function of nitrogen rate based on the yield response function.
\begin{align*} P_{corn} * \text{corn yield} \end{align*}Write a function that takes corn price and nitrogen rate as its arguments, calculate revenue (and yield as an intermediate step), and return revenue. In doing so, use the function you created in Exercise 2.
Work here
Answer
Here is the case for loops. Squaring the numbers one through five by hand is five lines, and that is perfectly fine. Now imagine a thousand integers. Nobody is writing a thousand lines, and if they did, nobody could read it or change it afterwards. But the real argument is the last bullet. The repetitive work you will actually meet is not squaring numbers. It is things like Monte Carlo simulation, where you run a genuinely complicated procedure thousands of times with slightly different inputs. Writing that out by hand is not merely tedious, it is impossible. Loops are how you describe the pattern once and let R do the repetition.
We sometimes need to run the same process over and over again often with slight changes in parameters.
In such a case, it is very time-consuming and messy to write all of the steps one by one.
For example, suppose you are interested in knowing the square of 1 through 5 with a step of 1 ([1,2,3,4,5]). The following code certainly works:
However, imagine you have to do this for 1000 integers.
Yes, you don’t want to write each one of them one by one as that would occupy 1000 lines of your code, and it would be time-consuming.
Things become even more challenging when you need to repeat far more complex processes, such as Monte Carlo simulations. To handle such repetitive tasks efficiently, it is important to learn how to write programs that automate these jobs using loops.
Strip a loop down and it has exactly two parts: the thing you do, which stays the same every time, and the thing that changes. In the squaring example the action is squaring, and the changing part is which number gets squared. Getting into the habit of naming those two parts separately, before you write any code, makes loops much easier to write. Then the syntax underneath is just that idea spelled out. For x in some list of values, do something with x. R runs the body once for each value, with x taking the next value each time round.
Looping means repeatedly executing the same process, with only the parameters changing.
In the example above:
Syntax
Two loops that do exactly the same thing. The first uses x as the loop variable, the second uses bluh bluh bluh, and the output is identical. That is the point of showing you the silly one. The loop variable is just a name you choose, with no meaning to R at all. Beginners often think there is something special about i or x. There is not. Pick a name that says what the value is, like year or county or nitrogen rate, and your loop becomes self-documenting. Also notice that both bodies call print. That is not decoration, and the exercise tab shows you why.
In each header, one colon five supplies the integers one through five in order. On each trip through the braces, the current value is squared by the caret two and passed to print, so you see one, four, nine, sixteen, and twenty-five on separate iterations. Renaming the loop variable requires changing both its name in the header and the matching name in the body. It does not change the sequence, the operation, or the output. Move to the Exercise tab to build the sequence with seq instead of the colon operator.
This does the same:
Your turn. Cube each element of a sequence that starts at five, goes up in fives, and stops at fifty. So you need seq to build the values and a for loop to walk through them. One thing to watch, and the callout under the answer spells it out. Inside a for loop, R does not automatically print. If you write i cubed on its own line, R computes it, throws it away, and shows you nothing at all. You have to wrap it in print to see the numbers. That trips up almost everyone the first time, because at the console the same expression would have printed quite happily.
In the answer, seq from five to fifty by five supplies the ten values, i receives them one at a time, and the caret three cubes the current value. The answer chunk is folded and marked not to evaluate, so opening it reveals the reference loop without executing it for you. Compare your live result with the answer’s printed sequence, and then move to lapply to see how returning a collection changes this pattern.
Write a for loop that cubes each element of a sequence starting at 5, increasing by 5, and ending at 50.
Work here
Answer
Why print()?
Inside a for loop, R does not auto-print. Writing i^3 on its own computes the value and discards it silently. If you want to see it, you have to ask.
lapply()Base R gives you a second way to loop, and it is the one you will end up using most. The syntax is lapply of A and B, where A is the values you want to loop over and B is what you want done to each one. B can be a function that already exists, like mean, or one you write on the spot. The note at the bottom repays reading twice. What lapply loops over depends on what you hand it: elements of a vector, elements of a list, and for a data frame, the columns. That last one surprises people until you remember a data frame is a list of equal-length columns.
Instead of using a for loop, we can use the lapply() function from the base package to loop.
Syntax
A is the list of valuesB is the function you would like to apply to each of the values in A
mean())Note:
A is a vector, lapply() works on each of the vector elementsA is a list, lapply() works on each of the list elements whatever they may beA is a data.frame, lapply() works on each of the columns (data.frame is a list of columns of equal length)Same job as the for loop, in one line. Read the backslash of x as function of x. It is just shorthand, and you will see it everywhere in modern R code. But the important sentence is the one in red at the bottom. lapply always returns a list. Always. That is what the l stands for. So even though you fed it a simple numeric vector and did simple arithmetic, what comes back is a list of five one-element pieces, not a vector. That is the single biggest difference from a for loop, and it is why you will so often see unlist or bind rows immediately after an lapply.
Read the whole call from left to right: one colon five is the input, the anonymous function receives one element at a time as x, and x caret two is the final expression in its body. R implicitly returns that final expression, so this short function does not need an explicit return call. Unlike the earlier for loop, you also do not need print, because lapply collects every returned value for you. The displayed list is therefore the reusable result of the loop, not merely a sequence of messages sent to the console.
This does the same thing as the for loop example we looked at earlier:
\(x): shorthand for function(x)The key difference from a for loop is the object class of the output after the loop.
Important: the output type of lappy() is always a list (that’s why it is called lapply)
A small habit that pays off constantly. Rather than writing the action inline, define it as a named function first, then loop over that. Here we define square it, and then lapply calls it. For a one-line action this looks like extra work, and honestly it is. But the moment the action is fifteen lines long, having it as a named function means you can test it on a single value, get it right, and only then put it in a loop. Debugging a broken function is easy. Debugging a broken function that is also inside a loop over a thousand values is not.
The first block autoruns so square it is ready when you reach the second block. In the lapply call, one through five is still the input sequence, and the shorthand function passes its current x into square it. The named function squares that value and returns it, then lapply stores each return value in its output list. You could test square it with three before looping and know that any remaining problem must be in the loop wiring rather than in the squaring rule itself.
It is often the case that you want to write a function of the action you intend to repeat first and then loop.
For example, for the loop of squaring numbers, you can first define a function that implements the action of squaring:
And then loop:
A very common situation: your function takes several arguments, but you only want to loop over one of them. Here square them add takes a and b, and we want to try many values of a while holding b at five. The lapply line shows the pattern. The value being looped is called x inside the shorthand function, we pass it as the first argument, and we pin b equals five explicitly. Notice that b is named and the looped value is not. That is the safe way round, and it connects directly to what we said about naming arguments. Get the pinning wrong and your loop values quietly land in the wrong slot.
The autorun block first defines the two-argument function, where a caret two plus b caret two is stored and returned. The next block supplies one through ten to the anonymous function, so those ten values become a in turn while b remains five for every call. The results therefore run from one squared plus twenty-five through ten squared plus twenty-five. As always, lapply collects those ten values as a list. Only a varies, even though the function itself still has two inputs.
Often times, you would like to loop over a single parameter of a function that has multiple arguments:
For example, you would like to fix the value of b at 5 while trying different values of a of the following function:
Then you can do this:
As you can see, this function tries each of 1:10 (called internally x), give it to square_them_add() as its first argument while b is fixed at 5.
Two exercises. The first is the lapply version of the cubing loop you just wrote, so you can feel the difference directly, including the fact that you get a list back rather than printed output. The second is the one to slow down on. You need a function that takes both a nitrogen rate and a corn price and returns revenue, and then you loop over corn prices while holding nitrogen at two hundred. Look carefully at how the answer names its arguments in the lapply call. That is not stylistic. Pass them positionally in the wrong order and the corn price gets treated as the nitrogen rate, the code runs perfectly happily, and you get nonsense.
For Exercise 1, seq from five to fifty by five creates ten inputs and the shorthand function returns x caret three for each, so the answer needs neither an explicit for-loop nor print. For Exercise 2, calc revenue first computes yield as one hundred twenty plus twenty-five times the natural log of N, multiplies by corn price, and returns revenue. Corn price seq contains two point five through four point zero in tenths. The final lapply sends each price in as x, explicitly fixes N at two hundred, and returns a list of sixteen revenues. The callout’s bad positional example, calc revenue of x and two hundred, is dangerous precisely because it is valid R: it silently assigns the varying prices to N and two hundred to corn price. The folded answer chunks do not evaluate, so use them to inspect the solution after attempting the live workspace rather than relying on them to create objects.
Use lapply() to cube each element of a sequence starting at 5, increasing by 5, and ending at 50.
Work here
Answer
Define a function that takes a nitrogen rate and corn price as arguments and calculate revenue. Yield can be estimated using the following equation:
\text{corn yield} = 120 + 25 \times log(\text{nitrogen rate})
At each value of the corn price sequence of seq(2.5, 4.0, by = 0.1), calculate the revenue using lapply() where nitrogen rate is fixed at 200 (lb/acre).
Work here
Answer
#--- takes BOTH arguments and returns revenue, as the question asks ---#
calc_revenue <- function(N, corn_price) {
yield <- 120 + 25 * log(N)
revenue <- corn_price * yield
return(revenue)
}
corn_price_seq <- seq(2.5, 4.0, by = 0.1)
#--- loop over corn price, hold nitrogen fixed at 200 ---#
lapply(corn_price_seq, \(x) calc_revenue(N = 200, corn_price = x))Name the arguments
Naming them (N = 200, corn_price = x) is what stops the looped value landing in the wrong slot. With positional arguments, calc_revenue(x, 200) would quietly compute yield at a nitrogen rate of 2.5 and multiply it by a corn price of 200.
So far every loop has walked down a single list of values. Real questions are rarely that tidy. The example here is a sensitivity analysis: you want to know how profitable corn production is across a range of corn prices and a range of nitrogen rates, which means every combination of the two, not each list separately. The trick is in the last two bullets, and it is worth remembering because it generalises to any number of variables. Build a data frame that holds all the combinations, one combination per row, and then loop over the rows. You have turned a two-dimensional problem back into a one-dimensional loop.
The example we have looked at is a very simple case where a loop is done over a single list of values
It is often the case that you want to loop over multiple variables.
Example
You are interested in understanding the sensitivity of the profitability of corn production with respect to corn price and nitrogen application rate.
So, you would like to loop over two sets of sequences of values:
How
The trick is to
create a data.frame of two (or as many variables as you would like to loop over) variables (corn price and nitrogen application rate), which stores all the permutations of the two variables
then loop over the rows of the data.frame
Concretely, we are asking how corn revenue responds to two things at once: the price of corn, between three and five dollars a bushel, and the nitrogen rate, between zero and three hundred pounds an acre. Keep those ranges in mind as we go through the four steps, because the numbers on the next few tabs come straight from them. And notice this grid is small, deliberately. Three prices and four nitrogen rates give twelve combinations, few enough that you can look at the whole table and check it by eye. Get the mechanics right on twelve rows before you turn it loose on twelve thousand.
We are interested in understanding the sensitivity of corn revenue to corn price and applied nitrogen amount.
We consider
Four steps, and they are worth naming before you look at the code. Step one builds the sequences and expands them into every combination with expand dot grid, giving you a data frame with one row per case. Step two writes a function that takes a row number, pulls the parameters out of that row, and computes what you want. Step three loops over the row numbers. Step four stacks the results back into one data frame with bind rows. That shape, expand then loop over rows then bind, is one you will reuse constantly, so learn it as a pattern rather than as four separate bits of code.
On Step 1, seq from three to five by one creates corn prices three, four, and five, while seq from zero to three hundred by one hundred creates nitrogen rates zero, one hundred, two hundred, and three hundred. Expand dot grid crosses corn price vec and nitrogen vec into twelve rows under the names corn price and nitrogen. The assignment is wrapped in parentheses so the table is both stored as parameters data and printed, and autorun makes those inputs available to later tabs.
On Step 2, gen rev corn of i treats i as a row position. Square-bracket indexing pulls corn price and nitrogen from that row by column name. The exponential response formula produces yield, price times yield produces revenue, and data frame packages the two inputs and their result into one row before return sends it back. Autorun defines the function without requiring you to visit the tab first.
On Step 3, one colon nrow of parameters data creates all twelve row numbers and lapply calls the named function once for each. The result rev data is therefore a list of twelve one-row data frames. Piping it to head shows the first six list elements for a quick structural check. On Step 4, bind rows of rev data stacks those matching one-row data frames and the surrounding parentheses both save and display final results. Flip across all four tabs and follow one parameter row from its construction to its final row in that combined table.
Get a sequence of values for corn price and nitrogen rate:
We then create a complete combination of the values using the expand.grid() function. The resulting object is a data.frame.
Define a function that
parameters_data to extract the parameters stored at the row numberThis function
i (act as a row number within the function)ith row of parameters_data data.frame of the resulting revenue, corn price, and nitrogen ratedata.frameDo a loop using lapply():
Combine the list of data.frames into a single data.frame using bind_rows() from the dplyr package.
This is the most practical advice in the chapter, so do not skip past it. Do not start by writing a function. Start by writing ordinary code that works for one specific case, here row one, with the number one written in literally. Run it, look at the result, satisfy yourself it is right. Only then turn it into a function by replacing that one with i. The reason is that debugging is much harder inside a function, where the intermediate objects disappear the moment it finishes. Working line by line at the top level, you can inspect everything as you go. Get it right first, then wrap it up.
The code on screen is Step 2 with the wrapper removed. The two square-bracket expressions read row one of parameters data, first from corn price and then from nitrogen. The next lines calculate yield and revenue, and data frame gathers the inputs and answer as data to return. At this stage you can print any one of those objects and catch a wrong extraction or formula immediately. Once every line works, place the code inside function of i, replace both literal row numbers with i, and return data to return. That mechanical conversion keeps the tested calculation unchanged.
Before defining a function, write code that works for one row.
We will work on a specific value of i. Here it is i = 1.
After you confirm the code you write gives you desired outcomes, make it a function by replacing 1 with i.
Now do the whole pattern yourself. Four steps, laid out for you: build the two price sequences, expand them into all combinations, write a function that takes a row number and returns profit for that combination, and loop. Notice it is profit here, not revenue, so you subtract the nitrogen cost as well as counting the corn income. And note that the nitrogen rate is fixed at two hundred throughout, so it is not one of the things you are looping over, even though it appears in both equations. Work through it in the order given, and use the tip from the last tab: get row one working before you write any function at all.
In the folded answer, the corn-price sequence runs from two point five to four point zero by zero point zero five, and the nitrogen-price sequence runs from zero point two to zero point six by zero point zero one. Expand dot grid stores every price pair in cases data. Inside get profit of i, square brackets extract both prices from row i; log of two hundred uses the fixed nitrogen rate in the yield equation, and n price times two hundred subtracts that fixed quantity’s cost. The returned one-row data frame carries both prices beside profit so the result remains identifiable. Finally, one colon nrow of cases data supplies every row number, the shorthand function calls get profit of i, and the pipe into bind rows converts the resulting list into one table. The answer chunk is folded and not evaluated, so it demonstrates the complete pattern without doing the exercise for your live workspace.
Find the profit of corn production at different price combinations of corn and nitrogen where nitrogen rate is fixed at 200 lb/acre.
seq(2.5, 4.0, by = 0.05)seq(0.2, 0.6, by = 0.01)\text{corn yield} = 120 + 25 \times log(\text{nitrogen rate})
\text{profit} = \text{corn price} \times \text{corn yield} - \text{nitrogen price} \times \text{nitrogen rate}
Answer
corn_price_vec <- seq(2.5, 4.0, by = 0.05)
#--- CHANGED: was by = 0.02, but the instruction asks for by = 0.01 ---#
n_price_vec <- seq(0.2, 0.6, by = 0.01)
cases_data <-
expand.grid(
corn_price = corn_price_vec,
n_price = n_price_vec
)
get_profit <- function(i) {
corn_price <- cases_data[i, "corn_price"]
n_price <- cases_data[i, "n_price"]
corn_yield <- 120 + 25 * log(200)
profit <- corn_price * corn_yield - n_price * 200
data_to_return <-
data.frame(
corn_price = corn_price,
n_price = n_price,
profit = profit
)
return(data_to_return)
}
lapply(
1:nrow(cases_data),
\(i) get_profit(i)
) %>%
bind_rows()Having spent all that time teaching you to loop, here is the twist. You should not have looped in any of those examples. Every one of them could be vectorized. A vectorized operation takes whole vectors as input and works on all the elements at once, which is what the addition example on screen is doing. Adding two thousand-element vectors is a single operation as far as your code is concerned. R is built for this, and the internals are compiled code, so it is not just tidier to write, it is genuinely faster. The next few tabs show you exactly how much faster.
The example creates x and y as the integers one through one thousand, so they have equal length and corresponding positions. The line z vec left arrow x plus y performs element-wise addition: the first element of x is added to the first of y, the second to the second, and so on, then the left arrow stores all one thousand sums in z vec. Here “in parallel” describes one vectorized expression over matching elements, not a future-based calculation distributed across CPU cores. Move to Compare to see the deliberately looped equivalent and confirm that its values match.
Actually, we should not have used a for loop or lapply() in any of the examples above in practice.
This is because they can be easily vectorized.
Vectorized operations are those that take vectors as inputs and work on each element of the vectors in parallel
Example
Two ways to add the same two vectors. The vectorized version is x plus y, and that is the entire program. The loop version walks the indices one at a time, pulling out element i from each vector, adding them, and building a list, which then has to be unlisted back into a vector. Then all dot equal confirms the two results are identical, and they are. So you get exactly the same answer either way. The question is what each one costs, and the sentence at the bottom is the hint. R is written to be good at the first kind of operation and comparatively bad at the second.
Vectorized
Non-vectorized (loop)
Compare
Both produce the same results. However, R is written in a way that is much better at doing vectorized operations.
Now we measure instead of asserting. Microbenchmark runs each expression a hundred times and reports the distribution, which is what you want, because a single timing is mostly noise. Look at the numbers when the block finishes. The vectorized version is dramatically faster. The explanation is in the second bullet, and it is worth understanding rather than memorising. Every trip round that loop, R has to work out what the objects are, check their types, allocate somewhere for the result, and so on. Vectorized code pays those costs once for the whole vector instead of once per element. That overhead is the entire difference.
The quoted labels inside microbenchmark name the two rows of timing output. The vectorized expression adds the existing x and y vectors directly. The non-vectorized expression constructs indices one through one thousand and calls a function for every index to add x bracket i and y bracket i. As the note says, this benchmark deliberately leaves off unlist, so the loop is not also charged for converting its list to a vector. Times equals one hundred repeats each candidate one hundred times, and unit equals milliseconds reports the summaries in milliseconds. Even with the conversion cost removed, the repeated R-level calls make the loop slower.
Let’s time them using the microbenchmark() function from the microbenchmark package.
Here, we do not unlist() after lapply() to just focus on the multiplication part.
The simplest possible illustration of the point. On top, lapply calling square it a thousand times. Underneath, square it applied to the whole vector at once. Same answers. But look at why the second one works. The body of square it is x squared, and squaring in R is already vectorized, so the function inherits that for free without you doing anything. That is the general lesson. A great many functions you write are automatically vectorized because everything inside them is, and you never needed the loop in the first place. Before you write lapply, try handing your function the whole vector and see what happens.
There is still a container difference to notice. Lapply of one colon one thousand and square it returns a list with one squared value per element, while square it of one colon one thousand returns one numeric vector because the caret operator accepts the whole vector. The numerical values agree, but the vectorized result is already in the simpler form most later calculations expect. That removes both the repeated function-call overhead and the need to unlist the answer.
Instead of this:
You can just do this:
Now the same idea on the real example. Compare this function with the one from the loop version. That one took a row number and went and looked up the parameters itself. This one just takes corn price and nitrogen directly as arguments and does the arithmetic. That change is precisely what makes it vectorizable. Then instead of looping over rows, we call it once inside mutate, handing it the two entire columns, and it computes every row’s revenue in one go. No loop, no list of data frames, no bind rows at the end. One line, and the result is already a proper data frame.
Inside gen rev corn short, exp of zero point four minus zero point zero two times nitrogen is evaluated element by element, the yield vector is multiplied by the matching corn-price vector, and return revenue sends that whole vector back. Autorun defines the function for the next block. Mutate of parameters data keeps the existing parameter columns and adds the returned vector as a new revenue column, while the left arrow stores the finished table as rev data two. The two input columns have equal length and correspond row by row, which is why each revenue is attached to the correct price and nitrogen combination.
Here is the vectorized version of the revenue sensitivity analysis:
Then use the function to calculate revenue and assign it to a new variable in the parameters_data data.
And the timing comparison for the realistic case. Run it and look at the gap. The vectorized version computes every combination in a single pass. The loop version calls a function once per row, and each of those calls does its own data frame lookup and builds its own little one-row data frame, which bind rows then stitches together. All of that is overhead you simply do not pay in the vectorized version. Keep this result in mind for the parallel section that follows, because the temptation there is to reach for more cores. Very often the better move is to stop looping at all.
Inside microbenchmark, the quoted labels identify the two timing rows. The vectorized candidate mutates parameters data once with the whole-column revenue function. The non-vectorized candidate uses seq len of nrow of parameters data to generate valid row indices, passes each to gen rev corn, and binds the returned data frames. Both assign their answer to rev data, so neither is timed for console printing. Times equals one hundred repeats each complete candidate one hundred times, and unit equals milliseconds reports milliseconds. The outputs are comparable because both produce the same columns for the same parameter combinations; the benchmark isolates the cost of how they get there.
Let’s compare the vectorized and non-vectorized version:
Parallelization means splitting work across several cores so that pieces run at the same time. Our focus is the easy case, which has the wonderful name embarrassingly parallel: every piece is completely independent, so no piece needs any other piece’s answer. Squaring integers qualifies, because computing one squared does not require knowing two squared. That independence is what makes parallelization straightforward, since you never have to think about ordering or about one worker waiting on another. And the last bullet is genuinely reassuring. Most of what you will actually want to speed up, simulations, or running a model over many datasets, does fall into this category.
Parallelization of computation involves distributing the task at hand to multiple cores so that multiple processes are done in parallel.
Our focus is on the so called embarrassingly parallel processes.
Embarrassingly parallel process: a collection of processes where each process is completely independent of any another (one process does not use the outputs of any of the other processes)
The example of integer squaring is embarrassingly parallel. In order to calculate 1^2, you do not need to use the result of 2^2 or any other squares.
Embarrassingly parallel processes are very easy to parallelize because you do not have to worry about which process to complete first to make other processes happen.
Fortunately, most of the processes you are interested in parallelizing fall under this category
The good news is that you have already learned the hard part. Parallelizing with the future dot apply package means writing future underscore lapply where you would have written lapply, and that is essentially the whole change. Look at the two lines at the bottom: same values, same function, same shape of result. Everything you know about lapply carries straight over, including the fact that you get a list back. There is one piece of setup to do first, which the next tab covers, and it is the piece people get wrong, so pay attention to it.
The installation line is a one-time setup that downloads future dot apply; library of future dot apply is the per-session step that makes future underscore lapply available. Those lines are marked not to evaluate in the rendered deck, so run them on your own computer when needed. In the comparison, both calls loop over one through one thousand and apply the anonymous squaring function. The serial result is stored in sq ls, and the future-based result in sq ls par. Replacing the function name does not itself choose the number or kind of workers, so move to Preparation for the plan that controls where those future calls run.
We will use the future_lapply() function from the future.apply package for parallelization.
Using the package, parallelization is a piece of cake as it is basically the same syntactically as lapply().
How
You can simply replace lapply() with future_lapply()!
Two things before you can parallelize. First, find out how many cores you have, with detect cores. Second, and this is the part that matters, tell R how to run the workers, using plan. We use multisession, which starts separate R sessions as workers and works on every platform. Now read the red callout, because it is the trap. There is another option, multicore, that uses forking. Forking does not exist on Windows, and it is switched off inside RStudio. If you ask for multicore where forking is unavailable, R does not warn you and does not switch to something else. It quietly runs everything on one core. Your code works, and you gain nothing.
On the Cores and plan tab, parallel double colon detect cores uses the double colon to call the function without attaching the parallel package. In plan of multisession, workers equals parallel double colon detect cores minus one, multisession chooses separate R processes and workers sets their count. Subtracting one leaves one detected core out of the worker pool so the computer retains some capacity for the operating system and your interactive session. The note underneath separates the plan from ordinary R code: setting a parallel plan does not make every calculation parallel. Only future-aware calls such as future underscore lapply use it, while ordinary work remains single-core.
Then flip to Which backend. Sequential means no parallel workers, multisession means portable separate sessions, and multicore means forked sessions with lower startup cost but the platform restrictions in the important callout. The critical behavior is already stated on screen: unavailable forking silently becomes sequential rather than switching to multisession. That is why multisession is the course default even on a machine where multicore might sometimes be faster.
You can find out how many cores you have available for parallel computation on your computer using the detectCores() function from the parallel package.
Before we implement parallelized lapply(), we need to declare what backend process we will be using by plan().
Note
Unless you tell R explicitly to parallelize things (like using future_lapply()), R always uses a single core by default. So, you do not have to change anything manually when you do not want to use multiple cores.
sequential: this is just a regular loop, no parallelizationmultisession: separate R sessions as workers. Works everywheremulticore: forked sessions. Faster to start than multisession, but only on Mac/Linux, and not inside RStudioWhy not multicore?
multicore relies on forking, which is unavailable on Windows and is switched off inside RStudio. When forking is unavailable, plan(multicore) does not switch to multisession for you — it quietly runs everything sequentially in a single process, with no error and no warning.
You would see your code run, get correct answers, and gain nothing at all. multisession is the safe default. Reach for multicore only when you know you are on Mac/Linux outside RStudio and startup cost matters.
So here it is. Future underscore lapply, squaring a thousand numbers, spread across your workers. Run it. The point of this tab is that nothing dramatic happens. You get a list of a thousand squares back, exactly as lapply would have given you, and it looks completely ordinary. That is the whole idea. Once the plan is set, the parallel version of your code is visually indistinguishable from the serial version, which means you can switch between them freely. Which raises the obvious question, and it is the one the next tab asks, so form an opinion before you click. Was any of this actually faster?
Now we time it, and the result is probably not what you expected. The parallelized version comes out slower than the plain lapply. Not slightly slower, noticeably slower. This is a real result, and it is worth sitting with for a moment rather than assuming something is broken on your machine. Nothing is broken. Parallelization is not free, and this particular example is constructed to show you the cost rather than the benefit. Have a guess at where the time is actually going before you move on, and then the next tab will tell you whether you were right.
The benchmark gives both expressions the same one thousand inputs and the same anonymous function, so the backend is the meaningful difference. Each result is assigned to sq ls, preventing console printing from becoming part of one timing but not the other. Times equals one hundred repeats both versions one hundred times and unit equals milliseconds reports milliseconds. The chunk is cached because starting workers and repeating the parallel benchmark is expensive during deck rendering; caching reuses the previously computed output until the code or its dependencies change. The timing you see is still the recorded result of the code, not a claim that parallel work is always slower.
microbenchmark(
#--- parallelized ---#
"parallelized" = {
sq_ls <- future_lapply(1:1000, function(x) x^2)
},
#--- non-parallelized ---#
"not parallelized" = {
sq_ls <- lapply(1:1000, function(x) x^2)
},
times = 100,
unit = "ms"
)Unit: milliseconds
expr min lq mean median
parallelized 1833.350588 1868.9578170 1891.6517737 1886.3670730
not parallelized 0.161581 0.1735325 0.2321408 0.1784525
uq max neval cld
1905.634100 2006.153001 100 a
0.189461 4.672688 100 bHere is the explanation. Handing a job to another core is itself work. The values have to be sent over, the worker has to be told what to do, and the answers have to come back. With the default scheduling, that communication overhead is generally paid once per future chunk, with about one chunk per worker, while the function-call and result-building overhead remains for every element. Squaring a number takes almost no time at all, so the overhead completely swamps the saving. The rule that follows is in the last bullet, and it is the one to remember. Parallelization pays when each individual iteration is slow. If your iterations are fast, adding cores makes things worse.
This is because communicating jobs to each core takes some time as well.
So, if each of the iterative processes is super fast (like this example where you just square a number), the time spent on communicating with the cores outweighs the time saving due to parallel computation.
Parallelization is more beneficial when each of the repetitive processes takes long.
So let us do it properly, with work that is actually slow. Monte Carlo simulation is close to the ideal use case: you run the same procedure many times on freshly generated data, and each run is completely independent of the others. The specific question here is a familiar one from econometrics. If an independent variable is correlated with the error term, does that bias your OLS estimate? We know the answer is yes, and that is deliberate. The point is not the finding, it is watching how long the simulation takes and how much of that time parallelization can give you back.
One of the very good use cases of parallelization is MC simulation
We will run MC simulations that test whether the correlation between an independent variable and error term would cause bias (yes, we know the answer).
Three steps, repeated a thousand times. Step one generates fifty thousand observations, and the structure matters. There is a shared term mu that goes into both x and the error v. That shared component is exactly what creates the correlation between the covariate and the error, which is the whole point of the exercise. Step two runs OLS and keeps the coefficient on x. Step three repeats the whole thing a thousand times, so we can look at the distribution of those estimates rather than a single draw. And note the last line: each repetition is independent, so this is embarrassingly parallel.
y = 1 + x + v
where + \mu \sim N(0,1) + x \sim N(0,1) + \mu + v \sim N(0,1) + \mu.
The \mu term causes correlation between x (the covariate) and v (the error term).
estimate the coefficient on x via OLS, and return the estimate.
repeat this process 1,000 times to understand the property of the OLS estimators under the data generating process.
This Monte Carlo simulation is embarrassingly parallel because each process is independent of any other.
Here are steps one and two as a function for one repetition, and notice it follows the advice from earlier in the chapter exactly. Step three happens later when you apply this function a thousand times. It takes i, a single index, and everything else it needs it creates internally. The sample size is set inside, the data is generated fresh on each call, the regression is run, and it returns just one number, the coefficient on x. Returning one number rather than the whole model object matters more than it looks. When you run this a thousand times in parallel, everything you return has to be sent back from the worker, and shipping a thousand full model objects around would cost you real time.
The index i is accepted so an apply function can call M C sim once per repetition, even though the data-generating calculation does not otherwise use its value. N left arrow fifty thousand fixes the sample size. The three r norm of N calls generate standard-normal draws, with the same generated mu vector added to both x and v; y left arrow one plus x plus v then implements the displayed model. Data table of y and x creates the regression data, and l m of y tilde x with data equals data estimates an intercept and slope from those named columns. Finally, reg dollar-sign coef bracket x selects only the fitted slope by name. Selecting by name makes the intended coefficient explicit and avoids depending on its numeric position in the coefficient vector.
Here is a function that implements steps 1 and 2 for one repetition. The later apply call implements step 3:
#--- one repetition of steps 1 and 2 ---#
MC_sim <- function(i) {
N <- 50000 # sample size
#--- steps 1 and 2: ---#
mu <- rnorm(N) # the common term shared by both x and u
x <- rnorm(N) + mu # independent variable
v <- rnorm(N) + mu # error
y <- 1 + x + v # dependent variable
data <- data.table(y = y, x = x)
#--- OLS ---#
reg <- lm(y ~ x, data = data) # OLS
#--- return the coef ---#
return(reg$coef["x"])
}Now the numbers, and this is the payoff for the whole section. A single run takes a small fraction of a second, which tells you each iteration is doing real work, unlike the squaring example. Then the sequential version of all one thousand runs, and then the parallel version with parallel-safe random-number streams. Compare those last two figures. That is a substantial speed-up, achieved with a small change to one line of code. Two caveats worth holding on to. These timings depend on how many cores this machine has, so yours will differ. And the gain is nowhere near the number of cores, because of the handover cost we just discussed.
In every block, tic starts the stopwatch and toc stops it and reports elapsed time. The first call uses index one only to satisfy the function interface and stores its single coefficient as single res. The sequential block sends indices one through one thousand to lapply, while the parallel block sends the same indices and function to future underscore lapply, so their M C results objects have the same list structure. Future dot seed equals true asks the future framework to create valid, reproducible random-number streams for the parallel iterations instead of letting workers accidentally reuse or mishandle random state. Both thousand-run chunks are cached so the rendered lecture can reuse these expensive results. That cache changes rendering time, not the algorithm being compared.
Single run:
Not parallelized (sequential):
Parallelized:
One alternative worth knowing about if you are on a Mac or Linux. The parallel package’s mclapply does the same job and is just as easy, since its syntax also matches lapply, with the number of cores set through mc dot cores. And pbmclapply is the same thing with a progress bar, which sounds trivial but is genuinely valuable when a job runs for an hour and you want to know whether it is halfway or nearly done. The catch is in the tab title. Both rely on forking, so neither works on Windows. If you are writing code other people will run, stay with future underscore lapply.
The first library call attaches the base-recommended parallel package, and m c lapply of one colon one thousand, M C sim, and m c dot cores equal detect cores minus one applies the same simulation to all one thousand indices while leaving one detected core outside the worker count. The second library attaches p b m c lapply and changes only the function name, preserving the inputs, simulation, and core setting while adding progress feedback. The entire chunk has evaluation turned off because the deck must remain renderable on platforms where fork-based execution is unavailable. Treat it as code to copy to a suitable Mac or Linux session, not as a benchmark produced by this slide.
For Mac or Linux users, parallel::mclapply() is just as compelling (or pbmclapply::pbmclapply() if you want to have a nice progress report, which is very helpful particularly when the process is long).
It is just as easy to use as future_lapply() because its syntax is the same as lapply().
You can control the number of cores to employ by adding mc.cores option. Here is an example code that does the same MC simulations we conducted above:
A new problem, and a very practical one. You have run a nitrogen trial, fitted a yield response curve, and now you want the nitrogen rate that maximises profit. The profit expression is on screen: corn price times yield, minus nitrogen price times the amount you applied. The distinction in the last line is the one to hold on to for the rest of this section. N is a decision variable, the thing you get to choose. The two prices are parameters, things the world hands you. Optimization means choosing the decision variable for given parameters, and the interesting work comes from repeating that across many parameter combinations.
The yield function is one hundred twenty plus twenty-five times the natural log of nitrogen rate, so nitrogen must be positive wherever you evaluate it. That is why the grids on the following tabs begin at zero point one rather than zero. In the objective, P C is dollars per bushel and multiplies yield in bushels per acre to produce revenue per acre. P N is dollars per pound and multiplies nitrogen in pounds per acre to produce nitrogen cost per acre. Subtracting cost from revenue gives the profit being maximized, with the price parameters held fixed during each optimization.
Suppose you have run a randomized nitrogen experiment for corn production on a field, collected data, and run a regression to find the following quantitative relationship between corn yield (bu/acre) and nitrogen rate (lb/acre):
\text{corn yield} = 120 + 25 \times log(\text{nitrogen rate})
You are interested in finding the best nitrogen rates that maximize profit at different combinations of corn and nitrogen prices for this field.
Max_{N}\;\; P_C \times Y(N) - P_N \times N
Here, N is a decision variable, and P_C and P_N are parameters.
Grid search is the least clever method available and often the right one. You simply evaluate the objective at a lot of values of the decision variable and keep the best. That is it. It is inefficient, because you compute a great many values you end up discarding, and it does not scale: with one or two decision variables it is fine, with ten it is hopeless. But for a problem like this it is reliable, easy to reason about, and impossible to get subtly wrong in the way a numerical optimizer can. The table shows the idea, profit computed at every nitrogen rate from nearly zero up to three hundred.
Read the pipeline below the definition. Seq from zero point one to three hundred by zero point one creates three thousand candidate nitrogen rates and starts above zero because the yield formula contains log of N. Data table of N stores those candidates in a column. The first mutate evaluates yield at every candidate, and the second computes profit using a fixed corn price of three point five and nitrogen price of zero point four. The tenth-of-a-pound spacing controls the grid’s resolution: a finer step can locate the maximum more precisely but creates more rows to evaluate. This output-context cell runs automatically and presents the resulting table, so the next tabs can focus on turning the same calculation into reusable code.
Grid search is a very inefficient yet effective tool for finding solutions to optimization problems as long as the dimension of the optimization problem is low (1 or 2 variables).
Grid search basically evaluates the objective function (profit here) at many levels of the decision variables (nitrogen here) and pick the one that maximizes the objective function (profit).
Example
Now in code. First a function taking a nitrogen rate and the two prices, returning profit, which is just the equation from the setup tab written out. Then we build a data frame of nitrogen rates from nought point one up to three hundred in steps of a tenth, and compute profit at every one of them with mutate. Notice what we are not doing here. There is no loop. Get profit is vectorized because everything inside it is, so handing it the whole column works directly. That is the lesson from the vectorization section showing up in real work, and it is why this runs instantly on three thousand rows.
The argument names mirror the mathematical notation: N is the candidate rate, P C the corn price, and P N the nitrogen price. Inside, P C times yield is revenue and P N times N is nitrogen cost, so the subtraction returns profit. Autorun makes get profit available to the later tabs. In the second block, tibble creates the nitrogen column and get profit of N, three point five, and zero point four holds both prices fixed while the whole N column varies. Mutate adds the resulting profit column, the left arrow saves the table as data main, and the surrounding parentheses display it at the same time. This table is the grid search: each row is one candidate and its objective value.
Let’s define a function that takes N, P_C, and P_N values and returns profits.
Let’s create a sequence of N values and then calculate profit at each level of N.
Here is the same information as a picture, and it is worth a moment. Profit rises steeply at low nitrogen rates, flattens out, and then begins to fall as the cost of the extra nitrogen outweighs the extra yield it buys. The dashed red line marks the maximum. Two things to take from the shape. The peak is quite flat, which in practice means being somewhat wrong about the optimal rate costs you very little, and that is genuinely reassuring for a farmer. And the curve is smooth with a single peak, which is exactly the situation where grid search cannot go wrong.
The first line finds the maximizing row or rows by comparing every profit with max of profit, then dollar-sign N extracts the corresponding rate into opt N. G g plot of data main supplies the grid-search table to all layers. In geom line, aes of y equals profit and x equals N maps the two columns to the axes and connects the evaluated candidates. Geom v line places the reference line at opt N; line type two makes it dashed and color equals red makes it red. Finally, annotate text writes the label at an explicit location, forty units to the left of the optimum and at a y value of four hundred fifty, with text size three. The output-context cell runs automatically so what you see is the rendered diagnostic figure rather than another result table.
Here is the visualization of profit-N relationship:
Two ways to pull the best row out. The first filters for the row where profit equals the maximum, which reads exactly like what you mean. The second sorts the whole dataset by profit and takes the last row. Both give you the same answer, but read the warning underneath, because this deck previously had it backwards. The filter version is roughly three times faster, and it has to be. Filtering scans the column once to find the largest value, whereas arranging sorts everything first, and sorting is fundamentally more expensive than scanning. Use arrange when you want the top several rows, or you actually need the ordering. For a single maximum, filter.
The dplyr double-colon prefixes call each function from dplyr without relying on the package search path. In the first block, max of profit is calculated from data main, and filter retains every row tied at that value rather than arbitrarily choosing one. In the second, arrange of data main by profit uses its default ascending order, the pipe passes that sorted table forward, n evaluates to the number of rows, and slice of n selects the final one. That distinction matters if the grid contains a tie: filtering returns all maximizing rows, while sorting and slicing returns one last row according to the current order. The code track setting of zero point five gives the editable code and its live output equal shares of the side-by-side cell; it changes the display layout, not either method’s result. The warning then tells you why the clearer filter is also the cheaper choice here.
Once the profit-N relationship is found, we can use dplyr::filter() combined with max() to identify the optimal N rate.
Alternatively, you can sort the data by profit in the ascending order (default) and pick the last row using dplyr::slice(n()).
Convenient, but slower
This second method is slower, not faster — about 3x on this data.
filter() only has to scan the column once to find the maximum. arrange() has to sort the whole dataset first, and sorting costs more than scanning. The gap widens as the data grows.
Use arrange() when you want the top few rows, or the ordering itself. To pick out a single maximum, filter() is both clearer and cheaper.
You can now find the optimal nitrogen rate for one price combination. The obvious next step is to do it for many, and here is the strategy, which should look familiar because it is the same three-step pattern from the loop section. Build a data frame of all the price combinations you care about. Write a function that takes a row number, pulls that row’s prices out, and runs the grid search you just wrote. Then loop over the row numbers. Notice that the grid search itself does not change at all. It just gets wrapped. That is what makes this pattern so reusable.
Now that you have written codes to find the optimal N at a given combination of corn and nitrogen prices.
We can move on to the next step of finding the optimal N rates at many various combinations of corn and nitrogen prices.
Here is the coding strategy:
Define a set of all the combinations of corn and nitrogen prices you want to analyze as a data.frame.
Define a function that extract corn and nitrogen prices from the parameter data.frame and find the optimal N rate at the given combination of prices.
Loop over the price combinations (loop over the rows of the data.frame created in step 1).
Step one, the parameter grid. Three corn prices and three nitrogen prices, expanded into all nine combinations, one per row. Two small things. We pipe the result through tibble, which does not change the data but gives you nicer printing, particularly once the grid gets long. And keep the size in mind. Nine rows is small enough to check by hand, which is exactly what you want while you are still getting the code right. When this works, widening the sequences is a one-character change, and none of the rest of the code has to care.
Specifically, seq from two point five to four point five by one gives corn prices two point five, three point five, and four point five, while seq from zero point two to zero point six by zero point two gives nitrogen prices zero point two, zero point four, and zero point six. Naming the arguments inside expand dot grid names the resulting columns P C and P N, which matches the arguments of get profit. The assignment is wrapped in parentheses so price parameters is stored and displayed together. Autorun is important here because later tabs read this exact nine-row object even if you have not manually executed Step 1.
Here, we define a data.frame of parameters to be explored. We will be looping over the rows of the data.frame.
Step two is the function, and look at what is inside it. It pulls the corn and nitrogen price out of row i, builds the nitrogen grid, computes profit at every rate using the get profit function from earlier, and filters down to the best row. Then the two mutates on the end are easy to overlook but important. They attach the prices onto the result. Without them you would get back a set of optimal nitrogen rates with no record of which price combination each one belongs to. When you are looping and stacking results, always carry the parameters through alongside the answer.
The expressions price parameters bracket i comma blank bracket dollar-sign P C and dollar-sign P N first select row i and then extract the named price from that one-row result. N data supplies the same three thousand positive candidate rates used earlier. In the pipeline, mutate hands the whole N column plus the two scalar prices to vectorized get profit, filter where profit equals max of profit keeps the maximizing candidate or candidates, and the two later mutates repeat the fixed prices beside the answer. Autorun defines get opt N for Step 3. The separate get opt N of one block is the single-case test from the chapter’s debugging advice: it should return the optimum for the first price row before you trust the function inside a loop.
Now, we will define a function that extract a combination of corn and nitrogen prices from price_parameters (extract a row from price_parameters), and then find the optimal N.
Check if this function works:
Step three, the loop, and we use future underscore lapply so the price combinations can run across cores. Then bind rows stacks the nine one-row data frames into a single table, the same finish as the earlier loop example. One honest caveat, in the note underneath. These slides run R inside your browser, which is single-threaded, and we never called plan here, so on this page it runs one combination after another regardless. Copy it onto your own machine, set a plan first, and the identical line parallelizes. That interchangeability is precisely the point, but you will not see the speed-up here.
The input one colon nrow of price parameters is the sequence of all nine row numbers, and passing get opt N by name tells the apply function what to run for each. The returned opt N all l s is a list because future underscore lapply preserves lapply’s output contract. Double brackets one, spoken as the first list element, extracts the first returned data frame itself so you can inspect its columns and values. Bind rows of opt N all l s then stacks every element, stores the combined table as opt N all, and the surrounding parentheses display it. The callout is therefore about execution mode only: sequential WebR and planned multisession R produce the same list and final table.
Loop:
In this browser, this runs sequentially
These slides run R in your browser through WebR, which is single-threaded, and plan() was never called here. future_lapply() still works — it just falls back to running one job after another.
Copy this onto your own computer, call plan(multisession, workers = ...) first, and the same line parallelizes with no other change. That interchangeability is the point.
Combine the list of data.frames into a single data.frame using bind_rows().
Now the same problem without looping at all, and the trick is a change of perspective. In strategy one the nitrogen grid lived inside the function and got rebuilt for every price combination. Here we put nitrogen into the parameter grid itself, alongside the two prices, so a single data frame holds every combination of all three. Then profit is one vectorized calculation over the whole thing, and finding the best rate becomes a grouped operation. No function, no loop, nothing to bind together at the end. Read the three steps and notice that nitrogen rate has been promoted from something the function handles to just another column.
Instead of writing a loop like above, we can actually vectorize the process. Here are the steps:
Define a set of all the combinations of nitrogen rate, corn price, and nitrogen price you want to analyze as a data.frame.
Calculate profits for all the combinations of nitrogen rate, corn price, and nitrogen price inside the data.frame
Find the optimal N rate for each combination of corn price and nitrogen price
Here is that combined grid. Three corn prices, three nitrogen prices and three thousand nitrogen rates, expanded into every combination, which gives twenty-seven thousand rows. Compare that with strategy one, where the grid held only the nine price combinations and the nitrogen rates were generated inside the function, over and over again. Everything is now in one table up front. This is also where the memory question starts to bite, and we come back to it shortly. This grid is still manageable, but the row count is the product of all three lengths, so it grows very quickly once the sequences get larger.
The three seq calls define those exact dimensions, including a zero point one nitrogen step that avoids evaluating the logarithm at zero. Expand dot grid with P C, P N, and N names all three columns and forms their Cartesian product. Piping to tibble improves printing, and arrange by P C and P N places rows with the same price pair together without changing which combinations exist. The result is stored under the distinct name eval parameters, so it cannot overwrite Strategy 1’s smaller price parameters table. Parentheses display the assigned grid, and autorun makes it ready for the calculation on the next tab.
Here, we define all the combinations of nitrogen rate, corn price, and nitrogen price you want to analyze as a data.frame.
Two lines to finish. First, profit for every row at once with mutate, handing get profit three whole columns. Second, group by the two prices and, within each group, take the row with the highest profit. That grouping is doing the same job the loop did in strategy one: it treats each price combination separately. But there is no function to write and nothing to stack at the end, because the answer is already a single data frame. Compare this with the three tabs it took to do the same thing by looping, and you can see why vectorizing is worth reaching for first.
The first parenthesized expression keeps all twenty-seven thousand candidate rows, adds the vector returned by get profit of N, P C, and P N as profit, stores the result as profit data, and prints it. In the second pipeline, group by P C and P N creates nine independent price groups. Arrange by profit orders rows by profit within those groups, and slice of n takes each group’s final, highest-profit row. The output therefore has one selected nitrogen rate for every corn-price and nitrogen-price pair. Here sorting makes the selection easy to read, although the earlier Best N tab showed that filtering directly at the group maximum is cheaper when you only need maxima.
Now, we will calculate profit for all the rows in eval_parameters.
Now, we can identify the optimal N rate at each of the corn and nitrogen combinations:
So which should you use? The deciding factor is memory. Strategy one avoids materializing the full nitrogen-by-price grid, but each active worker can hold a working nitrogen grid and the result list retains every returned optimum row. Its memory use therefore grows with the worker count and the number of price combinations. Strategy two builds every combination at once, which is why it is fast, and also why it can exhaust your RAM. Once you run out of RAM your machine starts swapping to disk and performance collapses, far worse than the loop would ever have been. So the rule is simple. If the full grid fits comfortably in memory, vectorize. If it does not, break the job into chunks and loop over those. The callout shows you how to watch memory usage on both platforms.
For Strategy 1, reducing the number of workers can lower peak memory because fewer temporary nitrogen grids exist at the same time, while processing price combinations in batches limits how much returned data is retained. For Strategy 2, estimate the product of the dimension lengths before expanding the grid, since every added dimension multiplies the row count. To observe the consequence while code runs, the callout directs Mac users to Activity Monitor under Applications and Utilities. Windows users can press the Windows key and R, enter resmon, and inspect memory in Resource Monitor. Watch the memory pressure during a representative batch rather than waiting for the computer to become unresponsive.
Which strategy you should take depends on the size of your computer’s RAM.
Going over the RAM memory limit will suffocate your computer, which leads to a substantial loss in computing performance.
Vectorized version is more memory-hungry:
If you can fit the entire dataset in the RAM memory, then take Strategy 2. Otherwise, break up the entire task into pieces like Strategy 1.
Keep track of RAM memory usage
Mac users: go to Applications \rightarrow Utilities \rightarrow Activity Monitor
Windows users: press Windows Key + R \rightarrow type “resmon” into the search box
One more wrinkle, and it is the one that makes this realistic. Until now the yield response has been the same everywhere. In an actual field it is not. Soil varies, and here yield depends on electrical conductivity and slope as well as nitrogen, so the optimal rate differs from one part of the field to another. Look carefully at the vocabulary in the last line, because it now has three categories rather than two. Nitrogen is still the decision variable, the prices are still parameters, but slope and EC are attributes: fixed characteristics of a place that you cannot choose and cannot change.
In the displayed yield equation, EC divided by forty and one plus slope scale the twenty-five-times-log-N response term. The objective substitutes that location-specific yield into corn revenue and subtracts nitrogen cost, just as before. For any one plot and one price pair, slope, EC, P C, and P N are held fixed while you search across positive N values. The labels and units below the equation keep those roles separate: N is pounds per acre, P C is dollars per bushel, and P N is dollars per pound. The next tabs focus on preserving each plot’s legitimate attribute pair when constructing that search grid.
Suppose you have run a randomized nitrogen experiment for corn production on a field, collected data, and run a regression to find the following quantitative relationship between corn yield (bu/acre) and nitrogen rate (lb/acre):
\text{corn yield} = 120 + (EC/40) \times (1 + slope)\times 25 \times log(\text{nitrogen rate})
You are interested in finding the best nitrogen rates that maximize profit for different parts of the field at a given corn and nitrogen price combination.
Max_{N} P_C \times [120 + (EC/40) \times (1 + slope) * 25 \times log(\text{nitrogen rate})] - P_N \times N
N: nitrogen rate (lb/acre)slope: the slopeEC: electrical conductivityP_C: corn price ($/bu)P_N: nitrogen price ($/lb)Here, N is a decision variable, slope and EC are attributes, and P_C and P_N are parameters.
A two-plot field, kept deliberately tiny so you can see everything at once. Plot one has zero slope and an EC of forty. Plot two has a slope of nought point two and an EC of thirty. That is the entire dataset. The objective underneath is to find the best nitrogen rate for each plot at a given price combination. Keep those two pairs in your head for the next two tabs, because the whole point of what follows is which slope and EC values legitimately go together. Two plots means exactly two valid combinations, and any method that produces more than two has invented something.
The three vectors passed to data frame have matching length two, so their positions define the rows: the first plot ID travels with the first slope and first EC, and the second values travel together likewise. The left arrow stores those paired attributes as field data, while the surrounding parentheses display both rows. Autorun ensures the object exists for the two strategy tabs. The optimization target is one answer per plot ID, not one answer per independently imagined slope and EC value, which is why retaining that identifier and the row pairing matters.
Data
Consider a 2-plot field like below for the sake of tractability:
Objective
You want to find the optimal nitrogen rate for each plot for a given combination of corn and nitrogen prices.
The obvious approach is the one we have used all along: expand grid over everything, slope, EC and nitrogen rate. Run it and count. Expand grid gives you every slope crossed with every EC, so you get four combinations of soil properties. But only two of them exist. There is no part of this field with a slope of nought point two and an EC of forty, and none with zero slope and an EC of thirty. Those rows are fabrications. Here it merely wastes effort, but on a real field with many plots you would be computing optimal rates for soils that do not exist and might quietly average them in later.
The arguments slope equals field data dollar-sign slope and EC equals field data dollar-sign EC hand expand dot grid two separate vectors, so the function has no knowledge that their positions were paired by plot. It crosses both slope values with both EC values and then crosses those four attribute pairs with every value in the previously defined N seq. With three thousand nitrogen candidates, that produces twelve thousand rows instead of the six thousand rows required for two real plots. It also drops plot ID, making the fabricated combinations harder to detect. The parenthesized assignment stores and displays this oversized table as eval data one, and autorun lets you inspect the problem immediately.
You can expand on all the variables, nitrogen rate (decision variable), slope and EC (attributes), and corn and nitrogen prices (parameters):
slope-ec combinations of c(0, 30) and c(0.2, 40)The fix is to stop treating slope and EC as independent lists, because they are not. They arrive together, as rows of the field data. Expand grid dot df from the reshape package does exactly that. It crosses whole rows of one data frame with rows of another, rather than crossing individual columns. So the slope and EC pairings stay intact, each plot keeps its own soil properties, and you only get combinations that actually exist. That distinction, crossing rows rather than crossing columns, is worth remembering. It comes up whenever some of your variables are attributes that travel together instead of dimensions you are free to vary independently.
The double colon in reshape double colon expand dot grid dot d f calls the function directly from reshape without attaching the whole package. Field data supplies two intact rows containing plot ID, slope, and EC. Data frame of N seq turns the three-thousand-value nitrogen sequence into a one-column data frame, and the function crosses each real plot row with every nitrogen row. The automatic output therefore has six thousand rows, preserves plot ID, and never creates the two impossible soil pairings. Compared with Strategy 1, the computational saving is useful, but the more important benefit is that the data now represent cases that could actually occur.
Instead of applying expand.grid() on all the three vectors using expand.grid(), we can use expand.grid.df() from the reshape package as follows.
This generates all unique combinations of rows from the two data frames, without creating any observations that are not possible in reality.