Instructions
First, load the library:
library(networkscaleup)We will simulate a data set from the Binomial formulation in Killworth, P. D., Johnsen, E. C., McCarty, C., Shelley, G. A., and Bernard, H. R. (1998).
set.seed(1998)
N_i = 50
N_k = 5
N = 1e5
sizes = rbinom(N_k, N, prob = runif(N_k, min = 0.01, max = 0.15))
degrees = round(exp(rnorm(N_i, mean = 5, sd = 1)))
ard = matrix(NA, nrow = N_i, ncol = N_k)
for(k in 1:N_k){
ard[,k] = rbinom(N_i, degrees, sizes[k] / N)
}
## Create some artificial covariates for use later
x = matrix(sample(1:5, size = N_i * N_k, replace = T),
nrow = N_i,
ncol = N_k)
z = cbind(rbinom(N_i, 1, 0.3), rnorm(N_i), rnorm(N_i), rnorm(N_i))We should also prepare the data for modeling by scaling all covariates to have standard deviation 1 and mean 0. Additionally, we need to define which columns of \(Z\) belong to \(z_{subpop}\) and \(z_{global}\).
x = scale(x)
z = scale(z)
z_subpop = z[,1:2]
z_global = z[,3:4]PIMLE
The plug-in MLE estimator from Killworth, P. D., Johnsen, E. C., McCarty, C., Shelley, G. A., and Bernard, H. R. (1998) is a two-stage estimator that first estimates the degrees for each respondent \(d_i\) by maximizing the following likelihood for each respondent: \[\begin{equation} L(d_i;y, \{N_k\}) = \prod_{k=1}^L \binom{d_i}{y_{ik}} \left(\frac{N_k}{N} \right)^{y_{ik}} \left(1 - \frac{N_k}{N} \right)^{d_i - y_{ik}}, \end{equation}\] where \(L\) is the number of subpopulations with known \(N_k\). For the second stage, the model plugs in the estimated \(d_i\) into the equation \[\begin{equation} \frac{y_{ik}}{d_i} = \frac{N_k}{N} \end{equation}\] and solves for the unknown \(N_k\) for each respondent. These values are then averaged to obtain a single estimate of \(N_k\).
To summarize, stage 1 estimates \(\hat{d}_i\) by \[\begin{equation} \hat{d}_i = N \cdot \frac{\sum_{k=1}^L y_{ik}}{\sum_{k=1}^L N_k} \end{equation}\] and then these estimates are used in stage 2 to estimate the unknown \(\hat{N}_k\) by \[\begin{equation} \hat{N}_k^{PIMLE} = \frac{N}{n} \sum_{i=1}^n \frac{y_{ik}}{\hat{d}_i} \end{equation}\]
These estimates are obtained using the following call to the
killworth function.
pimle.est = killworth(ard,
known_sizes = sizes[c(1, 2, 4)],
known_ind = c(1, 2, 4),
N = N, model = "PIMLE")
#> Warning in killworth(ard, known_sizes = sizes[c(1, 2, 4)], known_ind = c(1, :
#> Estimated a 0 degree for at least one respondent. Ignoring response for PIMLE
plot(degrees ~ pimle.est$degrees, xlab = "Estimated PIMLE degrees", ylab = "True Degrees")
abline(0, 1, col = "red")
round(data.frame(true = sizes[c(3, 5)],
pimle = pimle.est$sizes))
#> true pimle
#> 1 1817 1645
#> 2 14134 14632Note that the function provides a warning saying that at least \(\hat{d}_i\) was 0. This occurs when a respondent does not resport knowing anyone in the known subpopulations. This is an issue for the PIMLE since a 0 value is in the denominator for \(\hat{N}_u^{PIMLE}\). Thus, we ignore the responses from respondents that correspond to \(\hat{d}_i = 0\).
MLE
Next, we analyze the data from the Killworth, P. D., McCarty, C., Bernard, H. R., Shelley, G. A., and Johnsen, E. C. (1998) MLE estimator. This is also a two-stage model, which an identical first stage, i.e. \[\begin{equation} \hat{d}_i = N \cdot \frac{\sum_{k=1}^L y_{ik}}{\sum_{k=1}^L N_k}. \end{equation}\] However, the second stage estimates \(\hat{N}_k\) by maximizing the Binomial likelihood with respect to \(\hat{N}_k\), fixing \(d_i\) at the estimated \(\hat{d}_i\). Thus, the estimate for the unknown subpopulation size is given by \[\begin{equation} \hat{N}_k^{MLE} = N \cdot \frac{\sum_{i=1}^n y_{ik}}{\sum_{i=1}^n \hat{d}_i}. \end{equation}\]
These estimates are also obtained using a single call to the
Killworth function.
mle.est = killworth(ard,
known_sizes = sizes[c(1, 2, 4)],
known_ind = c(1, 2, 4),
N = N, model = "MLE")
plot(degrees ~ mle.est$degrees, xlab = "Estimated MLE degrees", ylab = "True Degrees")
abline(0, 1, col = "red")
round(data.frame(true = sizes[c(3, 5)],
pimle = mle.est$sizes))
#> true pimle
#> 1 1817 1526
#> 2 14134 12907Note that there is no warning here since the denominator depends on the summation of \(\hat{d}_i\).
