# Two-level hierarchical model using time-series cross sectional data?

A question from someone who is fairly new to hierarchical modeling, and I'm looking for the best approach within R, preferably with the package lme4, MCMCpack, or rjags using a BUGS document. I'm unsure of the best approach, so I would like some guidance.

I'm interested in creating a two-level hierarchical model with data that is cross-sectional, time series, and at the individual level merged with data from the group level. Let me explain the two datasets that were merged:

The group level dataset shows the number of police on duty within 100 different counties. This data was collected 10 different times, so there are 10 sets of 100 counties. I merged this group level county data with individual level demographic and crime data (merged by county of the individual).

The individual-level data also has an indicator (1 or 0) for each individual on whether or not they reported a crime during that period. This individual level data was also collected at 10 points in time -- so it matches up with the group level data -- but it is cross-section, not panel data (different people each period). This is my dependent variable, so I'm looking for a logit or probit approach.

Basically, I would like to create a hierarchical model with two levels: county and time, where the period variable (1-10) is nested within the counties (1-100). This seems fairly straight forward -- a two level nested model -- but my approaches up to this point have failed.

Based on a book (Gelman and Hill) recommended by a colleague, I feel like I have the understanding to program the basic hierarchical models in BUGS and lme4, but the book does not go into detail on the more complicated models like nested over time, and other references have not been useful.

Below is a truncated sample of what my data looks like in R. Any advice, recommendations on packages to use, and sample code for modelling is much appreciated!

counties <- c(1,1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,2,3,3,3,3,3,3,3,3,3) #only 3 counties for explanation purposes
police <- c(1,22,4,56,3,32,12,8,43,5,45,34,33,21,62,22,3,12,19,29,11,8,32,33,18,12,12) #number of police per county
personID  <- seq(1,27) # only 27 people for explanation purposes
period <- c(1,1,1,2,2,2,3,3,3,1,1,1,2,2,2,3,3,3,1,1,1,2,2,2,3,3,3) #only 3 periods for explanation purposes
age <- c(45,55,23,67,21,34,39,48,52,45,32,71,55,56,19,34,48,56,77,33,22,21,44,64,51,55,60) #an individual level predictor not in the hierarchy
crime <- c(1,1,0,0,1,1,0,0,1,0,0,1,1,0,1,0,0,1,1,0,0,0,0,1,0,1,0) #Dependent Var: Did individual report a crime? Yes/No

sample <- matrix(c(personID, period, age, crime, counties, police), nrow=27, ncol=6)
colnames(sample) <- c("personID", "period", "age", "crime", "counties", "police")

> sample
personID period age crime counties police
[1,]        1      1  45     1        1      1
[2,]        2      1  55     1        1     22
[3,]        3      1  23     0        1      4
[4,]        4      2  67     0        1     56
[5,]        5      2  21     1        1      3
[6,]        6      2  34     1        1     32
[7,]        7      3  39     0        1     12
[8,]        8      3  48     0        1      8
[9,]        9      3  52     1        1     43
[10,]       10      1  45     0        2      5
[11,]       11      1  32     0        2     45
[12,]       12      1  71     1        2     34
[13,]       13      2  55     1        2     33
[14,]       14      2  56     0        2     21
[15,]       15      2  19     1        2     62
[16,]       16      3  34     0        2     22
[17,]       17      3  48     0        2      3
[18,]       18      3  56     1        2     12
[19,]       19      1  77     1        3     19
[20,]       20      1  33     0        3     29
[21,]       21      1  22     0        3     11
[22,]       22      2  21     0        3      8
[23,]       23      2  44     0        3     32
[24,]       24      2  64     1        3     33
[25,]       25      3  51     0        3     18
[26,]       26      3  55     1        3     12
[27,]       27      3  60     0        3     12

• Hi, you may want to look at this doingbayesiandataanalysis.blogspot.com.au/2012/11/… and this stats.stackexchange.com/questions/28099/… to get you started, and before somebody comes along to give you a fuller answer. – Matt Albrecht Nov 28 '12 at 4:49
• Thanks for the response. I'm going to go through these in detail and will let you know -- it looks like something I've seen before though, except for the split-plot example on the blog. I've been comfortable with the Gelman book and the lmer function for about 1.5 years, so I'm not completely clueless. I've just can't seem to find a clear example using time-series cross-sectional data. Maybe I should post some of my coding attempts as well. Cheers. – Captain Murphy Nov 29 '12 at 0:42
• Hi @CaptainMurphy. Did you find a suitable solution to this code? I have a similar situation, where I have a multinomial response variable measured for each individual at x different time points for each individual (which varies between individuals). So the first level is the response variable through time, and the second level is to look at response clustered within individuals. I have both metric and nominal predictor variables too. If you managed to find a solution, I would appreciate any advice. Thanks – user3237820 Jul 26 '15 at 10:42
• @user3237820 sorry but no. I moved on from this project without a great solution. – Captain Murphy Aug 27 '15 at 18:24

I'm not 100% sure, so check carefully, and you'll probably want to put priors on the variances instead of the 1.0E-12, but maybe something like this?

model {
for( i in 1 : nData ) {
crime[i] ~ dbern( mu[i] )
logit(mu[i]) <- base + b0[counties[i], period[i]] + b1[counties[i], period[i]] * police[i]
}
base ~ dnorm( 0 , 1.0E-12)

for (j in 1 : nCounties){
for (k in 1 : nPeriod) {
b0[j, k] ~ dnorm(b0h[j], 1.0E-12)
b1[j, k] ~ dnorm(b1h[j], 1.0E-12)
}
}

for ( j in 1 : nCounties ) {
b0h[j] ~ dnorm( 0 , 1.0E-12 )
b1h[j] ~ dnorm( 0 , 1.0E-12 )
}
}