Semi-parametric analysis of overdispersed count and metric data with varying follow-up times: Asymptotic theory and small sample approximations.

Count data are common endpoints in clinical trials, for example magnetic resonance imaging lesion counts in multiple sclerosis. They often exhibit high levels of overdispersion, that is variances are larger than the means. Inference is regularly based on negative binomial regression along with maximum-likelihood estimators. Although this approach can account for heterogeneity it postulates a common overdispersion parameter across groups. Such parametric assumptions are usually difficult to verify, especially in small trials. Therefore, novel procedures that are based on asymptotic results for newly developed rate and variance estimators are proposed in a general framework. Moreover, in case of small samples the procedures are carried out using permutation techniques. Here, the usual assumption of exchangeability under the null hypothesis is not met due to varying follow-up times and unequal overdispersion parameters. This problem is solved by the use of studentized permutations leading to valid inference methods for situations with (i) varying follow-up times, (ii) different overdispersion parameters, and (iii) small sample sizes.

INTRODUCTION

Metric data and especially count data are common endpoints in clinical trials. Examples include relapses and magnetic resonance imaging (MRI) lesion counts in relapsing‐remitting multiple sclerosis (MS), exacerbations in chronic obstructive pulmonary disease (COPD), and hospitalizations in heart failure. For several of these the negative binomial distribution has been suggested to be an appropriate model accounting for between‐patient heterogeneity in event rates manifesting in overdispersion, that is variances exceeding the means. For instance, Wang, Meyerson, Tang, and Qian (2009) suggested the negative binomial model for the analyses of relapses, and Sormani et al. (1999, 2001, 2005) and Van den Elskamp, Knol, Uitdehaag, and Barkhof (2009) for various types of MRI lesion counts in MS. Based on two large‐scale COPD trials, Keene, Calverley, Jones, Vestbo, and Anderson (2008) assessed various models and recommended the negative binomial model for application. In the situations described above, commonly analyses methods (e.g. PROC GENMOD in SAS) are applied based on large sample properties of underlying Maximum‐Likelihood‐Estimates (MLE) and the assumption of a common overdispersion parameter across treatment groups. Such distributional assumptions, however, can hardly be verified; especially in case of small to moderate sample sizes (Aban, Cutter, & Mavinga, 2009). Even if the distribution is correctly specified the MLEs of the overdispersion parameters are biased (Link & Sauer, 1997; Lord, 2006; Paul & Islam, 1995; Saha, 2011; Saha & Paul, 2005) that may lead to wrong conclusions. Moreover, it is quite common that varying follow‐up times occur, see for example, Chen et al. (2013), McCullagh and Nelder (1989). All of the above‐ mentioned characteristics may not only be shared by count data, but also by metric data measured on an arbitrary scale. Simultaneously accommodating all of these complications in an accurate statistical inference method in a unified way is a rather challenging task. To the best of our knowledge no suitable methods currently exist that can simultaneously handle heteroscedastic data (counts) with varying follow‐up times.

It is the aim of the present paper to develop valid inference procedures for the analysis of such data in general models allowing for possibly time‐varying follow‐up times and different overdispersion parameters in a nonparametric way. This is accomplished by newly derived unbiased estimators (based on the methods of moments) for the (count) rates and their variances. The rigorous study of their large sample properties then leads to asymptotically correct tests and confidence intervals for treatment effects using critical values from the standard normal distribution.

With small samples the use of normal quantiles for inference can lead to liberal or conservative decisions whereas permutation tests offer an opportunity to derive quantiles from appropriate reference distributions. In particular, the application of studentized permutation procedures is tempting since they have been shown to control the type‐I‐error rate very accurately in various situations (Chung & Romano, 2013; Chung & Romano, 2016; Janssen, 1997; Konietschke & Pauly, 2014; Pauly, Brunner, & Konietschke, 2015). The problem in this particular situation is that with varying follow‐up times and unequal overdispersion parameters the usual assumption of independently identically distributed (iid) observations in the groups is not met. This issue can be solved by applying more general theorems on permutation statistics by Janssen and Pauls (2003) and Janssen (2005) and Pauly (2011). Even though data may not be exchangeable under the null hypothesis, the derived permutation methods are asymptotically correct in that they control the type I error rate or the coverage probability for hypothesis tests and confidence intervals, respectively.

The paper is organized as follows: The statistical model and point estimates are given in Section 2. Unbiased variance estimators are provided in Section 3. In Section 4, test procedures and confidence intervals are derived. Permutation‐based small sample size approximations and simulation results are presented in Section 5. Finally, two illustrative data examples are analyzed in Section 6. The paper closes with a discussion of the proposed methods in Section 7. All proofs are given in the supplement to this paper.

STATISTICAL MODEL, POINT ESTIMATES, AND MULTIVARIATE NORMALITY

We consider a general semi‐parametric two‐sample layout with independent random variables with Here, the index i represents the treatment groups ( control, and treatment), and k the subject within treatment group i with individual follow‐up time , and the expectation of group i. Note that the variance may depend on , for example if follows a Negative Binomial distribution (in this special case ), a Poisson distribution (), or an Exponential distribution (here ). We further assume that the fourth moments exist and are bounded, that is for a constant and .

The design is allowed to be completely heteroscedastic, that is every observation might have a different expectation and variance. All statistical procedures for the analysis of iid observations are inappropriate for statistical inference in model (1). Let denote the total sample size, the total follow‐up times in group i, , and let denote the total follow‐up times across both treatment groups. The unknown rate parameters can be estimated without bias by and can be interpreted as a weighted mean of the data. The variance of is given by For the derivation of asymptotic results for the rate estimates (2), the following mild regularity conditions on sample sizes and follow‐up times are required: Assumption (4) ensures that the follow‐up times appear on a fixed time interval of interest, while Assumptions (5)–(7) guarantee the existence of limiting variances of the point estimates, see Theorem 2.1 below. In particular, it follows immediately, that the estimator is consistent as and , respectively. However, the variance defined in (3) represents an unknown weighted sequence of the quantities , which depends on both the follow‐up times and sample sizes. Thus, it cannot be represented by model constants. In order to derive inference methods for the general hypothesis , however, the estimator needs to be multiplied by adequate known coefficients, such that converges to a specific variance constant, which is, asymptotically, unaffected by the follow‐up times and sample sizes. The result along with the multivariate normality of the estimator of are given in the next theorem.

Under Assumptions (4), (6), and (7), is a diagonal limiting covariance matrix.

Note that the diagonal covariance matrix neither depends on the sample sizes , nor on the time‐varying coefficients . The matrix, is, however, unknown in practical applications, and needs to be estimated. An unbiased and ‐consistent estimator is derived in the next section.

ESTIMATION OF THE VARIANCE

Moment‐based estimators for variances denote, roughly speaking, the squared deviation from the mean. In model (1), however, no uniquely defined mean exists. In particular, the variance is a sum of variances, and is not defined as a fixed variance constant. Therefore, the usual sample variance moment‐based estimator is biased, a rather inappropriate characteristic of a variance estimator. Below, we derive an unbiased and consistent moment‐based estimator of .

Define the random variables , and note that for all , and . The variables describe the deviation of to its estimated expectation. An unbiased moment‐based estimator can now be derived by considering the squared deviation from along with a bias correction. Define and consider The estimator is not a usual sample variance estimator, since it only involves sums of the follow‐up times as weighting factors. However, it describes the mean squared deviation from the observations to their estimated mean . Further let denote the diagonal matrix with diagonal elements and , respectively. It is shown in the next theorem, that is an unbiased estimator of and that is ‐consistent.

For each the estimator is an unbiased estimator of . Moreover, the estimator is ‐consistent, that is

A detailed proof is given in the supplementary material.

We note that the variance estimator may become negative in “severe” situations, that is if any is way larger than the others. In this case we suggest to use the asymptotically unbiased version of instead.

The asymptotic normality of the point estimates and the consistent variance estimates can now be used for the derivation of test procedures and confidence intervals.

TEST PROCEDURES AND CONFIDENCE INTERVALS

In this section, different procedures for testing the null hypothesis as well as confidence intervals for the treatment effect will be discussed, where is continuously differentiable in . Let denote the gradient of h with estimator . It follows from the multivariate delta‐method that where The variance is unknown, and must be estimated in practical applications. However, is a linear combination of the individual variances , respectively. It follows immediately, that a consistent estimator is given by Based on the asymptotic normality of and Slutsky's Theorem, it thus follows that where . For large sample sizes, the null hypothesis will be rejected at a two‐sided significance level α, if , where denotes the ‐quantile of the standard normal distribution. Asymptotic ‐confidence intervals for θ are obtained from

SMALL SAMPLE APPROXIMATIONS AND SIMULATION RESULTS

Extensive simulations were conducted to investigate the accuracies of the test procedures derived in Section 4 for small sample sizes with regard to (i) controlling the type‐1 error rate at the nominal significance level (), (ii) their powers to detect certain alternatives , and (iii) the coverage probabilities of the corresponding confidence intervals in (16). All simulations were conducted with R environment, version 2.15.2. (R Development Core Team, 2010), each with simulation runs.

In all simulations, we focus on testing the hypothesis corresponding to the function . The test statistic is given by which yield to asymptotically valid tests for . Moreover, confidence intervals can be derived from (16), respectively. Simulation studies indicate, however, that the statistic in (18) tends to result in rather liberal conclusions for small sample sizes (). Therefore, we propose a studentized permutation approach to approximate its sampling distribution for small sample sizes. This will be explained in the next section.

A studentized permutation approach

Permutation tests are widely known to be robust and exact level α tests when the data are exchangeable. Exchangeability implies, however, that variances across the groups are identical. As mentioned above, the data are allowed to be completely heteroscedastic in model (1). Roughly speaking, a usual permutation test would fail to test the null hypotheses formulated above. However, asymptotic permutation tests can be obtained, if appropriate studentized statistics are permuted, which will now be briefly explained: It turns out that the test statistic follows, asymptotically, a standard normal distribution under the null hypothesis. A permutation or resampling test would now lead to accurate results (at least asymptotically), if the conditional permutation distribution of the test statistic , say , would generally mimick the null distribution of the test statistic. That is, both distributions should at least coincide asymptotically. If that is the case, critical values (or P‐values) could be computed from the permutation distribution instead of the standard normal distribution for making inferences. Therefore, the goal of the following investigations is to show that the permutation distribution of , , is indeed the standard normal distribution. In order to do so, some notations and ideas about the permutation schemes are necessary:

Let denote the pooled sample, and let denote the corresponding vector of the pooled follow‐up times . For a fixed, but random permutation π of , let and denote the permuted data and corresponding follow‐up times, respectively.

Permuting and using the same random permutation π, the permuted values and are not necessarily independent, which is a rather (at least technically) undesirable property in this context. We therefore propose to permute and independently. This is similar to two sample problems with right‐censored survival data, where it is also recommended that the permuted failure times do not occur in general with their corresponding censoring indicators, see Janssen and Mayer (2003) as well as Brendel, Janssen, Meyer, and Pauly (2014). To this end, we consider another random permutation of that is independent of π and calculate the permuted estimators and = . Note that the possible number of random permutation is considerably increased when permuting both and independently.

It turns out that the distribution of the test statistic differs in the general model (1) from its permutation distribution, and a valid level α test can not be achieved in this setup. Therefore, we consider the distribution of the test statistic defined in (15) and of the studentized quantity The conditional limiting distribution of given the data will be derived in the next theorem.

Let as given in (19) and denote by the standard normal distribution function. If , then we have convergence under the null as well as under the alternative with

Theorem 5.1 states that the limiting standard normal distribution of does not depend on the distribution of the data, particularly, it is achieved for arbitrary , that is it even holds under the alternative.

Let , where denotes the ‐quantile from the studentized permutation distribution of . In the next theorem, we will show that both the conditional and unconditional tests are asymptotically equivalent, which means, that both tests have, asymptotically, the same power to detect certain alternatives.

Suppose that the assumptions of Theorem 5.1 are fulfilled.

Under the null hypothesis , the studentized permutation test is asymptotically exact at α level of significance, that is , and asymptotically equivalent to , that is

The permutation test is consistent, that is we have convergence

In particular, Theorem 5.1 states that the distributions of the pivotal quantity and of the studentized permutation statistic asymptotically coincide. Under the assumptions of Theorem 5.1, approximate ‐confidence intervals for θ can be obtained from

Simulation results

In a negative binomial‐‐model we investigate the empirical control of the preassigned type‐1 error rate at the usual two‐sided significance level of the statistic in (18) using the standard normal approximation as given in (15), and the permutation test using the quantiles of the conditional distribution of in (19) as critical values. As a further competing procedure, we estimate the variances using maximum likelihood methods. In this ‐model the variance is given by the weighted sequence of the quantities , respectively. An intuitive plug‐in estimation approach is achieved by replacing the unknown parameter by from above and by a consistent maximum‐likelihood estimator (ML) , for example by using see, for example Schneider, Schmidli, and Friede (2013). This estimation approach, however, has the disadvantage that neither nor are unbiased estimators of or , respectively, resulting in biased variance estimators. The variance estimators used in are finally replaced by , and the corresponding Wald‐statistic, which is asymptotically equivalent to the Likelihood‐ratio test, denoted by LRT.

Type‐1 error rate simulations

We explore the behavior of the test statistics for smaller and larger effect rates λ1 and λ2 as well as smaller and larger overdispersion parameters ϕ1 and ϕ2 .

All simulation designs are motivated by the examples presented in Section 6. A major assessment criterion for the accuracy of the procedures is their behavior when increasing sample sizes are combined with increasing variance parameter constellations (positive pairing) or with decreasing variances (negative pairing). We investigate balanced situations with sample size vector and unbalanced situations with sample size vector . The sample sizes are increased by adding a constant m to the components of the vectors or , respectively. The different simulation settings are displayed in Table 1. Each simulation setting represents a different design with an increasing sample size m, where , see Table 1.

Table 1

Simulated designs, where and ,

Setting	λ_1=λ_2	Sizes	Overdisp.	Interpretation
1	1.5	n=n_1+m	ϕ=ϕ_1	Balanced/equal overdispersion
2	1.5	n=n_2+m	ϕ=ϕ_1	Unbalanced/equal overdispersion
3	1.5	n=n_1+m	ϕ=ϕ_2	Balanced/unequal overdispersion
4	1.5	n=n_2+m	ϕ=ϕ_2	Unbalanced/unequal overdispersion (positive pairing)
5	1.5	n=n_2+m	ϕ=ϕ_3	Unbalanced/unequal overdispersion (negative pairing)
6	10	n=n_1+m	ϕ=ϕ_4	Balanced/equal overdispersion
7	10	n=n_2+m	ϕ=ϕ_4	Unbalanced/equal overdispersion
8	10	n=n_1+m	ϕ=ϕ_5	Balanced/unequal overdispersion
9	10	n=n_2+m	ϕ=ϕ_5	Unbalanced/unequal overdispersion (positive pairing)
10	10	n=n_2+m	ϕ=ϕ_6	Unbalanced/unequal overdispersion (negative pairing)

Here , , , , , and denote vectors of overdispersion parameters and means that every component of , that is each group size, is increased by m.

Data were generated from , where denotes the realization from a uniformly distributed random variable , respectively. For each simulation setting, the same generated follow‐up times were used for the simulation runs, but they were newly generated for each design. The number of random permutations was set to . The simulated type‐1 error rates for a significance level assuming uniformly distributed follow‐up times are displayed in Figure 1.

Figure 1

Type‐I error level (α = 5%) simulation results (y‐axis) of the statistics in (18), permutation test in (19) and ML‐based statistics for different distributions, sample size increments (x‐axis), where denote the realizations from . The simulation settings are described in Table 1

It turns out that in case of small effect rates () and small overdispersion parameters the statistics based on the normal approximation as well as the LRT statistics based on ML tend to be slightly liberal. It can be readily seen from Figure 1 that the permutation tests control the type‐1 error rate best, even for extremely small sample sizes. In case of larger effect rates and overdispersion parameters the distribution of the data is much more skewed. In these situations the procedures based on the normal approximation and ML tend to considerably overreject the null hypothesis . Remarkably, the estimated type‐1 error rates are even larger than 20% and 10%, respectively in Settings 6–10 (see Figure 1). In comparison, the permutation technique greatly improves the finite sample performance of all asymptotic procedures, and is therefore recommended in practical applications.

In order to investigate the impact of the underlying distributions of the follow‐up times, we resimulate the same designs with exponentially distributed follow‐up times . The results are displayed in the supplementary material. It can be seen that the shape of the underlying follow‐up times distributions slightly affect the behavior of the statistics in all scenarios. This is intuitively clear, since the different follow‐up times particularly influence the variance of the effect estimators, and increase the variance with wider ranging follow‐up times or certain amount of skewness. Therefore, all procedures tend to be slightly more liberal when wide ranging follow‐up times and small sample sizes are apparent. This can be particularly seen by the permutation test. The liberality, disappears with increasing sample sizes.

Power comparisons

The type‐1 error rate simulation results presented in Section 5.2.1 indicate a quite liberal behavior of the methods and ML‐based statistics under certain parameter constellations and small sample sizes. All methods tend to accurate conclusions with large sample sizes. The liberality of these methods increases the “power” of the methods to detect alternatives in small sample size settings. In an additional simulation study, not presented here, it turned out, that with large sample sizes, that is when all competing methods are accurate, their powers are all very similar.

Simulated coverage rates of the confidence intervals

Next we investigate the empirical coverage probabilities of the corresponding confidence intervals. Data were generated by and for varying , , and different overdispersion parameters. For illustration purposes, we only display the results using uniformly distributed follow‐up times, different overdispersion parameters and and rate . The results are displayed in Figure 2. It is readily seen that the competing procedures tend to be rather liberal, while the empirical coverage probabilities of the permutation‐based confidence intervals are closer to the nominal level of 95%. The quality of the approximation depends on sample sizes and the actual levels of heteroscedasticity across the groups and their allocations. If the larger sample has a smaller variance than the smaller sample (), the confidence intervals tend to be slightly liberal for small samples. However, this issue vanishes with increasing sample sizes.

Figure 2

Empirical coverage probabilities of nominal 95% confidence intervals of the corresponding confidence intervals given in (16), permutation‐ based confidence intervals given in (20) and ML‐based LRT statistics for different distributions and rate increments (x‐axis) and unequal overdispersion parameters (), where denote the realizations from

Simulation results for general metric data

As mentioned in the Introduction and in the description of model (1), data is not required to be count data and thus, numerical investigations of the behavior of the studentized permutation test are intriguing. We therefore investigate the empirical control of the type‐1 error rate of the studentized permutation test in completely heteroscedastic designs with metric data following exponential or χ2‐distributions. The method will be compared with using the standard normal approximation. Exponentially distributed variables were generated by , , and χ2‐variables were generated by , respectively.

The results are displayed in Table 2 and show that the studentized permutation approach controls the nominal type‐1 error rate very well and greatly improves the standard normal approximation.

Table 2

Type‐I error level (α = 5%) simulation results of the statistics in (18) and the permutation test in (19) using χ2‐square and exponentially distributed data in different designs, where denote the realizations from

		X_ik∼χt_ik2		X_ik∼Exp(t_ik·1/2)
n 1	n 2	T(L)(π,π^′)	T_(L)	T(L)(π,π^′)	T_(L)
7	7	0.0567	0.1168	0.0440	0.0937
7	15	0.0479	0.1063	0.0373	0.0884
12	12	0.0521	0.0862	0.0500	0.0807
12	20	0.0473	0.0825	0.0364	0.0644
17	17	0.0498	0.0757	0.0355	0.0565
17	25	0.0521	0.0769	0.0540	0.0822
27	27	0.0535	0.0694	0.0854	0.1058
27	35	0.0522	0.0684	0.0454	0.0618
32	32	0.0544	0.0698	0.0469	0.0575
32	40	0.0494	0.0634	0.0526	0.0623

TWO ILLUSTRATIVE EXAMPLES

Pediatric MS with disease onset under the age of 16 is uncommon and qualifies as a rare disease. Differences in clinical presentation before and after puberty have been reported (Huppke et al., 2014). Randomized controlled trials in pediatric MS have been very rare (Unkel et al., 2016), but are becomming more common now (Rose & Müller, 2016). We consider a randomized controlled trial assessing efficacy and safety of interferon beta‐1a compared to no treatment in pediatric MS reported by Pakdaman, Fallah, Sahraian, Pakdaman, and Meysamie (2006). In this trial, 16 patients were randomized to verum or control. Relapse rates and new T2 lesions were both considered as endpoints. The estimated rates and overdispersion parameters are given in Table 3. As a second example, we consider the Acyclovir trial reported by Lycke et al. (1996). In this experiment, Acyclovir treatment was used in a randomized, double‐blind, placebo‐controlled clinical trial with parallel groups to test the hypothesis that herpes virus infections are involved in the pathogenesis of MS. In total, adult patients were recruited, whereas were randomized to placebo or active treatment, respectively. The data (relapse counts) can be found in Figure 1 in the original publication (Lycke et al., 1996). As a secondary analysis of this trial, the relapse counts from patients that showed a progressive course during the trial were excluded from the statistical analysis. In this situation, patients have different follow‐up times and estimators must be weighted accordingly.

Table 3

Estimated rates and overdispersion parameters (Variance / Mean Ratio) for the two example studies

Endpoint	Group	Estimated rate λ^_i	Sample variance	Estimated overdispersion
Pediatric MS trial (N=16)
T2 lesions	Control	11.875	13.268	1.117
	Active	10.625	16.839	1.585
Relapses	Control	4.5	6.571	1.460
	Active	2.375	0.268	0.113
Acyclovir trial (N=60)
Relapses	Control	3.133	6.602	2.107
	ACYC	2.067	3.030	1.466
Acyclovir trial (N=60; Secondary analysis)
Relapses	Control	3.205	6.602	2.060
	ACYC	2.118	3.172	1.498

The estimated rates and overdispersions being defined as variance‐to‐mean ratios are given in Table 3. It can be readily seen from Table 3 that the overdispersion parameters seem to differ between the treatment groups, and even underdispersed counts are apparent. The effect of the different overdispersion parameters on the behavior of the statistical methods has been analyzed in detail in extensive simulation studies in Section 5.2.

Both motivating examples discussed above used over‐ and underdispersed counts as outcomes. Here, we present the results based on standard methods including normal approximation and maximum‐likelihood as well as the new developed methods. The test statistic being used is given by where denotes the estimated variance of the effect estimator using a MLE estimator of the overdispersion parameter ϕ, which is assumed to be identical across both treatment groups.

As competing methods, we also analyze the data using both a Negative Binomial Regression‐ and Poisson Regression using SAS PROC GENMOD.

Thus, the illustrative examples include constant as well as varying follow‐up times, and even the analyses with constant follow‐up times still presents a challenge since the sample sizes are with 16 and 60 very and moderately small, and the overdispersion is fairly pronounced, in particular for the MRI lesion counts and relapses. The effect estimates, standard errors, test statistics, P‐values as well as 95%‐confidence intervals are displayed in Table 4.

Table 4

Statistical analysis of the examples using : Approximate method, Effect (), Standard Error (SE), Test Statistic (= Effect / SE), and 95% confidence intervals

Method	Effect	SE	Statistic	P‐value	95% CI
T2 lesions
Normal (15)	0.111	0.174	0.638	0.524	(−0.231; 0.453)
LRT (21)	0.111	0.162	0.686	0.493	(−0.207; 0.429)
LRT.Pool (22)	0.111	0.161	0.691	0.489	(−0.204; 0.427)
Perm (19)	0.111	0.174	0.638	0.545	(−0.269; 0.510)
NB‐Reg	0.111	0.161	0.691	0.489	(−0.204; 0.428)
Pois‐Reg	0.111	0.149	0.745	0.456	(−0.181; 0.405)
Relapses
Normal (15)	0.639	0.216	2.964	0.003	(0.216; 1.062)
LRT (21)	0.639	0.302	2.116	0.034	(0.047; 1.231)
LRT.Pool (22)	0.639	0.284	2.254	0.024	(0.083; 1.195)
Perm (19)	0.639	0.216	2.964	0.026	(0.116; 1.162)
NB‐Reg	0.639	0.284	2.254	0.024	(0.096; 1.215)
Pois‐Reg	0.639	0.284	2.254	0.024	(0.096; 1.215)
Acyclovir relapses
Normal (15)	0.416	0.215	1.939	0.052	(−0.004; 0.837)
LRT (21)	0.416	0.228	1.824	0.068	(−0.031; 0.863)
LRT.Pool (22)	0.416	0.231	1.805	0.071	(−0.036; 0.868)
Perm (19)	0.416	0.215	1.939	0.054	(−0.007; 0.842)
NB‐Reg	0.416	0.231	1.805	0.071	(−0.035; 0.870)
Pois‐Reg	0.416	0.164	2.544	0.011	(0.098; 0.741)
Acyclovir relapses (Secondary analysis)
Normal (15)	0.414	0.218	1.904	0.057	(−0.012; 0.841)
LRT (21)	0.414	0.230	1.798	0.072	(−0.037; 0.866)
LRT.Pool (22)	0.414	0.233	1.781	0.075	(−0.076; 0.845)
Perm (19)	0.414	0.218	1.904	0.062	(−0.022; 0.845)
NB‐Reg	0.415	0.233	1.780	0.075	(−0.040; 0.874)
Pois‐Reg	0.422	0.165	2.553	0.011	(0.101; 0.750)

It can be readily seen from Table 4, that the estimated standard errors of the effect estimates for the T2 lesions are likely, and therefore all methods results in the same conclusion. Only the estimated standard error being computed via a Poisson‐Regression tends to be smaller. This occurs because the Poisson‐Regression sets the overdispersion to be zero, by default. A significant effect at 5% level can not be detected with any method (P > 0.05). The relapse rates are significantly different at 5%‐level of significance. It can be seen, however, that the estimates of the standard errors significantly differ from the moment‐based unbiased variance estimators (SE = 0.216 vs. SE = 0.302 using ML). Therefore, the P‐values based on ML estimates are larger than using the moments‐based estimator and standard normal distribution (P = 0.003 vs. P = 0.034). However, since sample size is rather small, the permutation approach is the most robust method in this setup, and results in a P‐value of P = 0.026. Since both over‐ and underdispersed counts were observed, the ML.Pool, the negative binomial, and poisson regression are tend to provide identical results.

The results obtained for the Acyclovir trial, however, differ significantly. First, both treatment groups show a different overdispersion. Therefore, the SE obtained by a Poisson‐Regression is way smaller than with all other methods, and thus results in a significant treatment effect at 5% level of significance. Comparing the other estimation approaches it can be seen that the ML‐based estimation approaches (assuming negative binomial distribution) of the SE tend to be larger than the unbiased methods‐of‐moments based methods. The largest SE is estimated via ML.Pool (which is identical to a NB‐Regression). The estimated standard error based on the unbiased variance estimate is given by SE = 0.215. Therefore, the P‐values range from 0.052 through 0.071. Due to the moderate sample size of , both the normal and permutation approximation tend to provide similar P‐values with and , respectively. The secondary analysis of the the Acyclovir trial shows similar results to the above. This occurs because only the relapse counts from four of the 60 patients were excluded from the analysis. However, slightly different effect estimates coming from the Negative Binomial and Poisson Regression can be seen. This occurs, because in case of unequal follow‐up times the rates are estimated using maximum likelihood estimation methods, which are not identical to moment (mean‐based) methods.

DISCUSSION

In this paper, inference methods for testing hypotheses formulated in terms of the effect rates of overdispersed counts were developed without assuming a specific data distribution and/or different overdispersion parameters. They are based on the asymptotic properties of novel unbiased estimators of the count rates and their variances. In order to provide valid methods for small sample sizes, resampling methods have been derived. Although data is in general not exchangeable, following the ideas of Neuhaus (1993), Janssen (1997, 2005), and Chung and Romano (2013), studentized permutation techniques could be applied. Simulation studies indicate, however, that the procedures control the nominal level reasonably well even with .

Furthermore, in clinical trials, the computation of confidence intervals for the treatment effects is important, following the ICH E9 guideline for randomized clinical trials: “Estimates of treatments shall be accompanied by confidence intervals, whenever possible& (ICH E9 Guideline 1998, chap. 5.5, p. 25). For instance, Saha (2013) investigates different methods for the computation of confidence intervals for the mean difference in the analysis of overdispersed count data (assuming constant follow‐up times ). In this paper, these procedures were generalized for possibly time‐varying and overdispersed count data and equipped with the studentized permutation approach. Extensive simulation studies show that the new methods improve the existing methods in terms of coverage probability and type‐I‐error rate control. Furthermore, we only considered one possible unbiased estimator of the rates by , which is known as a weighted mean estimator. Another unbiased estimator is given by the unweighted mean , or least‐square based estimators , where and denote the vectors of follow‐up times and response per group i, respectively. Investigating and comparing those estimators and generalizations thereof is tempting and will be subject to future research.

In future investigations, the results shall be extended to more general models allowing for covariates (e.g. for baseline adjustment) and several samples. Furthermore, investigating the overlap of range‐preserving confidence intervals for the effects is an interesting attempt for making inferences (Noguchi & Marmolejo‐Ramos, 2016).

CONFLICT OF INTEREST

The authors have declared no conflict of interest.

Supporting Information

Click here for additional data file.