Colleagues,

I've been asked to help one of the professors in my department do some
secondary data analysis on data from a previous study. I'm the only one on
the team with any experience in spatial analysis, but I'm still new to it
myself (took 1 graduate course). Still, she's asked some interesting
questions and wants to know if the data we have can answer them. I'd like to
tap into the expertise available here on this list for advice and feedback
on the analysis approach I've looked at so far. 

Here's a brief description of the context. As part of a large HIV prevention
study, we've got a dataset where 774 persons were sampled from a set of
social venues in a single large Midwestern city in the US. The target
population of study was young (age 15-25) African American males who have
sex with males. Samples were taken during randomly selected 4-hour blocks of
time distributed across different venues and different days of the week.
Individuals in the venue were approached, screened to see if they were in
the target population, and if so then took a brief survey (i.e., this is a
venue-day-time sampling strategy). Due to the need for anonymity, the only
location data we have for each respondent is the zip code where he resided
at the time of the interview. 

There are a total of 57 zip codes in the city, though our sample only
contains respondents from 44 zip codes. Given the way the general African
American population is concentrated in certain areas and almost non-existent
in others, this was no surprise to us. As you might expect, the number of
respondents from each zip code varies considerably (see summary below). 

> summary(CITY.ZCTA.data$all.waves.cases)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   0.00    1.00    3.00   13.58   20.00   75.00 

For each individual, we have an indicator variable that tells us whether or
not he engaged in unprotected sex within the last 3 months (see summary
below), and another that tells us whether or not he had contact with an
outreach prevention worker. 

> summary(CITY.ZCTA.data$risky9.yes
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  0.000   0.000   1.000   3.456   6.000  24.000 

We'd like to be able to calculate a meaningful rate or proportion of how
many individuals in the target population have been engaging in risky
behavior in each zip code so that I can map that and also relate it to the
density of HIV prevention service organizations (number per sq. mile) in
each zip code. We'd also like to see how the proportion of individuals
who've had contact with outreach workers relates to service organization
density. 

Here are the dilemmas I've identified so far. I'd appreciate any advice you
could offer to help me out with doing a reasonable analysis.

1. With numbers of respondents so small for so many zip codes, do we even
have enough data for a meaningful analysis? The sampling design was not
designed to facilitate spatial analysis, so if this sort of secondary
analysis is not viable and defensible for peer-reviewed publication, we'd
rather know that before we invest much more time in it. 

2. Assuming we have enough data, I'm not sure whether it's better to use the
total respondents from each zip code as the denominator, or whether we
should instead use the US Census data on the total number of African
American males age 15-25 as the denominator. In the examples I've seen for
analyzing zonal data, it seems like the data sources are usually more
comprehensive than our sampling method and more likely to capture all
"events of interest" that should be in the numerator. For example, if one is
using number of SIDS deaths divided by number of births in a set of counties
for the proportion, you could get fairly comprehensive cause of death data
for infants from death certificates. Our sampling method obviously doesn't
give us a chance to detect all persons in the target population of each zip
code who have engaged in risky behavior, so it seems to me like the number
of respondents is the more appropriate denominator. Is this reasonable? Are
there any better alternatives?

3. If we use the sampled individuals from the zip code as the denominator,
I'm trying to decide if the very small numbers from some zip codes means
that we should switch to analyzing empirical Bayesian rates instead of using
the raw proportion of people in order to smooth out problems with the
stability of the rate when it's estimated from very small numbers. Is that
appropriate in this case?

4. I tried calculating Bayesian rates based on formulas from Bailey &
Gatrell's 1995 book on Interactive Spatial Data Analysis (pp 304-306), but I
get two values that strike me as odd. In one case, the Bayesian rate end up
as a negative number (-0.27), and in another, it ends up as a value of about
1.5. Both of those seem like they should be impossible values. Shouldn't
these estimates range from 0 to 1?

> summary(CITY.ZCTA.data$bayes.risky9)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
-0.2702  0.1982  0.2401  0.2470  0.2626  1.4880 13.0000


Here's the code I used to calculate the Bayesian rates (the CITY.ZCTA.data
object is a SpatialPolygonsDataFrame).

# Get gamma hat - the pooled city-wide estimate of risky behavior in target
population
# (see Bailey & Gatrell 1995, pg 306)

gamma.hat.risky9 <-
sum(CITY.ZCTA.data$risky9.yes)/sum(CITY.ZCTA.data$all.waves.cases)

# Get the observed rate for each ZCTA (see Bailey & Gatrell 1995, pg 304)

CITY.ZCTA.data$obs.rate.risky9 <-
CITY.ZCTA.data$risky9.yes/CITY.ZCTA.data$all.waves.cases

# Calculate phi hat (see Bailey & Gatrell 1995, pg 306)

phi.hat.risky9 <-
(sum(CITY.ZCTA.data$all.waves.cases*(CITY.ZCTA.data$obs.rate.risky9 -
gamma.hat.risky9)^2, na.rm=TRUE) / 
            sum(CITY.ZCTA.data$all.waves.cases)) -
(gamma.hat.risky9/mean(CITY.ZCTA.data$all.waves.cases))

# Calculate the weight for each observation (see Bailey & Gatrell 1995, pg
306)

CITY.ZCTA.data$w.hat.risky9 <- phi.hat.risky9 / (phi.hat.risky9 +
(gamma.hat.risky9 / CITY.ZCTA.data$all.waves.cases))

# Calculate the Bayesian estimates using the weighting formula (see Bailey &
Gatrell 1995, pg 304)

CITY.ZCTA.data$bayes.risky9 <-
CITY.ZCTA.data$w.hat.risky9*CITY.ZCTA.data$obs.rate.risky9 + 
                               ((1 -
CITY.ZCTA.data$w.hat.risky9)*gamma.hat.risky9)

All feedback and suggestions are welcome. 

Steven J. Pierce, M.S.
Doctoral Student in Ecological-Community Psychology
Department of Psychology
Michigan State University
E-mail: pierces1@msu.edu


	[[alternative HTML version deleted]]


