The assumption of the symmetry of the underlying distribution is important to many statistical inference and modeling procedures. This paper provides a test of symmetry using kernel density estimation and the Kullback-Leibler information. Based on simulation studies, the new test procedure outperforms other tests of symmetry found in the literature, including the Runs Test of Symmetry. We illustrate our new procedure using real data.

Keywords: test of symmetry, power of the test, overlap coefficients, kernel density estimation, kullback-leibler information

Assume that a random sample, $X_1,X_2\mathrm{....},X_n$ , is drawn from absolutely continuous distribution having known median, assumed to be 0. In the case of an unknown median, or if the center of the distribution is not known, then the data can be centered by a consistent estimate of the median. However, the implications of centering the data around a consistent estimator of the median on the asymptotic properties are not straightforward. Therefore, further investigations are needed to study the robustness of the proposed test of symmetry and compare it with other available tests of symmetry when the median is unknown. In this paper we will discuss only the case where the median of the underlying distribution is assumed known.

Consider testing for symmetry $H_0:f(x)=f(-x)$ $\text{versus }{H}_a:f(x)\ne f(-x);\text{for some }x.$ Let $f_1(x)=f(x)\text{ and }{f}_2(x)=f(-x)$ . Under the null hypothesis, $D(f_1,f_2)=0$ . An equivalent hypothesis for testing the symmetry is $H_0:D(f_1,f_2)=0$ $\text{versus }{H}_a:D(f_1,f_2)>0$ let $\widehat{D}$ be a consistent nonparametric estimator of $D(f_1,f_2)$ . Under the null hypothesis of symmetry and some regularity assumptions, which will be discussed later in this paper, we propose the following test of symmetry:

$z_0=\frac{\widehat{D}-0}{{\widehat{\sigma }}_{\widehat{D}}}\stackrel{L}{\to }N(0,1)$

For large n, where ${\widehat{\sigma }}_{\widehat{D}}$ is a consistent estimator of the standard error of $\widehat{D}$ . An asymptotic significant test procedure at level $\alpha$ is to reject $H_0$ if $z_0>z_{\alpha }$ , where $z_{\alpha }$ is the upper $\alpha$ percentile of the standard normal distribution.

Kernel estimation of $D(f_1,f_2)$

For the i.i.d. sample $X_1,X_2\mathrm{....},X_n$ , let ${\widehat{D}}_{11}(f_1,f_1)$ be an estimate of $D_{11}(f_1,f_1)$ . To address which estimator of $D_{11}(f_1,f_1)$ will be appropriate to our inference procedure we need to state some necessary conditions: C1: f is continuous. (Smoothness conditions) C2: f is k times differentiable. (Smoothness conditions) C3: $D_{11}([X],[X])<1$ , where [X] is the integer part of X. (Tail condition) C4: $In{f}_{f(x)>0}f(x)>0$ (Tail condition) C5: ${\displaystyle \int f{(\ln f)}^2}<\infty$ (Peak condition) (Note that, this is also a mild tail condition.) C6: f is bounded. (Peak condition)

Some suggested estimators for $-D_{11}(f_1,f_1)={\displaystyle \underset{-\infty }{\overset{\infty }{\int }}{f}_1(x)\ln \left(f_1(x)\right)dx}$ may be found in the literature. These include the plug-in estimates of entropy which are based on a consistent density estimate $f_n$ of f. For example, the integral estimate of entropy introduced by Dmitriev & Tarasenko.¹⁴ Joe¹⁵ considers estimating $-D_{11}(f_1,f_1)$ when $f_1$ is a multivariate pdf, but he points out that the calculation when ${\widehat{f}}_1$ is a kernel estimator gets more difficult when the dimension of the integral is more than two. He therefore excludes the integral estimate from further study. The integral estimator can however be easily calculated if, for example, ${\widehat{f}}_1$ is a histogram.

The re-substitution estimate is proposed by Ahmad & Lin¹⁶ as follows:

$$-{\widehat{D}}_{11}({\widehat{f}}_1,{\widehat{f}}_1)=-\frac{1}{n}{\displaystyle \sum _{i=1}^n\ln {\widehat{f}}_1(X_i),}$$ (3)

Where ${\widehat{f}}_1$ is a kernel density estimator? They showed the mean square consistency of (3), such that ${}_n{\underrightarrow{\lim }}_{\infty }E\left\{{({\widehat{D}}_{11}({\widehat{f}}_1,{\widehat{f}}_1)-D_{11}(f_1,f_1))}^2\right\}=0$ Joe¹⁵ considers the estimation of $-D_{11}(f_1,f_1)$ for multivariate pdfs by an entropy estimate of the re-substitution type (3), also based on a kernel density estimate. He obtained asymptotic bias and variance terms, and showed that non-unimodal kernels satisfying certain conditions can reduce the mean square error. His analysis and simulations suggest that the sample size needed for good estimates increases rapidly when the dimension of the multivariate density increases. His results rely heavily on conditions C4 and C6. Hall & Morton¹⁷ investigated the properties of an estimator of the type (3) both when $f_n$ is a histogram density estimator and when it is a kernel estimator. For the histogram estimation they showed that ${}_n{\underrightarrow{\lim }}_{\infty }{n}^{1/2}({\widehat{D}}_{11}({\widehat{f}}_1,f_1)-D_{11}(f_1,f_1))\sim N(0,{\sigma }^2)$ under certain tail and smoothness conditions with ${\sigma }^2=Var(\ln (f(X))$ .(4)

Other estimators using sampling-spacing are investigated by Tarasenko,¹⁸ Beirlant & van Zuijlen,¹⁹ Hall,²⁰ Cressie,²¹ Dudewicz & van der Meulen,²² and Beirlant.²³ Finally, other nonparametric estimator has been discussed by many authors including Vasicek,²⁴ Dudewicz & Van der Meulen,²² Bowman²⁵ and Alizadeh.²⁶ Among these various entropy estimators, Vasicek’s sample entropy has been most widely used in developing entropy based statistical procedures. However, deriving the asymptotic distribution for there is hard to establish. Therefore, in this paper we will adopt the kernel re-substitution estimate which is proposed by Ahmad & Lin.¹⁶

We will adopt the notation of Samawi et al.⁹ Our proposed test of symmetry is as follow: Let $X_1,X_2\mathrm{....},X_n$ be a random sample from absolutely continuous distribution $F(.)$ which is continuously differentiable with uniformly bounded derivatives and having known median.

Let K be a kernel function satisfying the condition ${\displaystyle {\int }_{-\infty }^{\infty }K(x)dx=1}$

For simplicity, the kernel K will be assumed to be a symmetric density function with mean 0 and finite variance; an example is the standard normal density. The kernel estimators for $f(w_i)\text{ and }f(-w_i),i=1,2,\mathrm{...},C$ , are:

${\widehat{f}}_K(-w_i)=\frac{1}{nh}{\displaystyle \sum _{j=1}^{n}K\left(\frac{-w_i-x_j}{h}\right)}$ (6)

and

${\widehat{f}}_K(w_i)=\frac{1}{nh}{\displaystyle \sum _{j=1}^{n}K\left(\frac{w_i-x_j}{h}\right),}$

Respectively, where $C$ is the number of bins and depends on the sample size. As in Samawi et al. [9], we suggest to take the integer of $C=\sqrt{n}$ . In addition, $h$ is the bandwidths of the kernel estimators satisfying the conditions that $h>0,h\to 0\text{ and (}nh\to \infty )$ as $n\to \infty$ . There are many choices of the bandwidths ( $h$ ). In our procedure we use the method suggested by Silverman²⁷ Using the normal distribution as the parametric family, the bandwidths of the kernel estimators are $h=0.9A{(n)}^{-1/5}$ , (8)

Where $A$ =min{standard deviation of ( $x_1,x_2\mathrm{....},x_n$ ), interquantile range of ( $x_1,x_2\mathrm{....},x_n$ )/1.349}. This form of (8) is found to be adequate choices of the bandwidth for many purposes which minimizes the integrated mean squared error (IMSE),

$IMSE={\displaystyle \int E{[{\widehat{f}}_K(x)-f(x)]}^{2}dx.}$

We will use the Samawi et al.⁹ suggestion to calculate the bins as follows: Let $R=range(x_1,x_2,\mathrm{...},x_n)$ , then bins will be selected as $w_i=w_{i-1}+{\delta }_x,$ where $i=2,\mathrm{...},C$ , $w_1=\min (x_1,x_2,\mathrm{...},x_n)$ and ${\delta }_x=\frac{R}{C}$ .

Using the above kernel estimator, the nonparametric kernel estimator of $D(f_1,f_2)$ under the null hypothesis is given by

$\widehat{D}={\displaystyle \underset{}{\overset{}{\int }}{\widehat{f}}_K(x)\ln \left(\frac{{\widehat{f}}_K(x)}{{\widehat{f}}_K(-x)}\right)dx},\text{ = }{\widehat{D}}_{11}({\widehat{f}}_K(x),{\widehat{f}}_K(x))-{\widehat{D}}_{12}({\widehat{f}}_K(x),{\widehat{f}}_K(-x)),$

Which can be approximated by?

$\widehat{D}=\frac{1}{C}{\displaystyle \sum _{i=1}^C\ln {\widehat{f}}_K(w_i)-\frac{1}{C}{\displaystyle \sum _{i=1}^C\ln {\widehat{f}}_K(-w_i)}}$

The approximate variance of $\widehat{D}$ is given by

$Var(\widehat{D})=\frac{Var({\displaystyle \sum _{i=1}^C\ln {\widehat{f}}_K(w_i))}}{C^2}+\frac{Var({\displaystyle \sum _{i=1}^C\ln {\widehat{f}}_K(-w_i))}}{C^2}.$

Asymptotic properties of $\widehat{D}$

The nonparametric kernel estimator of $D(f_1,f_2)$ ( $\widehat{D}$ ) is based on the univariate kernel for density estimation, $K:ℝ\to ℝ$ . The necessary regularity conditions imposed on the univariate kernel for density estimation are:

I. ${\displaystyle {\int }_{R}K(z)dz=1.}$

II. ${\displaystyle {\int }_R{z}^{\beta }K(z)dz=0\text{ for any }\beta =1,\mathrm{...},r}-1\text{, and }{\displaystyle {\int }_R|z{|}^{r}K(z)dz<\infty .}$

III. $R={\displaystyle {\int }_R{K}^2(z)dz<\infty .}$

IV. $h>0,h\to 0\text{ , (}nh\to \infty )\text{ and (}\frac{nh}{\log n}\to \infty )$

These conditions may be found in Silverman²⁷ (Chapter 3) or Wand & Jones [28] (Chapter 2).

To show consistency of , apply the kernel density asymptotic properties found in Silverman,²⁷ (Chapter 3) or Wand & Jones,²⁸ (Chapter 2). Under assumptions 1-4 and assuming that the density $f:ℝ\to ℝ$ is continuous at each $w_i$ , i=1, 2,… C,

$$Bias({\widehat{f}}_K(-w_i))=o{(1)}_{-}\text{ and }Bias({\widehat{f}}_K(w_i))=o{(1)}_{+}$$ (12)

$$Var({\widehat{f}}_K(-w_i))=\frac{f(-w_i)}{nh}{\displaystyle {\int }_{ℝ}{K}^2(z)dz+o(\frac{1}{n{h}_{}}}\text{) and }Var({\widehat{f}}_K(w_i))=\frac{f(w_i)}{nh}{\displaystyle {\int }_{ℝ}{K}^2(z)dz+o(\frac{1}{nh}}\text{), }$$ (13)

and for $h>0,h\to 0\text{ and (}nh\to \infty )$ as $n\to \infty$

${\widehat{f}}_K(-w_i){\to }^{P}f(-w_i)\text{ and }{\widehat{f}}_K(w_i){\to }^{P}f(w_i)$ If f(.) uniformly continuous, then the kernel density estimate is strongly consistent. Moreover, as in Ahmad & Lin,¹⁶ ${}_C{\underrightarrow{\lim }}_{\infty }E\left\{{({\widehat{D}}_{11}({\widehat{f}}_K(x),{\widehat{f}}_K(x))-D_{11}(f_K(x),f_K(x)))}^2\right\}=0,$ and hence ${\widehat{D}}_{11}({\widehat{f}}_K(x),{\widehat{f}}_K(x))\stackrel{P}{\to }{D}_{11}(f_K(x),f_K(x)),\text{ as }C\to \infty$ and . However, since $\widehat{D}={\widehat{D}}_{11}({\widehat{f}}_K(x),{\widehat{f}}_K(x))-{\widehat{D}}_{12}({\widehat{f}}_K(x),{\widehat{f}}_K(-x))$ therefore $\widehat{D}\stackrel{p}{\to }D(f(w),f(-w)),\text{ as }C\to \infty .$

To drive the asymptotic distribution of $\widehat{D}$ , we will define $D(f_1,f_2)$ as a functional

$D(f_1,f_2)={\displaystyle \underset{-\infty }{\overset{\infty }{\int }}{f}_1(w)\ln (}{f}_1(w))dw-{\displaystyle \underset{-\infty }{\overset{\infty }{\int }}{f}_1(w)\ln }{f}_2(w)dw={\displaystyle \underset{-\infty }{\overset{\infty }{\int }}\ln (}{f}_1(w))d{F}_1-{\displaystyle \underset{-\infty }{\overset{\infty }{\int }}\ln }{f}_2(w)d{F}_1$

Using the previously stated regularity conditions, some regularity conditions given by Serfing²⁹ and assuming that the Gteuax derivatives of the functional $D(f_1,f_2)$ exist, we can show that the partial influence function of the functional $D(f_1,f_2)$ [30] are as follows:

$L_1(w;F_1,F_1)=\ln f_1(w)-{\displaystyle \underset{-\infty }{\overset{\infty }{\int }}{f}_1(w)\ln }{f}_1(w)dw,$

and

$L_2(w;F_1,F_2)=\ln f_2(w)-{\displaystyle \underset{-\infty }{\overset{\infty }{\int }}{f}_1(w)\ln }{f}_2(w)dw.$

Note that ${\displaystyle \int L_1(w;F_1(w),F_1(w))}d{F}_1(w)=0\text{ and }{\displaystyle \int L_2(w;F_1(w),F_2(w))}d{F}_1(w)=0.$ Now using this functional representation of $D(f_1,f_2)$ , then as in Samawi et al.³⁰ and Serfing,²⁹

$\sqrt{C}(\widehat{D}-D(f_1,f_2))\stackrel{L}{\to }N(0,{\sigma }_{\widehat{D}}^2),$

where ${\sigma }_{\widehat{D}}^2={\displaystyle \int L_1^2(w;F_1,F_1)d{F}_1}+{\displaystyle \int L_2^2(w;F_1,F_2)d{F}_1}$

A consistent estimate for ${\sigma }_{\widehat{D}}^2$ is given by ${\widehat{\sigma }}_{\widehat{D}}^2=\frac{1}{C}{\displaystyle \sum _{i=1}^C{L}_1^2(w;{\widehat{F}}_1,{\widehat{F}}_1)}+\frac{1}{C}{\displaystyle \sum _{i=1}^C{L}_2^2(w;{\widehat{F}}_1,{\widehat{F}}_2)},$

Where, $L_1^{}(w_i;{\widehat{F}}_1,{\widehat{F}}_1)=\ln {\widehat{f}}_1(w_i)-{\widehat{D}}_{11}({\widehat{f}}_1(w_i),{\widehat{f}}_1(w_i))\text{ and }{L}_2^{}(w_i;{\widehat{F}}_1,{\widehat{F}}_2)=\ln {\widehat{f}}_2(w_i)-{\widehat{D}}_{12}({\widehat{f}}_1(w_i),{\widehat{f}}_2(w_i)),i=1,2,\mathrm{...},C,$

Where in our case $f_1(w_i)=f(w_i)\text{ and }{f}_2(w_i)=f(-w_i)$ .

For discussions about different methods addressing the issue of the performance of kernel density estimation at the boundary, see Hall & Park.³¹

As in Samawi et al.,⁹ to gain some insight of our procedure, a simulation study was conducted to investigate the performance of our new test of symmetry based on $\widehat{D}$ . We compared our proposed test of symmetry with the test proposed by McWilliams,⁵ Modarres & Gastwirth,³² Mira⁸ Bonferroni’s test, and Samawi et al.⁹ tests of symmetry.

As in McWilliams [5], the runs test is described as follows: For any random sample of size n, let $Y_{(1)},Y_{(2),\mathrm{...},}{Y}_{(n)}$ denote the sample values ordered from the smallest to largest according to their absolute value (signs are retained), and $S_1,S_{2,\mathrm{...},}{S}_n$ denote indicator variables designating the sign of the $Y_{(j)}$ values [ $S_j=1\text{ if }{Y}_{(j)}\text{ is nonnegative, 0 otherwise}$ ]. Thus, the test statistic used for testing symmetry is $R_{}^{*}$ = the number of runs in $S_1,S_{2,\mathrm{...},}{S}_n$ sequence= $1+{\displaystyle \sum _{j=2}^n{I}_j}$ , where

$I_j=\{\begin{array}{c}0\text{ if }{S}_j=S_{j-1}\\ 1\text{ if }{S}_j\ne S_{j-1}\end{array}$

We reject the null hypothesis if $R_{}^{*}$ is smaller than a critical value ( $c_{\alpha }$ ) at the pre-specified value of $\alpha$ . Moreover, Mira [8] Bonferroni’s test is ${\gamma }_1(F_n)=2({\overline{X}}_n-X_{s:n})$ , where $X_{s:n}=Median(X_1,X_2,\mathrm{...},X_n)$ . The process is to reject the null hypothesis if $|{\gamma }_1(F_n)|\ge \frac{a_n}{\sqrt{n}}{S}_c({\gamma }_1,F_n),$ where

$\begin{array}{l}{a}_n\to z_{1-\frac{\alpha }{2}}\text{ as }n\to \infty ,S_c^2({\gamma }_1,F_n)=4{\hat{\sigma }}^2+-4D_{n,c}{S}_{{\overline{\mu }}_F},\\ {\hat{\sigma }}^2=\frac{1}{n-1}\sum _{i=1}^n(X_i-{\overline{X}}_n{)}^2,S_{{\overline{\mu }}_F}={\overline{X}}_n-\frac{2}{n}\sum _{i=1}^n{X}_{i}I(X_i\le X_{s:n}),D_{n,c}=\\ \frac{n^{1/5}}{2c}(X_{\left[\right(n/2)+c{n}^{4/5}]:n}-X_{\left[\right(n/2)+c{n}^{4/5}+1]:n}),\text{and }c=\mathrm{0.5.}\end{array}$

The Modarres & Gastwirth³² test is the hybrid test of sign test in the first stage and a percentile-modified two-sample Wilcoxon see Gastwirth³³ test in the second stage. Finally, Samawi et al.⁹ test of symmetry is based on kernel estimate of the overlap measure.

In the following simulation, SAS version 9.3 {proc kde; method=srot} is used. As in McWilliams,⁵ the generalized lambda distribution see, Ramberg & Schmeiser³⁴ is used in our simulation with following set of parameters:

To generate the observations we used $x_i={\lambda }_1+\frac{1}{{\lambda }_2}(u_i^{{\lambda }_3}-{(1-u_i)}^{{\lambda }_4},i=1,\mathrm{...},m,$ where $u_i$ a uniform random number. The significance level used in the simulation is $\alpha =0.05,$ with sample sizes n=30, 50, and 100. To investigate the Type I error, the symmetric distributions used in the simulation are the first case of the generalized lambda and the normal. Our simulation is based on 5000 simulated samples. The 95% confidence intervals of the true probability of type I error under the null hypothesis with $\alpha =0.05$  are (0.04396, 0.05504).

Table 1.1 shows the estimated probability of type I error. Our test is an asymptotic test with a slight bias in D(., .) and in the variance estimation for small sample size. For sample sizes more than 30, the test seems to have an estimated probability of type I error close to the nominal value 0.05. However, Bonferroni’s test seems to be conservative test procedure, while Modarres, Gastwirth test is slightly conservative for small sample size. Table 1.2 and Table 1.3 show that using D(., .) based test is more powerful than McWilliams,⁵ Bonferroni’s, Modarres & Gastwirth³² and Samawi et al.⁹ tests in all of the presented cases. The efficiency increases as the sample size increases.

Note: The values of skewness $({\alpha }_{\text{3}}\text{)}$ and kurtosis $({\alpha }_4)$ are from McWilliams.⁵

We applied our new test procedure of symmetry to the base deficit (bd) data as in Samawi et al.⁹ The base deficit score refers to a deficit of "base" present in the blood. Base deficit scores were first established by Davis et al.³⁵ The base deficit score has been found correlated to many variables in the trauma population, such as, mechanism of injury, the presence of intra-abdominal injury, transfusion requirements, mortality, the risk of complications, and the number of days spent in the intensive care unit as indicated by Tremblay et al.³⁶ and Davis et al.³⁷   

The samples used in this illustration are part from the data collected based on a retrospective study of the trauma registry at a level 1 trauma center between January, 1998 and May, 2000. The primary concern was to determine at what point we can differentiate between life and death based on a base deficit score. A first step in this analysis is to determine if there is a difference in location for the base deficit score of those who survive and those who fail to survive. As is frequently the case in such studies, the underlying distribution is assumed “normal” or at least symmetric and a t-test or a nonparametric test would be performed without checking the assumptions. In either case a test of symmetry is almost never considered as a means of determining how one may proceed in the analysis. Based on the conclusions of a test of symmetry, the analyst can chose the most powerful test for location. The goal is to test the hypothesis that, on average, the base deficit score is the same for those who survive and those who fail to survive their injuries. The injuries of interest in this group of patients are either penetrating injury or blunt injury. However, before deciding on the test procedure, we need to check the assumptions of underlying distribution of the base deficit for both penetrating injury and blunt injury groups of patients. In particular, the assumption of symmetry of the underlying distribution needs to be verified. The data will be centered about the estimated measure of location to perform the tests of symmetry.

Figure 1.1 and Figure 1.2 show the box plot for penetrating injury and blunt injury groups for dead and alive patients respectively. Clearly there is some asymmetry on all four distributions. Also, Table 2.1 and Table 2.2 show summery statistics for penetrating injury and blunt injury groups for dead and alive patients respectively. Table 2.3 shows the overlap based test, the runs test and the proposed test of symmetry based on the Kullback-Leibler information of symmetry for the underlying distribution for patients discharged alive and dead patients of blunt trauma and penetrating trauma. We reject the assumption of symmetry for underlying distribution of these groups.

Figure 1.1 Box plot to base deficit for dead patients.

Figure 1.2 Box plot to base deficit for alive patients.

The proposed test of symmetry based on the Kullback-Leibler information, appears to outperform the other tests of symmetry in the literature in terms of power. Our test is more sensitive to detect a slight asymmetry in the underlying distribution than other tests proposed in the literature. Moreover, the kernel density estimation literature is very rich and many of the proposed methods and the improved methods are available on statistical software, such as SAS™, S-plus, Stata and R. Since based on the Kullback-Leibler information can be used in multivariate cases as well as in univariate cases, our proposed test of symmetry can be extended to multivariate cases for diagonal symmetry, conditional symmetry and other types of symmetry.