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
Best Answer
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?