Analysis for Ram

Load Library & Data

library(conover.test)
library(nloptr)
library(rgdal)

## Loading required package: sp

## Please note that rgdal will be retired during 2023,
## plan transition to sf/stars/terra functions using GDAL and PROJ
## at your earliest convenience.
## See https://r-spatial.org/r/2022/04/12/evolution.html and https://github.com/r-spatial/evolution
## rgdal: version: 1.6-5, (SVN revision 1199)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.5.2, released 2022/09/02
## Path to GDAL shared files: C:/Users/lfult/AppData/Local/R/win-library/4.2/rgdal/gdal
## GDAL binary built with GEOS: TRUE 
## Loaded PROJ runtime: Rel. 8.2.1, January 1st, 2022, [PJ_VERSION: 821]
## Path to PROJ shared files: C:/Users/lfult/AppData/Local/R/win-library/4.2/rgdal/proj
## PROJ CDN enabled: FALSE
## Linking to sp version:1.6-0
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading sp or rgdal.

library(psych)        
library(raster)      
library(RColorBrewer)
library(rgdal)       
library(rgeos) 

## rgeos version: 0.6-2, (SVN revision 693)
##  GEOS runtime version: 3.9.3-CAPI-1.14.3 
##  Please note that rgeos will be retired during 2023,
## plan transition to sf functions using GEOS at your earliest convenience.
##  GEOS using OverlayNG
##  Linking to sp version: 1.6-0 
##  Polygon checking: TRUE

library(spatialEco)

## 
## Attaching package: 'spatialEco'

## The following object is masked from 'package:raster':
## 
##     shift

library(tmap)        
library(tmaptools) 
library(kableExtra)
library(knitr)
library(leaflet)  
library(leaflet.extras)
library(totalcensus)

mydata=read.csv('c:/users/lfult/documents/ram/statedata.csv', stringsAsFactors = T)

Load Shapefile

setwd('C:/Users/lfult/Documents/Ram')
myshape=shapefile("cb_2018_us_state_500k.shp") 

Describe

options(scipen=99999)

describe(mydata[, c(1185:ncol(mydata)-1)])%>%kbl()%>%kable_classic(html_font='Cambria')

	vars	n	mean	sd	median	trimmed	mad	min	max	range	skew	kurtosis	se
Total	1	51	21554.1960784	23740.5854806	14194.0000000	16654.4634146	15266.3322000	910.0000000	101344.0000000	100434.0000000	1.8569017	2.9901988	3324.3469032
var	2	51	3933.3979804	6888.1558946	639.5214622	2392.7317781	886.4657914	7.4298818	33823.2045200	33815.7746382	2.3236794	5.7272700	964.5347515
xbar	3	51	18.2353605	20.0850977	12.0084602	14.0900706	12.9156787	0.7698816	85.7394247	84.9695431	1.8569017	2.9901988	2.8124762
xbar_biased	4	51	33.5989352	31.9428496	21.7262295	27.8628413	20.9738680	2.8980892	118.4034091	115.5053199	1.3432795	0.7485017	4.4728936
n	5	51	1182.0000000	0.0000000	1182.0000000	1182.0000000	0.0000000	1182.0000000	1182.0000000	0.0000000	NaN	NaN	0.0000000
n_biased	6	51	600.9803922	177.8128781	610.0000000	610.5609756	163.0860000	170.0000000	911.0000000	741.0000000	-0.4435440	-0.3596224	24.8987832
Delta	7	51	0.9706574	0.0314283	0.9796617	0.9759335	0.0202657	0.8414961	0.9987270	0.1572308	-1.8273924	3.8767779	0.0044008
Theta	8	51	0.2102130	0.1293119	0.2015234	0.2013180	0.1528013	0.0195444	0.5984396	0.5788953	0.6574610	0.1616377	0.0181073
Beta	9	51	0.2172816	0.1346885	0.2106093	0.2074573	0.1587638	0.0195693	0.6484961	0.6289268	0.7250731	0.4843069	0.0188602
BiasedTheta	10	51	1.5362662	0.4798542	1.4145480	1.4892150	0.4488485	0.8242908	2.9465289	2.1222381	0.7942653	-0.0377244	0.0671930
BiasedBeta	11	51	1.9722063	0.5708166	1.8527853	1.9465832	0.5530646	0.9057147	3.4527112	2.5469965	0.4145925	-0.5037366	0.0799303
BiasedDelta	12	51	0.7824102	0.0935657	0.8007953	0.7895562	0.0937878	0.5263758	0.9115819	0.3852061	-0.7047926	-0.1429572	0.0131018
sum.xi..1.	13	51	116.4208796	62.8603054	98.1464465	111.6726869	68.7487665	14.5341008	289.5494289	275.0153280	0.6118526	-0.3970701	8.8022034
Population2021	14	51	6510422.6274510	7394300.0767053	4506589.0000000	5005480.9024390	4089437.7888000	579483.0000000	39142991.0000000	38563508.0000000	2.4782690	6.9828900	1035409.1132737
Rate_100K	15	51	320.7400909	86.5440936	335.4873125	324.8457208	91.2386147	121.4798149	454.6724377	333.1926228	-0.4016157	-0.8298810	12.1185971
Region*	16	51	2.6078431	1.1149607	3.0000000	2.6341463	1.4826000	1.0000000	4.0000000	3.0000000	-0.2234952	-1.3422133	0.1561257
Division*	17	51	5.2745098	2.4501300	5.0000000	5.3658537	2.9652000	1.0000000	9.0000000	8.0000000	-0.2482803	-1.1128120	0.3430868

Region

library(dplyr)

## 
## Attaching package: 'dplyr'

## The following object is masked from 'package:kableExtra':
## 
##     group_rows

## The following object is masked from 'package:spatialEco':
## 
##     combine

## The following objects are masked from 'package:rgeos':
## 
##     intersect, setdiff, symdiff, union

## The following objects are masked from 'package:raster':
## 
##     intersect, select, union

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library(magrittr)

## 
## Attaching package: 'magrittr'

## The following object is masked from 'package:raster':
## 
##     extract

tempdata1=subset(mydata, select=c(State, Total, Rate_100K, Region, n,   n_biased,   xbar, xbar_biased,Delta, Theta, Beta, BiasedDelta, BiasedBeta, BiasedTheta )) 
tempdata1=tempdata1[order(tempdata1$Region),]
tempdata1%>%kbl()%>%kable_classic(html_font='Cambria')

	State	Total	Rate_100K	Region	n	n_biased	xbar	xbar_biased	Delta	Theta	Beta	BiasedDelta	BiasedBeta	BiasedTheta
13	IA	10538	329.5505	Midwest	1182	487	8.9153976	21.638604	0.9874450	0.1119330	0.1133562	0.8046602	1.4089343	1.1337134
15	IL	39381	310.4173	Midwest	1182	784	33.3172589	51.036990	0.9900199	0.3325094	0.3358613	0.8533957	3.4527112	2.9465289
16	IN	25959	380.9918	Midwest	1182	710	21.9619289	36.563380	0.9940843	0.1299204	0.1306935	0.8357833	2.5060841	2.0945433
17	KS	10167	346.0609	Midwest	1182	390	8.6015228	26.330769	0.9896722	0.0888348	0.0897619	0.8236600	1.1732698	0.9663754
23	MI	42613	424.5378	Midwest	1182	572	36.0516074	74.512238	0.9924816	0.2710492	0.2731025	0.8847195	2.2996835	2.0345749
24	MN	12806	224.2154	Midwest	1182	674	10.8341794	19.000000	0.9653617	0.3752776	0.3887430	0.7759486	2.6069990	2.0228972
25	MO	20708	335.6336	Midwest	1182	598	17.5194585	35.872910	0.9981831	0.0318307	0.0318886	0.8477168	1.0350956	0.8774680
29	ND	2232	286.9138	Midwest	1182	399	1.8883249	5.596491	0.9452482	0.1033893	0.1093779	0.6488844	1.8188668	1.1802342
30	NE	4827	245.8298	Midwest	1182	500	4.0837563	9.916000	0.9843328	0.0639811	0.0649995	0.7250597	1.5487860	1.1229623
36	OH	42073	357.6316	Midwest	1182	500	35.5947546	85.538000	0.9973082	0.0958152	0.0960738	0.8927804	2.1798665	1.9461420
42	SD	3214	358.6397	Midwest	1182	348	2.7191201	9.255747	0.9536443	0.1260468	0.1321738	0.7188339	1.4797483	1.0636933
49	WI	16486	280.3693	Midwest	1182	911	13.9475465	18.267837	0.9767236	0.3246482	0.3323850	0.7911921	1.3044559	1.0320752
7	CT	11034	304.5244	Northeast	1182	483	9.3350254	22.913044	0.9796617	0.1898589	0.1938004	0.8099953	1.3726053	1.1118038
20	MA	21035	300.9432	Northeast	1182	678	17.7961083	36.907080	0.9986617	0.0238158	0.0238477	0.8361503	2.5462602	2.1290562
22	ME	2989	217.0286	Northeast	1182	578	2.5287648	5.178201	0.9239784	0.1922407	0.2080575	0.6407807	1.7464087	1.1190650
31	NH	2972	214.1974	Northeast	1182	545	2.5143824	5.794495	0.9786677	0.0536376	0.0548068	0.6524712	1.8527853	1.2088891
32	NJ	35774	385.9964	Northeast	1182	730	30.2656514	49.065753	0.9968909	0.0941001	0.0943936	0.8640544	1.5784414	1.3638592
35	NY	77403	389.7924	Northeast	1182	811	65.4847716	97.672010	0.9969226	0.2015234	0.2021455	0.8975778	2.6535755	2.3817906
39	PA	50860	390.8682	Northeast	1182	755	43.0287648	67.365563	0.9904586	0.4105562	0.4145112	0.8818616	1.7689515	1.5599704
40	RI	3897	355.2464	Northeast	1182	407	3.2969543	9.626536	0.9730977	0.0886958	0.0911479	0.7229646	1.4927161	1.0791808
47	VT	910	140.6552	Northeast	1182	314	0.7698816	2.898089	0.9186838	0.0626039	0.0681452	0.5263758	2.8425970	1.4962742
2	AL	21133	418.4880	South	1182	655	17.8790186	32.271756	0.9863449	0.2441391	0.2475190	0.8350520	1.5713197	1.3121336
3	AR	13062	431.3565	South	1182	779	11.0507614	16.813864	0.9674090	0.3601549	0.3722881	0.7734325	1.9545547	1.5117162
8	DC	1392	208.1368	South	1182	412	1.1776650	3.385922	0.8414961	0.1866644	0.2218245	0.5559740	2.6116877	1.4520305
9	DE	3371	335.4873	South	1182	513	2.8519459	6.573099	0.9432715	0.1617867	0.1715166	0.6686915	1.8519161	1.2383605
10	FL	87141	399.2153	South	1182	819	73.7233502	110.829060	0.9978203	0.1606924	0.1610434	0.9100999	0.9057147	0.8242908
11	GA	42348	392.5462	South	1182	715	35.8274112	65.492308	0.9973703	0.0942139	0.0944623	0.8800465	1.8070296	1.5902702
18	KY	18094	401.5010	South	1182	652	15.3079526	27.751534	0.9846860	0.2344263	0.2380722	0.8169234	2.1508300	1.7570633
19	LA	18136	391.9519	South	1182	695	15.3434856	26.096403	0.9764683	0.3610585	0.3697596	0.8007953	3.0866050	2.4717390
21	MD	15578	252.2912	South	1182	802	13.1793570	19.709476	0.9833961	0.2188291	0.2225239	0.7829888	2.3685813	1.8545726
26	MS	13097	444.0284	South	1182	610	11.0803722	21.726229	0.9720622	0.3095608	0.3184578	0.7920624	2.3919640	1.8945846
28	NC	28945	273.9477	South	1182	615	24.4881557	47.292683	0.9918998	0.1983587	0.1999786	0.8520415	2.8799117	2.4538043
37	OK	18147	454.6724	South	1182	176	15.3527919	118.403409	0.9987270	0.0195444	0.0195693	0.9115819	1.2884707	1.1745465
41	SC	17869	344.0802	South	1182	599	15.1175973	29.859766	0.9812612	0.2832849	0.2886947	0.8234405	2.1114107	1.7386212
43	TN	28113	403.4383	South	1182	703	23.7842640	40.039829	0.9964745	0.0838517	0.0841483	0.8504132	1.5967276	1.3578783
44	TX	92159	311.7813	South	1182	840	77.9686971	109.800000	0.9933307	0.5199962	0.5234874	0.9062903	1.8558964	1.6819809
46	VA	23718	273.9633	South	1182	825	20.0659898	28.823030	0.9854391	0.2921792	0.2964965	0.8159022	2.5390911	2.0716500
50	WV	8083	452.6957	South	1182	656	6.8384095	12.748476	0.9693607	0.2095242	0.2161468	0.7502743	1.6512403	1.2388832
1	AK	1441	196.2729	West	1182	298	1.2191201	4.963087	0.9512042	0.0594879	0.0625396	0.6420678	1.5153181	0.9729369
4	AZ	29852	410.9085	West	1182	564	25.2554991	53.510638	0.9917298	0.2088675	0.2106093	0.8678458	1.8224349	1.5815926
5	CA	101344	258.9071	West	1182	909	85.7394247	112.929593	0.9957943	0.3605940	0.3621170	0.9099696	1.1710774	1.0656448
6	CO	14194	244.2484	West	1182	773	12.0084602	18.447607	0.9787091	0.2556709	0.2612328	0.7872528	1.6320888	1.2848665
12	HI	1758	121.4798	West	1182	346	1.4873096	5.109827	0.9004206	0.1481055	0.1644848	0.6303984	2.0467506	1.2902682
14	ID	5463	286.8750	West	1182	553	4.6218274	9.951175	0.9433108	0.2620079	0.2777535	0.7163596	1.9746340	1.4145480
27	MT	3705	334.9222	West	1182	529	3.1345178	7.017013	0.9246664	0.2361346	0.2553728	0.6762415	1.8791420	1.2707539
33	NM	9164	432.9428	West	1182	828	7.7529611	11.084541	0.9228115	0.5984396	0.6484961	0.7153445	2.6530290	1.8978296
34	NV	11976	380.6252	West	1182	591	10.1319797	20.282572	0.9747607	0.2557237	0.2623451	0.7895741	2.0937263	1.6531521
38	OR	8726	205.0137	West	1182	652	7.3824027	13.384969	0.9587416	0.3045865	0.3176940	0.7531805	1.7721549	1.3347525
45	UT	5341	159.9527	West	1182	586	4.5186125	9.138225	0.9160813	0.3791961	0.4139328	0.6979129	2.3571976	1.6451186
48	WA	16013	206.8664	West	1182	611	13.5473773	26.836334	0.9788164	0.2869824	0.2931933	0.8051148	2.9371596	2.3647505
51	WY	2023	349.1043	West	1182	170	1.7115059	12.323529	0.9681369	0.0545339	0.0563287	0.7507831	1.4360169	1.0781371

Kruskal Wallis, ~Region

mykw=aov(BiasedDelta~Region, data=mydata)
mykw

## Call:
##    aov(formula = BiasedDelta ~ Region, data = mydata)
## 
## Terms:
##                    Region Residuals
## Sum of Squares  0.0334837 0.4042435
## Deg. of Freedom         3        47
## 
## Residual standard error: 0.09274118
## Estimated effects may be unbalanced

pairwise.wilcox.test(mydata$BiasedDelta, mydata$Region, p.adjust.method = "bonferroni")

## 
##  Pairwise comparisons using Wilcoxon rank sum exact test 
## 
## data:  mydata$BiasedDelta and mydata$Region 
## 
##           Midwest Northeast South
## Northeast 1.00    -         -    
## South     1.00    1.00      -    
## West      0.52    1.00      0.23 
## 
## P value adjustment method: bonferroni

conover.test(mydata$BiasedDelta, mydata$Region)

##   Kruskal-Wallis rank sum test
## 
## data: x and group
## Kruskal-Wallis chi-squared = 4.4768, df = 3, p-value = 0.21
## 
## 
##                            Comparison of x by group                            
##                                 (No adjustment)                                
## Col Mean-|
## Row Mean |    Midwest   Northeas      South
## ---------+---------------------------------
## Northeas |   0.564036
##          |     0.2877
##          |
##    South |  -0.319019  -0.895125
##          |     0.3756     0.1876
##          |
##     West |   1.581509   0.886457   2.044827
##          |     0.0602     0.1899    0.0232*
## 
## alpha = 0.05
## Reject Ho if p <= alpha/2

boxplot(BiasedDelta~Region, data=mydata, notch=F, col=c('mistyrose', 'lightsteelblue','palegreen', 'plum'), horizontal=T, xlab='Delta Hat Biased', ylab='', cex.axis=0.6)

mykw1=kruskal.test(Rate_100K~Region, data=mydata)
mykw1

## 
##  Kruskal-Wallis rank sum test
## 
## data:  Rate_100K by Region
## Kruskal-Wallis chi-squared = 8.378, df = 3, p-value = 0.03881

pairwise.wilcox.test(mydata$Rate_100K, mydata$Region, p.adjust.method = "bonferroni")

## 
##  Pairwise comparisons using Wilcoxon rank sum exact test 
## 
## data:  mydata$Rate_100K and mydata$Region 
## 
##           Midwest Northeast South
## Northeast 1.000   -         -    
## South     0.702   0.233     -    
## West      1.000   1.000     0.081
## 
## P value adjustment method: bonferroni

conover.test(mydata$Rate_100K, mydata$Region)

##   Kruskal-Wallis rank sum test
## 
## data: x and group
## Kruskal-Wallis chi-squared = 8.378, df = 3, p-value = 0.04
## 
## 
##                            Comparison of x by group                            
##                                 (No adjustment)                                
## Col Mean-|
## Row Mean |    Midwest   Northeas      South
## ---------+---------------------------------
## Northeas |   0.495317
##          |     0.3113
##          |
##    South |  -1.578028  -1.973139
##          |     0.0606     0.0272
##          |
##     West |   1.215579   0.618516   2.935626
##          |     0.1151     0.2696    0.0026*
## 
## alpha = 0.05
## Reject Ho if p <= alpha/2

boxplot(Rate_100K~Region, data=mydata, notch=F, col=c('mistyrose', 'lightsteelblue', 'palegreen', 'plum'), horizontal=T, xlab='Deaths / 100K', ylab='', cex.axis=0.6)

Kruskal Wallis ~ Division

mykw=kruskal.test(BiasedDelta~Division, data=mydata)
mykw

## 
##  Kruskal-Wallis rank sum test
## 
## data:  BiasedDelta by Division
## Kruskal-Wallis chi-squared = 16.784, df = 8, p-value = 0.03244

pairwise.wilcox.test(mydata$BiasedDelta, mydata$Division, p.adjust.method = "bonferroni")

## 
##  Pairwise comparisons using Wilcoxon rank sum exact test 
## 
## data:  mydata$BiasedDelta and mydata$Division 
## 
##                  EastNorthCentral EastSouthCentral MiddleAtlantic Mountain
## EastSouthCentral 1.00             -                -              -       
## MiddleAtlantic   1.00             1.00             -              -       
## Mountain         0.39             1.00             0.87           -       
## NewEngland       1.00             1.00             0.86           1.00    
## Pacific          1.00             1.00             1.00           1.00    
## SouthAtlantic    1.00             1.00             1.00           1.00    
## WestNorthCentral 1.00             1.00             0.60           1.00    
## WestSouthCentral 1.00             1.00             1.00           1.00    
##                  NewEngland Pacific SouthAtlantic WestNorthCentral
## EastSouthCentral -          -       -             -               
## MiddleAtlantic   -          -       -             -               
## Mountain         -          -       -             -               
## NewEngland       -          -       -             -               
## Pacific          1.00       -       -             -               
## SouthAtlantic    1.00       1.00    -             -               
## WestNorthCentral 1.00       1.00    1.00          -               
## WestSouthCentral 1.00       1.00    1.00          1.00            
## 
## P value adjustment method: bonferroni

boxplot(BiasedDelta~Division, data=mydata, notch=F, col=rainbow(9), horizontal=T, xlab='Delta Hat Biased', ylab='', cex.axis=0.6)

mykw1=kruskal.test(Rate_100K~Division, data=mydata)
mykw1

## 
##  Kruskal-Wallis rank sum test
## 
## data:  Rate_100K by Division
## Kruskal-Wallis chi-squared = 24.75, df = 8, p-value = 0.001714

pairwise.wilcox.test(mydata$Rate_100K, mydata$Division, p.adjust.method = "bonferroni")

## 
##  Pairwise comparisons using Wilcoxon rank sum exact test 
## 
## data:  mydata$Rate_100K and mydata$Division 
## 
##                  EastNorthCentral EastSouthCentral MiddleAtlantic Mountain
## EastSouthCentral 1.00             -                -              -       
## MiddleAtlantic   1.00             1.00             -              -       
## Mountain         1.00             1.00             1.00           -       
## NewEngland       1.00             0.34             0.86           1.00    
## Pacific          0.29             0.57             1.00           1.00    
## SouthAtlantic    1.00             1.00             1.00           1.00    
## WestNorthCentral 1.00             0.22             0.60           1.00    
## WestSouthCentral 1.00             1.00             1.00           1.00    
##                  NewEngland Pacific SouthAtlantic WestNorthCentral
## EastSouthCentral -          -       -             -               
## MiddleAtlantic   -          -       -             -               
## Mountain         -          -       -             -               
## NewEngland       -          -       -             -               
## Pacific          1.00       -       -             -               
## SouthAtlantic    1.00       0.14    -             -               
## WestNorthCentral 1.00       0.36    1.00          -               
## WestSouthCentral 0.69       0.57    1.00          1.00            
## 
## P value adjustment method: bonferroni

boxplot(Rate_100K~Division, data=mydata, notch=F, col=rainbow(9), horizontal=T, xlab='Deaths / 100K', ylab='', cex.axis=0.6)

Merge Shape and Flat File

mydata$M=mydata$State
myshape@data$M=as.factor(myshape@data$STUSPS)
states=merge(myshape, mydata, by="M",type="right")
states=na.omit(states)

Plot Deaths, All

qpal<-colorBin("Reds",c(910, 101344),bins=6,pretty=TRUE, alpha=.1)
qpal2<-colorBin("Blues",c(121,455),bins=6,pretty=TRUE, alpha=.1)#all deaths
qpal3<-colorBin("Greens",c(0.53,0.91),bins=6,pretty=TRUE, alpha=.1)#N=137
qpal4<-colorBin("Oranges",c(0,11),bins=4,pretty=TRUE, alpha=.1)#N=25


leaf=leaflet() %>%
  addTiles(group = "OSM (default)") %>%
  addMapPane("borders", zIndex = 410) %>%
  
  #Base Diagrams
  addPolylines(data = states,color = "black",
               opacity = 1, weight = 1, group="Borders", options = pathOptions(pane="borders"))%>%
  fitBounds(-124.8, -66.9, 24.4,49.4) %>% setView(-98.6, 39.83, zoom = 4)%>%

  addPolygons(data=states,stroke = FALSE,fillOpacity = 1, smoothFactor = 0.2, 
              color=~qpal(states@data$Total),
              popup = paste("State: ", states@data$M, "<br>", 
                            "Deaths: ",states@data$Total), 
              group="Deaths") %>%
  
    addPolygons(data=states,stroke = FALSE,fillOpacity = 1, smoothFactor = 0.2, 
              color=~qpal2(states@data$Rate_100K),
              popup = paste("State: ", states@data$M, "<br>", 
                            "Death Rate: ",states@data$Rate_100K), 
              group="Death Rate") %>%
  
    addPolygons(data=states,stroke = FALSE,fillOpacity = 1, smoothFactor = 0.2, 
              color=~qpal3(states@data$BiasedDelta),
              popup = paste("State: ", states@data$M, "<br>", 
                            "Dispersion (biased): ",states@data$BiasedDelta), 
              group="Dispersion") %>%
  
  addLegend(data=states, 
            "topleft", opacity=1, pal = qpal, 
            values = ~Total,
            title = "Deaths")%>%
  
    
  addLegend(data=states, 
            "topright", opacity=1, pal = qpal2, 
            values = ~Rate_100K,
            title = "Death Rate/100K")%>%
  
    addLegend(data=states, 
            "bottomleft", opacity=1, pal = qpal3, 
            values = ~BiasedDelta,
            title = "Dispersion")%>%
  
  addLayersControl(
    baseGroups = c("Deaths", "Death Rate", "Dispersion"),
    overlayGroups = c("Borders"), options = layersControlOptions(collapsed = TRUE)
  )
  
leaf