rm(list=ls()) ##################################################### # Step 5 of the code by Schmertmann & Rau on life expectancy estimates # on the Kreis-level. # Purpose: Construct the hierarchical H and Z matrices for # the German county map. # # Level-j effects are Zj %*% v # The Rx1 vector of the additive impact on Kreis-level # parameters is Hj %*% Zj %*% v # License: GPL-2, https://www.gnu.org/licenses/old-licenses/gpl-2.0.html ##################################################### library(dplyr) ## geographic coding information (CRITICALLY DEPENDENT ON THE ORDER OF KREISE ## in the dths.expos data frame!! MUST MATCH!!) ## ## AGS1 is the state code ## AGS2 is the Regierungsbezirk code (0 = 'this state has no RBZs') ## AGS3 is the Kreis code load("../initialdata/geo.RData") tmp = geo.df %>% mutate( effect1 = ix1, parent1 = 1, effect2 = ix2, parent2 = (effect2 !=0) * ix1, effect3 = ix3, parent3 = (effect3 !=0) * ix2) nkreise = nrow(tmp) #---------------------------------- # LEVEL 1 #---------------------------------- sel = (tmp$effect1 > 0) parents = tmp$parent1[sel] effects = tmp$effect1[sel] H1 = matrix(0, nrow=nkreise, ncol=length(unique(effects)) ) ix = cbind( which(sel), effects) H1[ix] = 1 G1 = 1* (table(parents, effects) > 0) #---------------------------------- # LEVEL 2 #---------------------------------- sel = (tmp$effect2 > 0) parents = tmp$parent2[sel] effects = tmp$effect2[sel] H2 = matrix(0, nrow=nkreise, ncol=length(unique(effects)) ) ix = cbind( which(sel), effects) H2[ix] = 1 G2 = 1* (table(parents, effects) > 0) #---------------------------------- # LEVEL 3 #---------------------------------- sel = (tmp$effect3 > 0) parents = tmp$parent3[sel] effects = tmp$effect3[sel] H3 = matrix(0, nrow=nkreise, ncol=length(unique(effects)) ) ix = cbind( which(sel), effects) H3[ix] = 1 G3 = 1* (table(parents, effects) > 0) ### Z1 = eigen( crossprod(G1))$vectors[, -(1:nrow(G1))] Z2 = eigen( crossprod(G2))$vectors[, -(1:nrow(G2))] Z3 = eigen( crossprod(G3))$vectors[, -(1:nrow(G3))] save( G1, G2, G3, H1, H2, H3, Z1, Z2, Z3, file='../data/HZ matrices.Rdata' )