Suggestion of confidence interval methods for the Cronbach alpha in application to complex survey data
Section 2. Design-based confidence intervals for $\alpha$

• Previous
• Next

In this section, we discuss two methods to obtain the confidence interval for $\alpha ,$ the confidence interval based on the linearization method using the influence function (Deville, 1999; Demnati and Rao, 2004) and the coverage-corrected bootstrap method (Hall, Martin and Schucany, 1989). In this discussion, we consider strategies to deal with stratification, since stratification is a common feature in surveys and may decrease the magnitude of the variances for the statistics of interest (Lohr, 1999). We note that the sampling design for the NCS-R used stratification (more details in Section 4). Later, in Section 3, we show that the linearization will be sufficient for most practical cases (e.g., scales with ordinal responses); however, the coverage rate may not be satisfactory with some non-normal distributions. The coverage-corrected bootstrap method when applied to survey data is proposed as a possible alternative to the linearization method in those cases (Section 2.2).

2.1  Linearization

A symmetric confidence interval can be obtained based on the normal approximation of an estimator for a finite population (Hájek, 1981; Sen, 1995). The linearization method is applied for the variance estimation of complex statistics. In a survey sampling setting, we consider a population index set $U=\left\{1,\dots ,N\right\}$ with population size $N.$ A random sample $S$ of size $n$ is selected from $U$ by a sampling design $p\left(s\right)=\mathrm{Pr}\left\{S=s\right\}$ for all $s\subset U.$ The value $w_k$ denotes the sampling weight associated with the index $k\in s.$ For probability sampling, the sampling weight for index $k$ is the inverse of the first order inclusion probability, i.e., $w_k={\left[\mathrm{Pr}\left\{k\in s\right\}\right]}^{-1}.$ For each unit $k$ of the population $U,$ there is a point (or observation) $x_k$ of $R^p,$ a $p$ -dimensional real space. In a similar manner to Deville (1999), let us consider the population $U$ that is represented by the measure $M$ as having a mass of $1/N$ in each of the points $x_k.$ In this way, we have ${\displaystyle \int 1dM}=1$ and ${\displaystyle \int ydM}=N^{-1}{\displaystyle {\sum }_{k\in U}{y}_k}$ for any vector value $y_k=y\left(x_k\right),$ where we define the integral of a vector as the integral of each component of the vector. The measure $\widehat{M}$ is the estimator of $M$ allocating a weight $w_k/N$ to any point $x_k,$ $k\in s$ and 0 to any other points. Following some conventional notation (e.g., Cochran, 1977), let ${\displaystyle \int ydM}=\overline{Y}.$ Also let ${\displaystyle \int yd\widehat{M}}=\widehat{\overline{Y}}.$ The influence function of a “functional” $T$ is defined as

$$IT\left(M;x\right)=\underset{t\to 0}{\lim }\frac{T\left(M+t{\delta }_x\right)-T\left(M\right)}{t},$$

where ${\delta }_x$ denotes the added unit mass at point $x$ (Deville, 1999), and the functional $T$ (Krätschmer, Schied and Zähle, 2012) maps a measure to a set (e.g., the real line). The examples of the functional include $\overline{Y}$ and $\widehat{\overline{Y}}.$ Note that this classical definition of the influence function (Hampel, Ronchetti, Rousseeuw and Stahel, 1986; Davison and Hinkley, 1997) is slightly different from that of Deville (1999) where he defines a measure $M$ to satisfy ${\displaystyle \int ydM}={\displaystyle {\sum }_{k\in U}{y}_k}.$ Let us define the linearized value $z_k=IT\left(M;x_k\right).$ Let $T\left(\widehat{M}\right)$ indicate the substitution estimator of $T\left(M\right)$ by replacing $M$ by $\widehat{M}.$ Assume that the postulate of Deville (1999), i.e., $n^{-1/2}{N}^{-1}\left(\widehat{X}-X\right)$ has a zero-mean multi-normal distribution as a limit, where $X$ and $\widehat{X}$ are the population total and the total estimator for general observation $x_k,$ and $N$ and $n$ tend toward infinity. This fact leads to ${\displaystyle \int xd\widehat{M}}-{\displaystyle \int xdM}=O_p\left(n^{-1/2}\right).$ Assuming that $T$ can be derived for any direction of an increase, a similar argument to Deville (1999) gives rise to the result

$$\frac{T\left(\widehat{M}\right)-T\left(M\right)}{N^c}=\frac{1}{N^{c+1}}{\displaystyle \sum _{k\in U}{z}_k\left(w_k-1\right)+o\left(n^{-1/2}\right)},(2.1)$$

for some positive value $c.$ Equation (2.1) results in the asymptotic variance of $T\left(\widehat{M}\right)$

$$\text{Avar}\left\{T\left(\widehat{M}\right)\right\}\text{=}\mathrm{var}\left(\widehat{\overline{Z}}\right).(2.2)$$

If $T={\displaystyle \int xdM}\left(x\right),$ then the influence function at $x_k$ $\left(k\in U\right)$ is

$$IT\left(M;x_k\right)=\underset{t\to 0}{\lim }\frac{{\displaystyle {\sum }_{i\in U}{x}_i/N+t{x}_k-{\displaystyle {\sum }_{i\in U}{x}_i/N}}}{t}=x_k.(2.3)$$

For a complex statistic as the functions of simple statistics, we have the influence function

$$I\left(f\left(T\right)\right)=D\left(f\right)IT,(2.4)$$

where $f$ is a differentiable function on the space of values for $T$ and $D\left(f\right)$ is the matrix of the partial derivatives of $f$ (Deville, 1999). In many cases, the linearized value $z_k$ includes parameters to be estimated. Let ${\widehat{z}}_k$ indicate the approximation of $z_k$ using some statistics estimated by the sample. Deville (1999) notes that with a fixed and finite number of estimated parameters, the variance estimators based on ${\widehat{z}}_k$ and $z_k$ are equivalent by an asymptotically negligible quantity.

Now, we obtain the linearized value for $\alpha$ as follows. Consider a data set

$$X={\left(x_1,\dots ,x_n\right)}^T\text{},$$

where $x_k$ is a $p$ -variate observation indicating $p$ items in an instrument and $n$ is the sample size. Let ${\sigma }_{ij}$ and ${\widehat{\sigma }}_{ij}$ $\left(i,j=1,\dots ,p\right)$ denote the ${\left(i,j\right)}^{\text{th}}$ elements of $\Sigma$ and $\widehat{\Sigma }$ as defined in Section 1, respectively. Specifically, we define ${\sigma }_{ij}={\displaystyle {\sum }_{k=1}^N\left({}_{i}x{}_k-{}_i\overline{X}{}_k\right)\left({}_{j}x{}_k-{}_j\overline{X}{}_k\right)/\left(N-1\right)}$ (Lohr, 1999), where ${}_{i}x{}_k$ and ${}_i\overline{X}{}_k$ are $i^{\text{th}}$ element of $x_k$ and its population mean, respectively. For simple random sampling without replacement (SRSWOR), we define ${\widehat{\sigma }}_{ij}={\displaystyle {\sum }_{k\in s}\left({}_{i}x{}_k-{}_i\widehat{\overline{X}}{}_k\right)\left({}_{j}x{}_k-{}_j\widehat{\overline{X}}{}_k\right)/\left(n-1\right)},$ where $n$ is the size of the sample $s$ and ${}_i\widehat{\overline{X}}{}_k$ is the sample mean of ${}_{i}x{}_k$ (Lohr, 1999). For obtaining ${\widehat{\sigma }}_{ij}$ for more complicated sampling methods including unequal probability sampling, we refer to Swain and Mishra (1994) and Patel and Bhatt (2016). In survey sampling, the sampling weights are used for correcting the disproportionality of the sample regarding the target population of interest (Pfeffermann, 1993). With the complex sampling designs often used in practice, failure to consider the sampling designs may provide biased inferences. For more of a discussion of the role of sampling weights, we refer to Pfeffermann (1993). For the variance estimation of survey data, the linearization method can be applied as in formula (2.2) incorporating the sampling weights. Following conventional notations of the vectorization of a matrix, let $\text{vech}\left(A\right)$ be the column vector of nonduplicated elements of the matrix $A,$ $vec\left(A\right)$ be the column vector composed of the columns of $A.$ Let $\widehat{t}$ indicate the collection of statistics as the components of ${\text{vech}}^T\left(\widehat{\Sigma }\right)$ and $t$ indicate the collection of corresponding parameters. Specifically, we let $t=\left({\text{vech}}^T\left({\displaystyle {\sum }_{k\in U}{x}_k{x}_k^T/N}\right),{\left({\displaystyle {\sum }_{k\in U}{x}_k/N}\right)}^T\right).$ Also, let the matrix $K_p$ indicate a transition matrix that satisfies the relationship $\text{vech}\left(A\right)=K_p^T\text{vec}\left(A\right),$ which borrows the transition matrix expression from van Zyl et al. (2000). We propose a linearized value for $\text{Var}\left(\widehat{\alpha }\right)$ as

$$z_k=\frac{p}{p-1}\frac{1}{{\left(1^T\widehat{\Sigma }1\right)}^2}\left\{1^T\widehat{\Sigma }1{\text{vec}}^T\left(I_p\right)-\text{tr}\left(\widehat{\Sigma }\right)\left(2{\text{vec}}^T\left(11^T\right)-{\text{vec}}^T\left(I_p\right)\right)\right\}{K}_p\widehat{J}{u}_k,k=1,\dots ,n,(2.5)$$

where a Jacobean matrix $J=\partial {\text{vech}}^T\left(\Sigma \right)/\partial t,$ $\widehat{J}=\left[\partial {\text{vech}}^T\left(\Sigma \right)/\partial t\right]{}_{\Sigma =\widehat{\Sigma },t=\widehat{t}},$ $u_k={\left({\text{vech}}^T\left(x_k{x}_k^T\right),x_k^T\right)}^T$ and $I_p$ is the $p\times p$ identity matrix. We can now obtain the linearized value (2.5).

Derivation of (2.5): We consider the variance of $p/\left(p-1\right)\text{tr}\left(\widehat{\Sigma }\right)/1^T\widehat{\Sigma }1$ since its variance is the same as $\text{Var}\left(\widehat{\alpha }\right).$ Let ${\alpha }^{*}=p/\left(p-1\right)\text{tr}\left(\Sigma \right)/1^T\Sigma 1.$ Also, let ${\text{vec}}^T\left(\Sigma \right)=\left({\sigma }_{11},\dots ,{\sigma }_{p1},{\sigma }_{12},\dots ,{\sigma }_{p2},{\sigma }_{13},\dots ,{\sigma }_{pp}\right)$ and ${\text{vech}}^T\left(\Sigma \right)=\left({\sigma }_{11},\dots ,{\sigma }_{p1},{\sigma }_{22},\dots ,{\sigma }_{p2},{\sigma }_{33},\dots ,{\sigma }_{pp}\right),$ a $p^2$ -vector and a $p\left(p+1\right)/2$ -vector, respectively. Then, we have

$$\frac{\partial {\alpha }^{*}}{\partial t}=\frac{\partial {\alpha }^{*}}{\partial {\text{vech}}^T\left(\Sigma \right)}J=\frac{p}{p-1}\frac{\partial }{\partial {\text{vech}}^T\left(\Sigma \right)}\left\{\frac{\text{tr}\left(\Sigma \right)}{1^T\Sigma 1}\right\}J.(2.6)$$

Now, in (2.6), we can show

$$\begin{array}{ll}\frac{\partial }{\partial {\text{vech}}^T\left(\Sigma \right)}\left\{\frac{\text{tr}\left(\Sigma \right)}{1^T\Sigma 1}\right\}\hfill & =\frac{1}{1^T\Sigma 1}\frac{\partial \text{tr}\left(\Sigma \right)}{\partial {\text{vech}}^T\left(\Sigma \right)}-\frac{\text{tr}\left(\Sigma \right)}{{\left(1^T\Sigma 1\right)}^2}\frac{\partial 1^T\Sigma 1}{\partial {\text{vech}}^T\left(\Sigma \right)}\hfill \\ \hfill & =\frac{1}{{\left(1^T\Sigma 1\right)}^2}\left\{1^T\Sigma 1{\text{vec}}^T\left(I_p\right)-\text{tr}\left(\Sigma \right)\left(2{\text{vec}}^T\left(11^T\right)-{vec}^T\left(I_p\right)\right)\right\}{K}_p.(2.7)\hfill \end{array}$$

Using (2.6) and (2.7), we can obtain

$$\frac{\partial {\alpha }^{*}}{\partial t}=\frac{p}{p-1}\frac{1}{{\left(1^T\Sigma 1\right)}^2}\left\{1^T\Sigma 1{\text{vec}}^T\left(I_p\right)-\text{tr}\left(\Sigma \right)\left(2{\text{vec}}^T\left(11^T\right)-{\text{vec}}^T\left(I_p\right)\right)\right\}{K}_{p}J.(2.8)$$

Note that the expression (2.8) is a vector that consists of the derivatives of ${\alpha }^{*}$ with respect to the components of $t.$ Each element in (2.8) is multiplied by the influence function corresponding to the statistics in $t$ as in (2.4). This is accomplished by multiplying (2.8) by $u_k={\left({\text{vech}}^T\left(x_k{x}_k^T\right),x_k^T\right)}^T\text{},k=1,\dots ,n,$ which is obtained by using (2.3). Now substituting $\Sigma$ by $\widehat{\Sigma }$ leads to the linearized value (2.5).

The formula for the new value (2.5) is easily implemented in the computer code using commonly available computer software. The relevant R code is available in the Supplementary Material.

We note that, in application to survey sampling, the estimate $\widehat{\Sigma }$ should be obtained properly by incorporating the survey design. The variance is estimated by $\widehat{\text{Var}}\left(\widehat{\overline{Z}}\right),$ where $\widehat{\text{Var}}$ indicates an operation to obtain the variance incorporating the weights and survey design properly, e.g., the Sen-Yates-Grundy variance estimator (Sen, 1953; Yates and Grundy, 1953), an unbiased variance estimator for the Horvitz-Thompson estimator (Horvitz and Thompson, 1952) under designs with fixed sample sizes (e.g., Särndal, Swensson and Wretman, 1992) or the variance estimator for sampling with the replacement as a conservative approximation (Wolter, 1985). Specifically, in this paper, the variance for the NCS-R data is estimated as

$$\widehat{\text{Var}}\left(\widehat{\overline{Z}}\right)={\displaystyle \sum _{h=1}^H\widehat{{\text{Var}}_h}\left(\widehat{\overline{Z}}\right)},(2.9)$$

where $\widehat{{\text{Var}}_h}\left(\widehat{\overline{Z}}\right)$ indicates the design-specific variance estimator for stratum $h\left(h=1,\dots ,H\right).$ Once $z_k$ values are obtained, standard statistical software for survey sampling such as R package “survey” (Lumley, 2004) can be used for the calculation of (2.9).

Now, consider a case that $x$ is a random variable following a distribution and that an observation is a realization of the random variable; in addition, a sample of size $n$ is obtained according to the random variable. In this specific case, we do not consider the finite population, where the design-based variance estimation is suitable as shown in the previous discussion. In a random variable setting, let $\widehat{\overline{Z}}$ indicate the estimator with the measure $\widehat{M}$ as the empirical distribution function (Fernholz, 1991). Employing the concept of a robust statistical inference based on the influence function (Davison and Hinkley, 1997), the sample variance for the population can be calculated by

$$\widehat{\text{Var}}(\widehat{\overline{Z}})=n^{-1}{\displaystyle \sum _{k=1}^n{(z_k-\overline{z})}^2}/(n-1),(2.10)$$

where $z_k$ is the linearized value (2.5) obtained from the statistic $\widehat{\alpha }$ based on the sample (size n) and $\overline{z}$ is the sample mean of $z_k(k=1,\dots ,n).$ We also note that the formula (2.10) is not constructed for infinite populations in survey methodology, where the finite population is seen as a realization from an infinite population. In that case, the outcomes of a statistical model give rise to the values of the characteristics of interest in the finite population, thus the model-based variance estimation is appropriate (Binder and Roberts, 2009). The formula (2.10) can be used for a general data analytical setting, where observations are considered as realizations of a random variable.

2.2  The coverage-corrected bootstrap method

The linearization provides reasonable estimates for the confidence intervals; however, in some cases, the coverage rate may not be satisfactory when the underlying distributions are non-normal (see Section 3). In these cases, some computer-intensive approaches such as the double bootstrap method, which is also called the coverage-corrected bootstrap may be implemented (Hall et al., 1989). We primarily discuss the double bootstrap method instead of the typical “single” bootstrap method (DiCiccio and Romano, 1988) since we observe that the single bootstrap method may not be satisfactory with non-normal underlying distributions (e.g., lognormal distribution) in terms of the coverage rate (Table 3.3).

For adjusting the bootstrap weight, the rescaling method referred to as the Rao-Wu bootstrap (Rao and Wu, 1988) is a popular approach for analyzing a lot of survey data, e.g., from Statistics Canada surveys (Mach, Saïdi and Pettapiece, 2007). The Rao-Wu bootstrap method is based on the assumption of sampling with a replacement, but is often employed for sampling without a replacement as well, when the first-stage sampling fraction is negligible (Mach et al., 2007). Herein, we propose implementing the coverage-corrected bootstrap method using the weight adjustment from Rao and Wu (1988). Among the various bootstrap confidence interval techniques (e.g., for these varieties, see Hwang, 1995), we consider the percentile bootstrap interval, which is a strictly nonparametric bootstrap approach (Hall, Martin and Schucany, 1989).

The coverage rates of the bootstrap confidence intervals can be corrected by incorporating additional bootstrap procedures. Because of bootstrapping the bootstrap sample, this kind of a procedure is referred to as the double bootstrap method (Martin, 1992). It is known that this method reduces the coverage error of two-sided confidence intervals by a factor of the order $n^{-1}$ compared to the single bootstrap or normal-theory confidence intervals (Martin, 1992). Suppose $\widehat{l}$ and $\widehat{u}$ are the lower and upper bounds of the percentile bootstrap confidence interval using the original data. As proposed by Hall et al. (1989), the $100\left(1-q\right)\%$ coverage-corrected bootstrap confidence interval can be defined as $\left(\widehat{l}-\delta ,\widehat{u}+\delta \right),$ where a positive value of $\delta$ satisfies

$$1-q=\mathrm{Pr}\left\{\widehat{\alpha }\in \left({\widehat{l}}^{*}-\delta ,{\widehat{u}}^{*}+\delta \right)\right\}.(2.11)$$

The values ${\widehat{l}}^{*}$ and ${\widehat{u}}^{*}$ indicate the lower and upper bounds, respectively, of the confidence interval obtained by bootstrapping a resampled data set. The probability in the right-hand side of equation (2.11) is empirically evaluated as shown in the following steps.

We use $\left(\max \left(\widehat{l}-\delta ,0\right),\min \left(\widehat{u}+\delta ,1\right)\right)$ since the true $\alpha$ is assumed to be between 0 and 1.

In the analysis, we have to resample the data without disrupting the survey design structure. The bootstrap is carried out within each stratum, and all observations in the same cluster should be kept together in a resampled data set (Lohr, 1999). For each resampled data set, new weights need to be obtained (Rao, Wu and Yue, 1992). Specifically, let $n_h$ indicate the sample size of the primary sampling unit (PSU) in stratum $h\left(h=1,\dots ,H\right).$ Suppose we resample $n_h^{*}$ clusters for each stratum. Then, the rescaled weight for observation $k$ in the resample is

$$w_k^{\left(b\right)}=w_k\left\{\left(1-\sqrt{\frac{n_h^{*}}{n_h-1}}\right)+\sqrt{\frac{n_h^{*}}{n_h-1}}\frac{n_h}{n_h^{*}}{m}_k\right\},(2.12)$$

where $m_k$ is the number of repetitions of the PSU that observation $k$ belongs to and $w_k$ is the original weight of observation $k$ (Rao et al., 1992; Mach, Dumais and Robinson, 2005; Mach et al., 2007). When $n_h^{*}=n_h-1,$ the bootstrap weight becomes $w_k^{\left(b\right)}=w_k\left\{\frac{n_h}{n_h-1}{m}_k\right\},$ which is a conventional bootstrap weight (Lohr, 1999). This procedure is repeated to obtain a total of $B$ bootstrap samples. For the actual data analysis, we use $B=500$ following the common practice of Statistics Canada surveys (Canadian Community Health Survey - Annual Component, 2007). To obtain the estimates, the stratification or cluster structure is no longer considered since the bootstrap weights take into account the survey design structure (Lohr, 1999). The percentile interval will be obtained based on the $B$ values of the estimates of $\alpha .$ For each resample, $\alpha$ is estimated based on the sample variance and covariance matrix incorporating the weights. To obtain the coverage-corrected confidence interval, we carry out the additional bootstrap with each bootstrap sample in a similar manner to what was explained above. In the simulation and data analysis we use 200 bootstrap samples for the second round of bootstrapping. The relevant R code is provided in the Supplementary Material.

• Previous
• Article main page
• Next