# How do I find the mean of counts of reported events when the last category is “5 or more”

I know that we can find mean (No. of errors per patient) from 0*5+1*30+2*56+3*15+4*10

But for 5 or more, what do I need to do?

• @BruceET so, is it better to ignore the 5 or more and calculate for n = 116 – user125163 Sep 4 '19 at 2:03
• Maybe there is some convention for handing this in your field of application. I'd say mean is above 2.13. But it's risky: one of those 7 patients might have made 27 errors, so the mean might be considerably above 12.13. // People who really care about precision in such matters should take the mean before making an incomplete table, and this losing information. – BruceET Sep 4 '19 at 3:34
• ( @BruceET You might find this interesting also.) Related to mean and median estimated of binned quantities with 'X+' upper bounds: Angle, J. (2004). The Salamander: A Model of the Right Tail of the Wage Distribution Truncated by Topcoding [Working Paper]. Economic Research Service. If I recall aright, this paper is getting a tad dated, and more flexible methods have been developed for these types of estimation problems this decade. – Alexis Sep 4 '19 at 19:58

Comment: Something went wrong when I cut and pasted your computation from your question into R. So some numbers in my first Comment (just now deleted) were wrong. (It is best not to let such mistakes clutter the comment section.) My second comment has been edited to make the necessary corrections.

If the data values (numbers of errors) $$v_i$$ are in the vector $$(0, 1, 2, 3, 4, 5).$$ And the respective frequencies with which they occur (number of patients making those errors) $$f_i$$ are in the vector $$(5, 30, 56, 15, 10, 7),$$ then then average number $$A$$ of errors per patient is found in R as

$$A = \frac{\sum_i f_iv_i}{\sum_i f_i} = \frac{262}{123} = 2.130,$$

v = 0:5;  f = c(5, 30, 56, 15, 10, 7)
sum(f*v);  sum(f)
[1] 262
[1] 123
sum(f*v)/sum(f)
[1]  2.130081


The true sample mean, if you had the actual counts in the last cell, could have a slightly larger numerator. So the true sample mean could be somewhat larger than 2.13.

Speculation: Given the pattern of the sample values, where frequencies are largest around 1 and 2 errors per patient, it does not seem likely that the last cell contains (undisclosed) numbers that would make the mean much larger. (There is always the possibility that one of the patients had a huge number of errors that is an outlier. But, depending on the circumstance of that patient, you might want to omit that outlier from the computations and make a note of its removal.)

Note: Sometimes, it is possible to make a good guess what the distribution of errors is. If experience has shown that patients make 2 errors each, on average according to a Poisson distribution, then a simulated experiment with errors of 123 patients gives the following tabulation:

set.seed(2019)
table(rpois(123,2))

0  1  2  3  4  5  6  7
21 35 30 22  7  5  2  1


There are few (simulated) patients with '5 or more' errors, and the sample mean from this tabulation is about 1.89, which is not much different from the value 2 used to simulate the data.

A chi-squared goodness of fit test says your data are consistent with $$\mathsf{Pois}(\lambda=2.13).$$ P-value around 0.2--in spite of maybe too few 0's and too many 2's.

On average, in a sample of size $$n = 123$$ from this distribution one would expect about 8 patients with '5 or more' errors and essentially none with 8 or more. Of course none of this speculation gives direct knowledge of your data. But it may be reasonable to expect that the sample mean 2.13 is not far off.

• There's also weighted.mean(c(5, 30, 56, 15, 10, 7) ,0:5) – Acccumulation Sep 4 '19 at 22:46

You need to make additional assumptions. If you can make no assumptions, you can only find the lower bound for the mean (as in Kinformationst's answer). If you can put a strict upper bound on the number of errors you can also put an upper bound on the mean, (as in Anant Gupta's answer ). Additionally you can assume a specific distribution for the errors (e.g. Poisson as in BruceET's answer ). While BruceET describes a clever but slightly "dirty" solution (and there is nothing wrong with that), you can be more rigorous about this. The (a bit more) rigorous way is to fit a right-censored model - here "right-censored" is just a fancy name for the fact that the last category means "5 or more".

There are multiple R packages that can fit censored models for you, here is code using the VGAM package:

library(VGAM)
censoring_bound <- 5
#Simulate data
dd <- data.frame(y_raw = rpois(100, 2.3))
dd$$censored = dd$$y_raw >= censoring_bound
dd$$y <- ifelse(dd$$censored, censoring_bound, dd$$y_raw) dd$$status <- ifelse(dd$censored, 0, 1) #Here we have intercept only (~ 1), meaning we expect all the values to come from the same distribution fit <- vglm(SurvS4(y, status) ~ 1, family = cens.poisson, data = dd) #ML estimate of Poisson mean, exponentiated because of the log link exp(fit@coefficients["(Intercept)"]) #Compare to the mean of the observed numbers mean(dd$y)


With this approach you can also fit other families (e.g. negative binomial)

• Among the joys of this site: A really excellent idea for a better answer. (+1) – BruceET Sep 4 '19 at 22:36

If you have no knowledge of how the last category of error is distributed, then the only way would be to calculate the floor of the mean. If 5 or more was only 5, then the mean for errors per patient would be 262/123=2.13

so

$$\langle error_p patient\rangle \geqslant$$ 2.13

For completeness sake, here's a solution using a Bayesian model (using R and Stan). We're going to fit a censored Poisson to the data (see the documentation of Stan for more information). The model assumes that the values are right censored at $$5$$. To check whether the Poisson is a suitable model for the data, we're also going to generate a predictive distribution. The model uses a normal prior with mean $$2$$ and standard deviation $$1$$ for $$\lambda$$.

library(rstan)
library(vcd)

# The model
model <- "data {
int<lower=1> N;
int Y[N];
int cens[N];
real lambda_start_mu;
real lambda_start_sigma;
}

parameters {
real<lower = 0>lambda;
}

model {
lambda ~ normal(lambda_start_mu, lambda_start_sigma);

for(n in 1:N){
if (cens[n] == 0) {
target += poisson_lpmf(Y[n] | lambda);
} else if (cens[n] == 1) {
target += poisson_lccdf(Y[n] | lambda);
}
}
}

generated quantities {
int y_pred;
y_pred = poisson_rng(lambda);
}"

# Prepare the data
b <- rep(0:5, times = c(5, 30, 56, 15, 10, 7))
dat <- list(
Y = b
, N = length(b)
, cens = ifelse(b== 5, 1, 0)
, lambda_start_mu = 2
, lambda_start_sigma = 1
)

# Fit the model using Stan
fit <- stan(
model_code=model
, data = dat
, iter = 10000
, chains = 5
, warmup = 2000
, verbose = TRUE
, seed = 142857
)

# Plot the estimate and credible intervals (not shown)
plot(fit)

# Print the summary
print(
fit
, digits = 3
)

mean se_mean    sd     2.5%      25%      50%      75%    97.5% n_eff Rhat
lambda    2.214   0.001 0.135    1.958    2.122    2.212    2.304    2.487 14433    1
y_pred    2.219   0.008 1.503    0.000    1.000    2.000    3.000    6.000 38981    1
lp__   -199.276   0.006 0.721 -201.332 -199.440 -199.000 -198.818 -198.767 16511    1


The mean of the posterior distribution is $$\hat{\lambda}=2.21$$ and the 95% credible interval based on quantiles is $$(1.96; 2.49)$$. This is very comparable with solutions based on maximum likelihood estimation.

Let's use a rootogram using the predicted frequencies to check the model fit of a Poisson:

mcmc_chain <- as.matrix(fit)
preds <- mcmc_chain[, "y_pred"]

preds[preds >= 5] <- 5

obs <- c(table(dat$$x)) fitted <- prop.table(table(preds))*dat$$n

rootogram(
obs, fitted
, scale = "sqrt"
, type = "hanging"
, names = c(0:4, "5+")
)


The red points display the square root of the expected counts. The square roots of the observed counts are shown as bars. Ideally, all bars are close to the 0 line. We can see that the Poisson seems not entirely ideal: We would expect a little more $$0$$s and $$3$$s and a little less $$2$$s.

• Can you clarify if your model assumes right censoring at 5? It looks like yes, but not all of us are familiar with R or Stan. – Weiwen Ng Sep 5 '19 at 13:14

If we can take thresholds for the last column, our mean range will be

• 1st Range

• 2nd Range

We can see that the difference is around .13. One can then take a call to take

• The maximum of the two
• The mean of the two
• Of course, the mean could conceivably, if not likely, be much larger than 2.244. – BruceET Sep 4 '19 at 22:35

We can put a lower bound on the mean; if all of the remaining $$7$$ had $$5$$ errors, then the mean is $$2.130081$$. We can't put an upper bound, however: if the remaining $$7$$ all had $$1000$$ errors, for instance, then the mean would be $$58.7$$. However, we can look at what reasonable numbers would be. If the remaining people had an exponential drop-off of $$4$$ people making $$5$$ errors, $$2$$ making $$6$$, and $$1$$ making $$7$$, then the mean is $$2.162602$$ (while this is somewhat arbitrary, if the tail were decaying slower than a factor of $$02$$ after $$5$$, then it would have to be decaying faster than $$2$$ from $$4$$ to $$5$$). If we're really pessimistic and take $$1$$ person making $$5$$ errors, $$1$$ person making $$6$$, etc. up to $$1$$ person making $$11$$ errors, then the mean is $$2.300813$$. This an example of Knightian error: not only do we not know how much the error is, we don't know for certain what probability distribution the error is drawn from. So there's some degree of subjectivity as to what is appropriate as far as reporting the level of error. One would be on reasonable grounds to say that $$2.16$$ is a good guess, and it could be an overestimate by $$0.03$$, giving an estimate of $$2.16 \pm 0.03$$. If one wanted to be really conservative, one could note that the mean of $$2.3$$ is at the extreme end of plausibility, giving something along the lines of $$2.2 \pm 0.1$$. A danger in giving the results in this form, however, is the possibility that they would be understood as coming from some highly rigorous confidence interval calculation, which would be misleading (of course, most people's understanding of confidence intervals, even when calculated through "rigorous" methods, is rather inaccurate anyway).