First load the data sets

mypath <- "C://Users//np83zg//OneDrive - Aalborg Universitet//Skrivebord//denmark"
setwd(mypath)
library(data.table)
mydata1 <- as.data.table(read.table("dat1.txt"))
mydata2 <- as.data.table(read.table("dat2.txt"))
mydata3 <- as.data.table(read.table("dat3.txt"))
mydata4 <- as.data.table(read.table("dat4.txt"))

Then calculate density by municipality as municipality workers per square kilometer

setkey(mydata1,workplace,year)
mydata1[,density:=count/size]
mydata_2003 <- mydata1[year==2003,]
mydata_2008 <- mydata1[year==2008,]
mydata_2016 <- mydata1[year==2016,]

mydata_2016[,comsize:=sum(size),by=.(commutingzone)]
mydata_2016[,comcount:=sum(count),by=.(commutingzone)]
mydata_2016[,com_density:=comcount/comsize]
mydata_2016[,com_log_density:=log(comcount/comsize)]
mydata_2016[,com_m_wage:=sum(size*m_wage)/comsize,by=.(commutingzone)]
mydata_2016[,com_s_uni_edu:=sum(size*s_uni_edu)/comsize,by=.(commutingzone)]
mydata_2016[,com_m_jobtenure:=sum(size*m_jobtenure)/comsize,by=.(commutingzone)]
mydata_2016[,com_m_exper:=sum(size*m_exper)/comsize,by=.(commutingzone)]

index <- match(unique(mydata_2016$commutingzone), mydata_2016$commutingzone)
index
##  [1]  1 42 44 46 47 48 55 56 58 59 60 61 70 71 72 74 80 84 86 87 88 90 91 93 98
comtab <- mydata_2016[index,]
comtab <- comtab[order(comtab$com_density),]
comtab <- comtab[,c("commutingzone","workplace","comcount","comsize","com_density","com_m_wage",
                         "com_m_exper","com_s_uni_edu","com_m_jobtenure")]
## Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1
## Reading layer `mun_shape' from data source `C:\Users\np83zg\OneDrive - Aalborg Universitet\Skrivebord\denmark\mun_shape.shp' using driver `ESRI Shapefile'
## Simple feature collection with 98 features and 4 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: 8.074458 ymin: 54.55906 xmax: 15.15735 ymax: 57.75238
## geographic CRS: WGS 84

Commutingzones

Now add some data from the dataset and plot these instead. First illustrate the commuting zones

names(mun_shape)[4] <- "workplace"
mydata_2008 <- merge(mun_shape,as.data.frame(mydata_2008),by="workplace")
mydata_2016 <- merge(mun_shape,as.data.frame(mydata_2016),by="workplace")
mydata_2008$log_density <- log(mydata_2008$density)
mydata_2016$log_density <- log(mydata_2016$density)
mydata_2016$dlog_density <- 100*(log(mydata_2016$density) - log(mydata_2008$density))
mydata_2016$dm_wage <- (mydata_2016$m_wage - mydata_2008$m_wage)/mydata_2008$m_wage
mydata_2016$ds_uni_edu <- (mydata_2016$s_uni_edu - mydata_2008$s_uni_edu)
#cols <- sapply(mydata_2016, is.numeric)
#cols <- names(cols)[cols]
#comtab <- mydata_2016[, lapply(.SD, mean), .SDcols = cols,by=.(commutingzone)]
tmap_mode("plot")
## tmap mode set to plotting
tm_shape(mydata_2016) + tm_polygons("commutingzone",palette="Accent",
                                    n=length(unique(mydata_2016$commutingzone))) +
tm_shape(mydata_2016) + tm_borders("white", lwd = 1.2)

# Commutingzones Descriptive

commutingzone workplace comcount comsize com_density com_m_wage com_m_exper com_s_uni_edu com_m_jobtenure
2 825 434 118.9 3.65 198.97 6.56 0.21 4.49
8 741 945 113.5 8.33 196.86 6.40 0.25 3.49
17 550 11235 1281.8 8.77 193.34 6.71 0.22 4.36
19 482 2933 288.8 10.16 189.92 7.11 0.20 3.97
13 665 6063 502.8 12.06 201.10 7.19 0.23 4.79
24 360 11932 881.9 13.53 192.46 6.79 0.21 3.70
6 787 15836 1068.6 14.82 195.44 7.51 0.21 4.66
11 760 22186 1469.6 15.10 206.67 7.61 0.19 4.08
20 492 1494 90.1 16.58 193.18 6.83 0.24 3.87
5 773 6453 366.5 17.61 192.19 6.37 0.18 4.25
23 376 16134 899.7 17.93 194.07 7.33 0.27 4.00
22 400 12651 588.1 21.51 188.86 7.47 0.22 4.00
4 860 21417 924.2 23.17 195.96 7.48 0.25 4.01
7 779 53272 2092.0 25.46 201.95 7.46 0.27 4.31
10 657 52465 2054.8 25.53 203.69 6.63 0.24 4.16
12 661 29160 1039.2 28.06 195.90 6.68 0.25 4.45
16 561 63282 2089.3 30.29 200.13 6.73 0.23 4.01
1 810 143957 4743.5 30.35 199.92 6.62 0.25 4.02
3 813 20660 650.3 31.77 199.38 6.92 0.19 4.14
18 540 46785 1436.5 32.57 202.38 7.40 0.25 4.29
14 530 103855 2670.2 38.89 208.34 6.71 0.26 4.03
15 410 106663 2668.9 39.97 200.28 6.27 0.23 3.86
21 420 140404 2801.0 50.13 197.42 6.34 0.26 3.94
9 706 275837 4657.2 59.23 201.15 6.23 0.28 3.88
25 101 968892 7394.4 131.03 205.88 6.70 0.27 3.58

Municipality Descriptives

Then make a table with municipality statistics

mun commutingzone count size density m_wage m_exper s_uni_edu m_jobtenure
96 Jammerbugt 1 9476 863.9 10.97 194.50 6.11 0.23 3.79
94 Rebild 1 7993 621.3 12.86 198.29 5.83 0.23 3.73
90 Brønderslev 1 9256 633.0 14.62 195.90 5.89 0.26 4.23
92 Vesthimmerlands 1 12256 769.9 15.92 197.01 6.68 0.21 4.25
95 Mariagerfjord 1 14948 718.2 20.81 202.93 7.42 0.23 4.20
97 Aalborg 1 90028 1137.2 79.17 207.25 7.31 0.32 3.96
93 Læsø 2 434 118.9 3.65 198.97 6.56 0.21 4.49
91 Frederikshavn 3 20660 650.3 31.77 199.38 6.92 0.19 4.14
98 Hjørring 4 21417 924.2 23.17 195.96 7.48 0.25 4.01
86 Morsø 5 6453 366.5 17.61 192.19 6.37 0.18 4.25
88 Thisted 6 15836 1068.6 14.82 195.44 7.51 0.21 4.66
87 Skive 7 16047 683.3 23.48 196.35 7.13 0.21 4.27
89 Viborg 7 37225 1408.7 26.43 204.66 7.61 0.30 4.32
80 Samsø 8 945 113.5 8.33 196.86 6.40 0.25 3.49
74 Syddjurs 9 9444 689.8 13.69 195.47 5.55 0.27 3.79
75 Norddjurs 9 10779 720.9 14.95 196.49 6.67 0.22 3.87
77 Odder 9 5093 223.6 22.78 191.66 5.62 0.28 3.89
76 Favrskov 9 13473 540.3 24.94 201.51 5.36 0.28 3.90
79 Silkeborg 9 28914 850.5 34.00 204.98 6.65 0.29 4.12
78 Randers 9 30518 747.6 40.82 197.97 6.43 0.27 3.91
81 Skanderborg 9 19276 416.7 46.26 206.40 6.25 0.28 4.09
82 Aarhus 9 158340 467.8 338.48 214.29 6.79 0.37 3.37
83 Ikast-Brande 10 17610 733.4 24.01 213.00 7.30 0.23 4.01
70 Herning 10 34855 1321.4 26.38 198.52 6.27 0.25 4.24
84 Ringkøbing-Skjern 11 22186 1469.6 15.10 206.67 7.61 0.19 4.08
73 Struer 12 6019 246.2 24.45 196.06 7.16 0.24 4.39
71 Holstebro 12 23141 793.0 29.18 195.85 6.53 0.25 4.47
72 Lemvig 13 6063 502.8 12.06 201.10 7.19 0.23 4.79
85 Hedensted 14 14952 551.3 27.12 203.98 7.63 0.21 4.09
58 Billund 14 14666 540.3 27.14 223.44 6.23 0.24 4.84
69 Vejle 14 41527 1058.8 39.22 204.89 6.77 0.28 3.64
67 Horsens 14 32710 519.8 62.93 204.28 6.08 0.26 3.94
64 Vejen 15 14184 813.5 17.44 198.75 6.53 0.20 3.87
57 Haderslev 15 17025 815.9 20.87 194.74 5.76 0.23 3.85
47 Middelfart 15 12312 298.8 41.20 200.45 5.10 0.25 3.70
68 Kolding 15 40763 607.1 67.14 206.44 7.06 0.26 3.93
66 Fredericia 15 22379 133.6 167.51 215.17 6.90 0.25 3.93
62 Fanø 16 541 54.6 9.91 183.61 6.00 0.25 4.13
63 Varde 16 15049 1240.0 12.14 193.51 6.65 0.21 3.97
61 Esbjerg 16 47692 794.7 60.01 211.58 6.89 0.26 4.07
60 Tønder 17 11235 1281.8 8.77 193.34 6.71 0.22 4.36
65 Aabenraa 18 21561 940.7 22.92 202.95 7.46 0.25 4.11
59 Sønderborg 18 25224 495.8 50.88 201.29 7.28 0.26 4.63
55 Langeland 19 2933 288.8 10.16 189.92 7.11 0.20 3.97
56 Ærø 20 1494 90.1 16.58 193.18 6.83 0.24 3.87
54 Nordfyns 21 6304 452.3 13.94 193.19 6.89 0.25 3.93
48 Assens 21 9932 511.3 19.42 196.29 5.85 0.25 4.19
49 Faaborg-Midtfyn 21 12818 633.7 20.23 195.17 5.73 0.25 3.79
51 Nyborg 21 8115 276.7 29.33 194.76 6.36 0.26 3.90
50 Kerteminde 21 6924 205.8 33.64 204.48 6.04 0.21 4.02
53 Svendborg 21 17249 415.5 41.51 198.44 7.12 0.31 3.91
52 Odense 21 79062 305.7 258.63 206.51 6.74 0.33 3.85
46 Bornholms 22 12651 588.1 21.51 188.86 7.47 0.22 4.00
44 Guldborgsund 23 16134 899.7 17.93 194.07 7.33 0.27 4.00
42 Lolland 24 11932 881.9 13.53 192.46 6.79 0.21 3.70
39 Stevns 25 3780 250.2 15.11 195.91 6.64 0.23 3.12
45 Vordingborg 25 11086 619.5 17.90 197.38 7.16 0.26 3.26
41 Lejre 25 4772 238.9 19.97 202.82 7.08 0.29 3.72
35 Faxe 25 8323 405.0 20.55 198.60 6.31 0.23 3.65
33 Odsherred 25 7480 354.0 21.13 196.20 6.18 0.23 3.59
36 Kalundborg 25 14697 575.2 25.55 221.13 7.52 0.24 4.16
40 Sorø 25 8567 308.4 27.78 206.90 7.16 0.28 3.58
32 Gribskov 25 8122 279.5 29.06 196.04 6.78 0.25 3.19
43 Næstved 25 22477 676.3 33.24 198.63 6.87 0.28 3.41
34 Holbæk 25 19906 577.3 34.48 200.62 6.62 0.29 3.44
37 Ringsted 25 12262 294.6 41.62 207.34 6.69 0.22 3.86
38 Slagelse 25 25720 567.9 45.29 199.37 6.55 0.27 3.38
26 Frederikssund 25 11474 247.0 46.45 202.35 6.69 0.25 3.49
29 Halsnæs 25 5669 121.8 46.54 199.36 6.89 0.25 3.80
25 Egedal 25 7992 125.8 63.53 209.51 6.26 0.34 3.36
20 Fredensborg 25 8685 112.1 77.48 213.47 5.84 0.35 3.54
28 Køge 25 21795 256.5 84.97 207.59 6.28 0.25 3.94
31 Solrød 25 3987 40.1 99.43 198.35 6.26 0.26 3.46
5 Dragør 25 2056 18.3 112.35 195.17 5.02 0.27 3.19
22 Hillerød 25 24293 214.9 113.04 223.05 6.85 0.39 3.94
21 Helsingør 25 15785 118.9 132.76 200.26 6.35 0.28 3.54
30 Roskilde 25 32770 211.9 154.65 213.52 6.73 0.35 3.76
18 Furesø 25 9135 56.8 160.83 212.49 6.04 0.31 3.74
19 Allerød 25 11059 67.5 163.84 231.94 6.39 0.35 3.88
27 Greve 25 12852 60.4 212.78 203.35 6.00 0.20 3.33
23 Hørsholm 25 6707 31.3 214.28 219.88 6.14 0.36 3.23
15 Ishøj 25 6172 26.4 233.79 207.06 6.05 0.21 3.93
24 Rudersdal 25 19570 73.4 266.62 235.99 7.04 0.40 3.68
16 TÃ¥rnby 25 21389 66.0 324.08 219.49 7.29 0.17 4.25
17 Vallensbæk 25 3656 9.5 384.84 210.18 5.75 0.30 3.17
12 Høje-Taastrup 25 30249 78.4 385.83 226.81 6.60 0.25 3.60
13 Lyngby-Taarbæk 25 25929 38.8 668.27 235.12 6.28 0.45 3.70
10 Albertslund 25 16593 23.2 715.22 219.44 6.16 0.21 3.80
4 Brøndby 25 19402 21.0 923.90 230.92 5.67 0.27 3.89
11 Hvidovre 25 24025 23.0 1044.57 223.69 6.25 0.30 3.76
14 Rødovre 25 12690 12.1 1048.76 202.88 5.89 0.21 3.44
3 Ballerup 25 36325 33.8 1074.70 253.08 6.06 0.35 4.16
6 Gentofte 25 28908 25.6 1129.22 237.70 5.54 0.47 3.44
8 Glostrup 25 18117 13.3 1362.18 234.10 6.32 0.35 3.91
7 Gladsaxe 25 35683 24.9 1433.05 252.42 6.16 0.42 3.29
9 Herlev 25 18177 12.1 1502.23 228.30 6.13 0.38 3.81
2 Frederiksberg 25 31094 8.1 3838.77 215.61 5.31 0.43 3.38
1 Københavns 25 299462 74.7 4008.86 227.89 5.16 0.45 3.08 
tmap_mode("plot")
## tmap mode set to plotting
tm_shape(mydata_2016) + tm_polygons("m_wage",palette="GnBu") +
tm_shape(mydata_2016) + tm_borders("white", lwd = 1.2)

tmap_mode("plot")
## tmap mode set to plotting
tm_shape(mydata_2016) + tm_polygons("log_density",palette="YlOrBr") +
tm_shape(mydata_2016) + tm_borders("white", lwd = 1.2)

tmap_mode("plot")
## tmap mode set to plotting
tm_shape(mydata_2016) + tm_polygons("s_uni_edu",palette="Accent") +
tm_shape(mydata_2016) + tm_borders("white", lwd = 1.2)

tmap_mode("plot")
## tmap mode set to plotting
tm_shape(mydata_2016) + tm_polygons("m_exper",palette="GnBu") +
tm_shape(mydata_2016) + tm_borders("white", lwd = 1.2)

tmap_mode("plot")
## tmap mode set to plotting
tm_shape(mydata_2016) + tm_polygons("m_jobtenure",palette="GnBu") +
tm_shape(mydata_2016) + tm_borders("white", lwd = 1.2)

tmap_mode("plot")
## tmap mode set to plotting
tm_shape(mydata_2016) + tm_polygons("dlog_density",palette="YlOrBr") +
tm_shape(mydata_2016) + tm_borders("white", lwd = 1.2)

tmap_mode("plot")
## tmap mode set to plotting
tm_shape(mydata_2016) + tm_polygons("dm_wage",palette="GnBu") +
tm_shape(mydata_2016) + tm_borders("white", lwd = 1.2)

tmap_mode("plot")
## tmap mode set to plotting
tm_shape(mydata_2016) + tm_polygons("ds_uni_edu",palette="GnBu") +
tm_shape(mydata_2016) + tm_borders("white", lwd = 1.2)