Learning How to Use Shapefiles in R

Question

I am trying to learn more about "shapefiles" and how to plot them in R. Ideally, I would like to make a "leaflet map" (interactive or static) that looks something like this: https://i.cbc.ca/1.5587407.1591305080!/fileImage/httpImage/covid-19-cases-in-the-city-of-toronto-since-may-1.jpg

For this example, I found a publicly available dataset that contains information about Canadian Postal Codes - I created a new random variable that represents "number of hockey players" in each postal code:

# download folder from here : https://www.serviceobjects.com/blog/free-zip-code-and-postal-code-database-with-geocoordinates/ 
#extract canadian part (CanadianPostalCodes202204) - call it "canada"

canada = read.csv("CanadianPostalCodes202204.csv")

# extract Ontario
some_prov <- canada[which(canada$PROVINCE_ABBR == "ON"), ]

# create some random variable for this example

some_prov$number_of_hockey_players = as.integer(rnorm(nrow(some_prov), 87, 10))

   POSTAL_CODE         CITY PROVINCE_ABBR TIME_ZONE LATITUDE LONGITUDE number_of_hockey_players
11     L7C 2T6 CALEDON EAST            ON         5 43.81709 -79.81907                       88
17     M1K 1W6  SCARBOROUGH            ON         5 43.71187 -79.26526                       98
21     N5Y 4N3       LONDON            ON         5 43.02114 -81.22003                       72
22     N4N 1J3      HANOVER            ON         5 44.14945 -81.03273                       87
24     N2M 3R9    KITCHENER            ON         5 43.43482 -80.48112                       99
26     L1C 0E4  BOWMANVILLE            ON         5 43.92937 -78.68348                       90

I know how to make a general leaflet map for this data - this looks something like this:

library(leaflet)

some_prov %>% 
  leaflet() %>% 
  addTiles() %>% 
  addMarkers(clusterOption=markerClusterOptions())

I would now like to make a map similar to the COVID map above. Doing some research online, it seems that making such a map might require a "shapefile", or some information that tells the map where to put the desired boundaries. For example, I found the following website that might be useful (https://www.gis-blog.com/create-a-leaflet-map-using-r/ - I changed "USA" to "CANADA"):

#load leaflet package for R
library(leaflet)
library(raster)

#arbitrary data (code from above)
canada <- getData("GADM", country="canada", level=2)
canada$number_of_hockey_players <- as.integer(rnorm(n=nrow(canada), 150, 30))

#create a color palette to fill the polygons
pal <- colorQuantile("Greens", NULL, n = 5)

#create a pop up (onClick)
polygon_popup <- paste0("<strong>Name: </strong>", canada$NAME_1, "<br>",
                        "<strong>Indicator: </strong>", round(canada$number_of_hockey_players,2))

#create leaflet map
  map = leaflet() %>% 
      addProviderTiles("CartoDB.Positron") %>% 
      setView(-98.35, 39.7,
              zoom = 4) %>% 
      addPolygons(data = canada, 
                  fillColor= ~pal(number_of_hockey_players),
                  fillOpacity = 0.4, 
                  weight = 2, 
                  color = "white",
                  popup = polygon_popup)
map

But I am not sure if this is the best way to do things.

I think I found some publicly available shapefile over here : https://www12.statcan.gc.ca/census-recensement/2011/geo/bound-limit/bound-limit-2011-eng.cfm - but I am not sure how to correctly use these files to accomplish my goal. I am currently searching for other relevant shapefiles.

Can someone please show me how I can take this dataset that I created and make a map (interactive or static) similar to the COVID map above (with either postal code boundaries, city boundaries, etc.) ?

Thank you!