How do I calculate the Terrain Ruggedness Index for each county in the United States in R?

Just because something is published does not mean that it is necessarly correct. In this case aggregating the TRI to a county is certainly incorrect. The distributional qualities of the metric, in relation to inference, become meaningless. Given the linked journal, bad dogs! You are functionally taking the mean of a derivative metric that represents localized mean deviation.

I would highly recommend reading up on MAUP, perhaps starting with Cressie's "Change of support and the modifiable areal unit problem" and ecological fallacy in spatial data by reading Wakefield's "Spatial Aggregation and the Ecological Fallacy".

Since the basic idea here is to identify topographic variability within an experimental unit to indicate "ruggedness", one could address the underlying distributions directly. Since highly relieved areas would also be expected to exhibit highly skewed, standard Gaussian moments may not be adequate. You can step out into non parametric statistics such as Median Absolute Deviation from Median (MAD).

Here is an example of what I am getting at and some potential solutions.

Add libraries and data.

library(raster)
library(spatialEco)
library(elevatr)
library(USAboundaries)

counties <- as(us_counties(map_date = "1930-01-01", 
              resolution = "high", states = c("CA")),
              "Spatial")
elev <- get_elev_raster(counties, z=5)

First, let's calculate the pixel-level TRI and calculate the mean for each county. You can see that the variability is not correctly represented, at least not visually.

r.tri <- spatialEco::tri(elev) 
counties@data <- data.frame(counties@data, r.tri = extract(r.tri, 
                            counties, fun=mean))
spplot(counties, "r.tri") 

Now we can calculate the MAD by passing the function directly to raster::extract.

counties@data <- data.frame(counties@data, tri = extract(elev, 
                            counties, fun=tri))
spplot(counties, "rough")

We can also write a global approximation of TRI using the median and the deviation value. This actually looks fairly reasonable and is comparable to MAD. Although, it did pick up Frenso county as very high ruggedness (which spans the southern Sierra's) whereas MAD did not.

tri <- function(x, ...) {
  x <- x[!is.na(x)]
  return( sqrt(sum(((median(x) - x)^2))) )
}

counties@data <- data.frame(counties@data, tri = extract(elev, 
                            counties, fun=tri))
spplot(counties, "tri")