
```{r setup}
dataDir = file.path('..', 'data')
dir.create(dataDir, showWarnings=FALSE)
```

```{r boundaries}
SyearBnd = c(2003, 2006, 2010,2013,2016)
Syear = 2014:2018

bndUrl = 'https://ec.europa.eu/eurostat/cache/GISCO/distribution/v2/nuts/download/#nuts16'

bndDetailedUrl = paste0(
	'https://ec.europa.eu/eurostat/cache/GISCO/distribution/v2/nuts/download/ref-nuts-',SyearBnd,'-20m.shp.zip')
names(bndDetailedUrl) = SyearBnd

bndFiles = mapply(
	Pmisc::downloadIfOld,
	url=bndDetailedUrl, 
		MoreArgs=list(path=dataDir))


bndFilesNuts2 = grep("RG.*3035.*LEVL_3", as.vector(bndFiles), value=TRUE)
bndFiles3 = mapply(unzip,
	zipfile= bndFilesNuts2,
	MoreArgs = list(exdir = dataDir))
bndFilesNuts4 = grep("[.]shp$", bndFiles3, value=TRUE)
nuts3 = mapply(
	raster::shapefile,
	bndFilesNuts4)
names(nuts3) = gsub(".*[_]20M[_]|[_]3035.*", "",
	names(nuts3))
library(raster)
toKeep = extent(2308241, 7600000, 180000, 6000000)
nuts = mapply(
	function(xx) {
#	xx2= rgeos::gSimplify(xx, tol=0.0001, topologyPreserve=TRUE)
	xx2 = SpatialPoints(xx)
	xx3 = xx[!is.na(over(xx2, as(toKeep, 'SpatialPolygons'))), ]
	xx3[xx3$CNTR_CODE != 'IS', ]
	},
	nuts3
	)
nuts[['2013']]$id2016 = nuts[['2016']]@data[
	match(nuts[['2013']]$NUTS_NAME,
	nuts[['2016']]$NUTS_NAME), 'FID']

stuff = over(SpatialPoints(nuts[['2013']], proj4string = nuts[['2013']]@proj4string),
	nuts[['2010']])

nuts[['2013']]$id2010 = stuff$FID

# nuts 2 for germany
germany2 = grep("RG.*2010_3035.*LEVL_2", as.vector(bndFiles), value=TRUE)
germany2 = unzip(germany2, exdir=dataDir)
germany2 = raster::shapefile(grep("shp", germany2, value=TRUE))

stuff = over(SpatialPoints(nuts[['2013']], proj4string = nuts[['2013']]@proj4string),
	germany2)

nuts[['2013']]$id2_2010 = stuff$FID
```


```{r downloadData}

Surl = c(outcome='https://ec.europa.eu/eurostat/estat-navtree-portlet-prod/BulkDownloadListing?file=data/demo_r_magec3.tsv.gz',
	pop = 'https://ec.europa.eu/eurostat/estat-navtree-portlet-prod/BulkDownloadListing?file=data/demo_r_pjangrp3.tsv.gz',
	employment= 'https://ec.europa.eu/eurostat/estat-navtree-portlet-prod/BulkDownloadListing?file=data/nama_10r_3empers.tsv.gz',
	manure = 'https://ec.europa.eu/eurostat/estat-navtree-portlet-prod/BulkDownloadListing?file=data/aei_fm_ms.tsv.gz',
	gdp = 'https://ec.europa.eu/eurostat/estat-navtree-portlet-prod/BulkDownloadListing?file=data/nama_10r_3gdp.tsv.gz'
)

Sfiles = mapply(
	Pmisc::downloadIfOld,
	url=Surl, 
	MoreArgs=list(path=dataDir))
Stables = mapply(
	read.table,
	file = Sfiles,
	MoreArgs = list(header=TRUE, 
	sep='\t',na.strings=':',
	stringsAsFactors=FALSE)
)
euTable = function(x, subset = NULL, ...) {
popTableRows = scan(what='a', text=x[,1], 
	sep=',',
	quiet=TRUE)
popTableRows = matrix(popTableRows, nrow=nrow(x),
	byrow=TRUE)
colnames(popTableRows) = unlist(
	strsplit(names(x)[1], '[.]'))[1:ncol(popTableRows)]

x = cbind(
	as.data.frame(popTableRows),
	x[,-1])
for(D in names(subset)) {
	x = x[x[,D] %in% subset[[D]], ]
}
x = x[, setdiff(colnames(x), names(subset))]
result = reshape2::melt(x,
	id.vars = grep("^X[[:digit:]]+$", colnames(x), invert=TRUE, value=TRUE),
	measure.vars = grep("^X[[:digit:]]+$", colnames(x), value=TRUE),  ...)
result$year = as.numeric(gsub("^X", "", result$variable))
result
}
Stables2 = mapply(
	euTable,
	x = Stables,
	value.name = names(Stables),
	subset = list(list(), list(), 
		list(wstatus='EMP', nace_r2 = 'TOTAL'), 
		list(unit='PC', indic_ag='I07A_EQ_Y'), 
		list(unit='PPS_HAB_EU27_2020'))#PPS_HAB_EU'))
	)

# copy 2013 codes to manure
manure2010 = manure3 = manure2 = Stables2$manure
manure2$geo = nuts$'2013'@data[
	match(as.character(manure2$geo),
		as.character(nuts$'2013'$id2016)), 'FID']
manure2 = manure2[!is.na(manure2$geo),]
manure3$geo = nuts$'2013'@data[
	match(as.character(manure3$geo),
		as.character(nuts$'2013'$FID)), 'id2016']
manure3 = manure3[!is.na(manure3$geo),]
manure2010$geo = nuts$'2013'@data[
	match(as.character(manure2010$geo),
		as.character(nuts$'2013'$id2010)), 'FID']
manure2010 = manure2010[!is.na(manure2010$geo),]

# nuts2, for Germany
manureDE = nuts$'2013'@data[ ,c('FID','id2010','id2016','id2_2010')]
manureDE = merge(manureDE, Stables2$manure, by.x='id2_2010', by.y='geo')
manureDE$geo = manureDE$FID
manureDE = manureDE[!manureDE$geo %in% manure2010$geo, ]

Stables2$manure[grepl("FR825|FRL05", Stables2$manure$geo), ]
Stables2$manure = rbind(Stables2$manure, manure3)
Stables2$manure[grepl("FR825|FRL05", Stables2$manure$geo), ]
Stables2$manure = rbind(Stables2$manure, manure2)
Stables2$manure[grepl("FR825|FRL05", Stables2$manure$geo), ]
Stables2$manure = rbind(Stables2$manure, manure2010)
Stables2$manure[grepl("FR825|FRL05", Stables2$manure$geo), ]
Stables2$manure = rbind(Stables2$manure, manureDE[,colnames(Stables2$manure)])
Stables2$manure = Stables2$manure[!is.na(Stables2$manure), ]
Stables2$manure = Stables2$manure[grep("[:]", Stables2$manure$manure, invert=TRUE), ]
Stables2$manure = Stables2$manure[!duplicated(
	Stables2$manure[,c('geo','year')]), ]
Stables2$manure[grepl("FR825|FRL05", Stables2$manure$geo), ]

nuts$'2013'@data[grep("FR825",nuts$'2013'$FID),]
manure3[grepl("FR825|FRL05", manure3$geo), ]
bob = Stables2$manure[grepl("FR825|FRL05", Stables2$manure$geo), ]
bob[order(bob$year), ]


Stables4 = list()
for(Dtable in setdiff(names(Stables2), 'outcome')) {
	yearsHave = unique(Stables2[[Dtable]]$year)
	yearsDontHave = setdiff(Syear, yearsHave)
	distYears = outer(yearsDontHave, yearsHave, '-')
	distYears = yearsHave[apply(distYears, 1, which.min)]
	toBind = list()
	for(Dyear in seq(1, len=length(yearsDontHave))) {
		toBindHere = Stables2[[Dtable]][Stables2[[Dtable]]$year == distYears[Dyear], ]
		toBindHere$year = yearsDontHave[Dyear]
		toBind[[Dyear]] = toBindHere
	}
	Stables4[[Dtable]] = do.call(rbind, toBind)
}

Stables3 = Stables2
for(D in names(Stables4))
	Stables3[[D]] = rbind(Stables2[[D]], Stables4[[D]])

Stables3$manure[grepl("FR825|FRL05", Stables3$manure$geo), ]


# fill employment
Stables3$employment$employment = as.numeric(gsub("[:]|[[:alpha:]]", "", 
	Stables3$employment$employment))
Stables3$employment = Stables3$employment[Stables3$employment$year %in% Syear, ]
missingTable = table(
	Stables3$employment$geo, 
	is.na(Stables3$employment$employment))

theEmpMean = tapply(Stables3$employment$employment, Stables3$employment$geo, mean, na.rm=TRUE)
Stables3$employment[is.na(Stables3$employment$employment), 'employment'] =
	theEmpMean[Stables3$employment[is.na(Stables3$employment$employment), 'geo']]

# fill gdp


Stables3$gdp$gdp = as.numeric(gsub("[:]|[[:alpha:]]", "", 
	Stables3$gdp$gdp))
Stables3$gdp = Stables3$gdp[Stables3$gdp$year %in% Syear, ]



missingTable = table(
	Stables3$gdp$geo, 
	is.na(Stables3$gdp$gdp))

theGmean = tapply(Stables3$gdp$gdp, Stables3$gdp$geo, mean, na.rm=TRUE)
Stables3$gdp[is.na(Stables3$gdp$gdp), 'gdp'] =
	theGmean[Stables3$gdp[is.na(Stables3$gdp$gdp), 'geo']]


# correct UK, ireland regions
nuts[['2013']]@data[grep("Dumfries|Border", nuts[['2013']]$NUTS_NAME), ]
nuts[['2016']]@data[grep("Dumfries|Border", nuts[['2016']]$NUTS_NAME), ]
Stables3$outcome[grepl("IE0[14]1|UKM[39]2|UKM64", Stables3$outcome$geo) & Stables3$outcome$sex == 'F'
	& Stables3$outcome$age == 'Y_LT5',]
Stables3$gdp[grep("IE0[14]1|UKM[39]2|UKM64", Stables3$gdp$geo),]



# gdp uses 2016 codes, outocme uses 2013


nuts$'2013'@data[grep("^IE|UKM[39]2", nuts$'2013'$FID),]


	i2013 = nuts[['2013']][grep("IE", nuts[['2013']]$FID), ]
	plot(i2013, lwd=5)
	text(SpatialPoints(i2013), i2013$FID, pos=3, col='blue')

	i2016 = nuts[['2016']][grep("IE", nuts[['2016']]$FID), ]
	plot(i2016, lwd=1, border='red', add=TRUE)
	text(SpatialPoints(i2016), i2016$FID, pos=1, col='green')

toOver = SpatialPoints(i2016)
toOver@proj4string = i2013@proj4string
cbind(i2016@data, over(toOver, i2013))

nuts[['2013']]@data[
	match(
		c('South-East (IE)', 'South-West (IE)'),
		nuts[['2013']]$NUTS_NAME),  'id2016'] = c('IE052', 'IE053')

(needToSwap2013 = grep("UKM|IE",unique(setdiff(Stables3$outcome$geo, Stables3$gdp$geo)), value=TRUE))
(needToSwap2016 = nuts[['2013']]@data[nuts[['2013']]$FID %in% needToSwap2013,'id2016'])

Stables3$gdp$geo = as.character(Stables3$gdp$geo)
toSwap = Stables3$gdp$geo %in% needToSwap2016
gdpNew = Stables3$gdp[toSwap,]

nuts$'2013'@data[
	match(as.character(Stables3$gdp[toSwap, 'geo']),
		as.character(nuts$'2013'$id2016)), c('FID','id2016','NUTS_NAME')]
gdpNew$geo = nuts$'2013'@data[
	match(gdpNew$geo,
		as.character(nuts$'2013'$id2016)),

	'FID']
Stables3$gdp = rbind(Stables3$gdp, gdpNew)

grep("IE01.", Stables3$gdp$geo, value=TRUE)

grep("UKM|IE",unique(setdiff(Stables3$outcome$geo, Stables3$gdp$geo)), value=TRUE)

popT = Stables3$pop[Stables3$pop$age == 'TOTAL' & Stables3$pop$sex == 'T', ]
popT$pop = as.numeric(gsub("[[:alpha:]]|[:]", "", popT$pop))
Stables3$employment = merge(
	Stables3$employment,
	popT[,c('geo','pop','year')]
	)

Stables3$employment$pct = pmin(1,1000*Stables3$employment$employment / Stables3$employment$pop)
Stables3$employment = Stables3$employment[,setdiff(names(Stables3$employment), 'pop')]



xFinal = merge(
	Stables2$outcome[,setdiff(names(Stables2$outcome), c('unit','variable'))], 
	Stables3$pop[,setdiff(names(Stables2$pop), c('unit','variable'))],
	by = c('sex','age','geo','year'), all.x=TRUE)

xFinal = xFinal[nchar(as.character(xFinal$geo)) >= 5, ]
xFinal = xFinal[grep("XX$", xFinal$geo, invert=TRUE), ]

xFinal = xFinal[xFinal$age != 'TOTAL', ]
xFinal = xFinal[xFinal$age != 'UNK', ]
xFinal = xFinal[!is.na(xFinal$pop), ]
xFinal = xFinal[xFinal$pop != ': ', ]
xFinal = xFinal[xFinal$outcome != ': ', ]

xFinal$pop = as.numeric(gsub("[[:alpha:]]+$", "", xFinal$pop))
xFinal = xFinal[xFinal$pop > 0, ]


nuts$'2013'@data[grep("FR825",nuts$'2013'$FID),]
bob = xFinal[grepl("FR825|FRL05", xFinal$geo) & xFinal$age == 'Y20-24' & xFinal$sex == 'F', ]
bob[order(bob$year), ]

for(D in setdiff(names(Stables3), c('outcome','pop'))) {
	xFinal = merge(
		xFinal, 
		Stables3[[D]][,
			setdiff(names(Stables3[[D]]), c('unit','variable'))],
		by = c('geo','year'),
		all.x=TRUE, all.y=FALSE
		)
}
names(xFinal) = gsub("^pct$", "pctEmp", names(xFinal))
xFinal$id0 = substr(xFinal$geo, 1, 2)

bob = xFinal[grepl("FR825|FRL05", xFinal$geo) & xFinal$age == 'Y20-24' & xFinal$sex == 'F', ]
bob[order(bob$year), ]


```



```{r merge}


for(D in c('2013','2016')) {
	xFinal[[paste0('nuts', D)]] = match(as.character(xFinal$geo),
		as.character(nuts[[D]]$NUTS_ID))
}




xFinal$deaths = as.numeric(gsub("[[:alpha:]]", "", xFinal$outcome))
xFinal$manure= as.numeric(gsub("[[:alpha:]]", "", xFinal$manure))


xFinal$ageLow = as.numeric(gsub("^(_|[[:alpha:]])+|[-].*", "", xFinal$age))
xFinal = xFinal[xFinal$ageLow >= 20 & xFinal$ageLow < 40 & xFinal$sex != 'T', ]

xFinal$outcome = as.numeric(gsub("[[:alpha:]]", "", xFinal$outcome))

xFinal$in6 = xFinal$geo %in% nuts[['2016']]$FID
xFinal$in3 = xFinal$geo %in% nuts[['2013']]$FID
xFinal$inOnly6 = xFinal$in6 & !xFinal$in3
xFinal$id0 = substr(xFinal$geo, 1, 2)
xFinal$inNeither = !xFinal$in6 & !xFinal$in3
xFinal[xFinal$inNeither, 'inOnly6'] = NA


xFinal$sex = factor(xFinal$sex)
xFinal$age = factor(xFinal$age)

xFinal = xFinal[!is.na(xFinal$gdp), ]

nuts = nuts[c('2010','2013','2016')]

bob = xFinal[grepl("FR825|FRL05", xFinal$geo) & xFinal$age == 'Y20-24' & xFinal$sex == 'F', ]
bob[order(bob$year), ]

save(xFinal, nuts, file='nuts.RData')
```

```{r euroSpace}
load('nuts.RData')

euro = xFinal[!xFinal$inNeither & xFinal$year == 2016, ]

euro$logPop = log(euro$pop)

theRates = glm(deaths ~ offset(logPop) + 0+age:sex, family='poisson', data=euro)
euro$expected =  predict(theRates, euro[,c('logPop','age','sex')], type='response')

sum(euro$expected)
sum(euro$deaths)

euroMelt = reshape2::melt(euro, 
  id.vars = c('geo','gdp','pctEmp','manure','sex', 'inOnly6'),
  measure.vars = c('expected','deaths')
  )
xDf = reshape2::dcast(
  euroMelt,
  geo + gdp + pctEmp + manure + inOnly6~ variable + sex,
  fun.aggregate = sum
  )
xDf$smr =  xDf$deaths_F /xDf$expected_F
xDf$logExpected_F = log(xDf$expected_F)

geoIn3 = unique(xDf[!xDf$inOnly6, 'geo'])
geoIn6 = unique(xDf[xDf$inOnly6, 'geo'])

xDf$id0 = substr(xDf$geo, 1, 2)

nuts3 = nuts[['2013']]
nuts3 = nuts3[nuts3$FID %in% geoIn3, ]
nuts6 = nuts[['2016']]
nuts6 = nuts6[nuts6$FID %in% geoIn6, ]

euroNuts = raster::bind(nuts6, nuts3)
saveRDS(euroNuts, 'euroNuts.rds')
x = merge(euroNuts, xDf, by.x='FID', by.y='geo')
x = x[!x$CNTR_CODE %in% c("CY",'TR'), ]
names(x) = gsub("FID", 'geo', names(x))
saveRDS(x, 'euroSpace.rds')
```







deprivation
https://ec.europa.eu/eurostat/web/products-datasets/-/tgs00104

https://ec.europa.eu/eurostat/web/population-demography-migration-projections/data/database


deaths by age and sex, nuts 2 and 3

https://ec.europa.eu/eurostat/web/products-datasets/-/demo_r_magec3

boundaries
https://ec.europa.eu/eurostat/web/gisco/geodata/reference-data/administrative-units-statistical-units/nuts#nuts16

Uk mortality counts

'https://www.ons.gov.uk/peoplepopulationandcommunity/birthsdeathsandmarriages/deaths/datasets/deathsregisteredbyareaofusualresidenceenglandandwales'

Live births

'https://www.ons.gov.uk/peoplepopulationandcommunity/birthsdeathsandmarriages/livebirths/datasets/birthsbyareaofusualresidenceofmotheruk'

'https://www.ons.gov.uk/file?uri=/peoplepopulationandcommunity/birthsdeathsandmarriages/livebirths/datasets/birthsbyareaofusualresidenceofmotheruk/2016/areaofusualresidencetables2016.xls'
'https://www.ons.gov.uk/file?uri=/peoplepopulationandcommunity/birthsdeathsandmarriages/livebirths/datasets/birthsbyareaofusualresidenceofmotheruk/2010/areaofusualresidencetable_tcm77-233219.xls'

Boundaries

'http://geoportal.statistics.gov.uk/datasets/local-authority-districts-december-2018-uk-bfe/data'

'https://opendata.arcgis.com/datasets/e50a6790161743038d586e1a9494d9ff_1.zip?outSR=%7B%22wkid%22%3A27700%2C%22latestWkid%22%3A27700%7D'


```{r downloadUk, eval=FALSE}

dataDir = 'C:/Users/brownpa/store/portugal'
dir.create(dataDir, showWarning=FALSE)
Sfiles = c(
	'2007'='data-tables--2007.zip',
	'2008'='data-tables--2008.zip',
	'2009'='data-tables--2009.zip',
	'2010'='deathsarea201_tcm77-236752.xls',
	'2011'='deathsarea2011_tcm77-295437.xls',
	'2012'='deathsarea2012_tcm77-352910.xls',
	'2013'='deathsarea2013_tcm77-400185.xls',
	'2014'='deathsarea2014tcm77431322.xls',
	'2015'='deathsregisteredbyareaofusualresidence2015.xls',
	'2016'='deathsregisteredbyareaofusualresidence2016.xls',
	'2017'='deathsregisteredbyareaofusualresidence2017.xls')

Surl = paste0(
	'https://www.ons.gov.uk/file?uri=/peoplepopulationandcommunity/',
	'birthsdeathsandmarriages/deaths/datasets/',
	'deathsregisteredbyareaofusualresidenceenglandandwales/',
	names(Sfiles), '/',
	Sfiles
	)


Sdownload = Pmisc::downloadIfOld(
	Surl,
	path=dataDir,
	mode='wb'
	)
StoRead = grep("[.]xls$", Sdownload, value=TRUE)

mSheets = mapply(
	readxl::excel_sheets,
	path = StoRead
	)

table1a = lapply(mSheets, grep, value=TRUE, pattern='Table 1a')
table1a = names(table1a)[which(unlist(lapply(table1a, length))>0)]
table1aDat = mapply(readxl::read_excel, path=table1a, MoreArgs = list(sheet='Table 1a', 
	skip=11, col_names=FALSE))

table1aHead = mapply(readxl::read_excel, path=table1a, MoreArgs = list(sheet='Table 1a', 
	skip=5, col_names=FALSE))

table1aHead2 = list()
for(D in names(table1aHead)){
	xx = table1aHead[[D]]
	xx = as.matrix(xx)[1:5,]
	xx[is.na(xx)] = ""
	popCol = grep("Population", xx[1,])
	xx[1,seq(from=popCol, by=1, len=3)] = 'Population'
	deathCol = grep("Deaths", xx[1,])
	xx[1,seq(from=deathCol, by=1, len=3)] = 'Deaths'

	smrCol = grep("Age.stand", xx[3,])
	if(!length(smrCol)) {
		smrCol = grep("Age.standard", xx[2,])
		xx[2,smrCol] = ''
	}

		xx[3,seq(from=smrCol, by=1, len=3)] = 'SMR'

	table1aHead2[[D]] = gsub("[[:space:]]", "", apply(as.data.frame(xx), 2, paste, collapse=''))
	table1aHead2[[D]] = gsub("1(.?)000", "1000", table1aHead2[[D]])
	table1aHead2[[D]] = gsub("Number.thousands.", "", table1aHead2[[D]])
}


codes = as.data.frame(table1aDat[[grep("2017", names(table1aDat))]][,1:2])
codes = codes[grep("^[[:punct:]]?E0[[:alnum:]]+$", codes[,1]), ]
names(codes) = c('id','name')

table1aa = list()
for(D in names(table1aDat)) {
	xx = as.data.frame(table1aDat[[D]])
	xxHead = table1aHead2[[D]]

	matchHead = match(names(xxHead), names(xx))
	names(xx)[matchHead] = xxHead

	if(any(names(xx)=='Areacodes')) {
		xxres = xx[xx$Areacodes %in% codes$id, ]
		names(xxres)[1:2] = names(codes)
	} else {
		xxres = merge(codes, xx, by.x='name', by.y=1, all.x=TRUE, all.y=FALSE)
	}
	year = gsub("^.*20", "20", D)
	year = gsub("([[:alpha:]]|[[:punct:]]).*", "", year)
	year = gsub('^201$', '2010', year)
	xxres$year = as.numeric(year)
	table1aa[[D]] = xxres
}

toKeep = names(which(table(unlist(lapply(table1aa, names)))  == length(table1aa)))

table1c = do.call(
	rbind,
	lapply(table1aa, function(xx) xx[,toKeep])
	)
rownames(table1c) = 1:nrow(table1c)

table1c= table1c[,grep("year|name|^id$|Infant[(]|Neonatal[(]", names(table1c))]



mSheets = lapply(mSheets, grep, value=TRUE, pattern="Information|Contents|Notes|Terms|Related", invert=TRUE)


mTables = mapply(
	readxl::read_excel,
	path = StoRead,
	MoreArgs = list(sheet=1:7)
	)

```




# populations

'https://www.ons.gov.uk/peoplepopulationandcommunity/populationandmigration/populationestimates/datasets/lowersuperoutputareamidyearpopulationestimates'

# deaths

'https://www.ons.gov.uk/peoplepopulationandcommunity/birthsdeathsandmarriages/deaths/adhocs/007931numberofdeathsregisteredbyspecificcausegroupssexandagespecifiedlowerlayersuperoutputarealsoagroupsinenglandandwales2001to2016'

# boundaries

'https://data.gov.uk/dataset/fa883558-22fb-4a1a-8529-cffdee47d500/lower-layer-super-output-area-lsoa-boundaries'


