# Modeling count data, using count data as regressors

I want to study the relation between prevalence of high-status members of the holders of a family name in the past and the status of current holders of the same family name (I don't assume any direct ancestral lineage. I'm interested in the predictive power of a family name in current social status, given an approximation of the social status of the holders of the same name in the past).

I believe I have a kind of error-in-variables, and a strange multi-level modelling problem. Any help to formulate my problem would be much appreciated.

$Y_i$: My dependent variable. Number of people with family name $i$ observed in a group of high-status people (e.g., graduates of a top university with the family name $i$, between the years 1990-2000).

$E_i$: Population size of the family name $i$ for the given time period (e.g., total number of people at graduate-level age between 1990-2000).

$X_i$: Independent variable of interest. This is also a a count data that represents the number of people with $i$ in another time period in another group of high-status people (e.g., number of elected officials in high positions with the family name $i$

$Fi$: Population for $i$ in the earlier time period (e.g., total number of people eligible for election between 1900-1910).

Let's assume the rates are actually very low (e.g., $Y_i$ << $E_i$ and $X_i$ << $F_i$). So I'm modeling $Y_i$ in a Poisson regression, using $E_i$ as the exposure variable:

$$E[Y_i|X_i] \sim Poisson(\lambda_i E_i)$$ where $\lambda_i$ is the population normalized rate.

In R, I use

m0 = glmer(y ~ offset(log(e)) + x, family='poisson')


If the coefficient of $x$ turns out to be statistically significant and in the expected direction, I'll conclude that family name has a predictive power given its past and try to explain and research possible mechanisms that explain this relation.

However, this model doesn't capture the fact that $X_i$ itself is a noisy observation/proxy for social status. In this model, the number of people with a given family name $i$ in the older time period would have an effect on the uncertainty in $X_i$. With finite samples and low rates, the count data is most of the time 0 (underestimating the underlying status) or 1 (overestimating the underlying status).

1. Shouldn't my model take into account that family names with large $F_i$ have a more confident estimation of older-generation high-status? How can I incorporate this into my model?

Since $X_i$ is also a count data, I thought maybe I can fit a first-stage Poisson model and used intercepts of $i$ as the predictor variable:

m1 = glmer(x ~ offset(log(f)) + (1|i), family='poisson')


where $i$ is the unique identifier for the corresponding surname. Then for each $i$ I'd obtain a different intercept and could use that as the predictor variable in the next step.

2. But again, how do I propagate the uncertainty on these estimates to the next step?

3. Do I need to employ a multi-level approach? That also doesn't seem right because I don't really have different levels, all my observations are count data at the same unit with different populations.

4. Would instrumental variables approach be suitable in this case? I could find other proxies of high-status in the past for the holders of a given family name. For example, number of doctors, lawyers, people who owned land, who were noble, who were mentioned in historical texts, etc. These are presumably correlated with the underlying social status, but hopefully the error terms are uncorrelated between each other. If that's a reasonable assumption what are the practical next steps I can take?

• It sounds more like you have 2 binomials than 2 count variables. Also, in R, glm by default uses family=gaussian, so you'd need to re-specify that to do what you want, & glm can't take (1|i), you need glmer in lme4 for that. Sep 15 '16 at 0:38
• Thanks. I corrected the R code -- which is what I use. You are right that I have two binomials, but with N >> p (N is on the order of national population, number of successes (Np) is mostly 0 or 1, I thought treating it as a count data was better approach. Sep 15 '16 at 0:47
• If N>>p, then sure. But if you're counting the number in Harvard's graduating class, or something like that, it's pretty finite & binomial might be better. I guess it depends on exactly what the counts are. Sep 15 '16 at 0:52
• I'm considering number of Spearmans in Harvard's graduation class versus number of all Spearmans in the country (assuming the student pool of Harvard is US only). $p$ would be the probability that a Spearman is a graduate of Harvard. Sep 15 '16 at 0:58
• OK, I was thinking it was the number of Spearmans in Harvard's graduating class as a proportion of everyone in Harvard's graduating class. Sep 15 '16 at 1:15

You are comparing two rates. You should reorganize your data, $$X_i$$ isn't really a predictor, it is also a response observed under other conditions (another time period.)

So let $$C=(Y_1, \dotsc, Y_n, X_1, \dotsc, X_n)^T$$, $$P=(E_1, \dotsc, E_n, F_1, \dotsc, F_n)^T$$ and let $$\text{period}$$ be the time period indicator. $$n$$ is the number of family names. Then a very simple model is

mod0 <- glm(C ~ offset(log(P)+period, family=poisson)


but this does not represent the heterogeneity of effects of different names. So with $$\text{names}$$ as an $$n$$-level factor, used as random, we get

library(lme4)
mod1 <- glmer(C ~ offset(log(P)) + period+ (1|names), family='poisson')


corresponding to your proposed model. Your null hypothesis is that the coefficient of $$\text{period}$$ is zero.