Understanding the Factors Influencing the Number of Extracurricular Clubs in American High Schools

2.1 Sampling of Schools

We first present the characteristics of schools, followed by a description of our school sampling methodology. Figure 1 shows the cumulative distribution of school count and student count as a function of school size. The “all schools” curves represent the 21,064 schools. The “sampled schools” curves represent 229 schools that were randomly selected from the 21,064 schools for this study. We will describe the selection of these schools later.

Figure 1

Cumulative distribution of school count and student count as functions of school size. While most schools are small, most students attend larger schools. To avoid a uniform sample dominated by small schools, we intentionally oversampled larger ones.

As an example of how to read the figure, the data point in the green circle on the “by school count (all schools)” curve means that 50% (see the Y-axis) of all schools have 357 (see the X-axis) or fewer students. Similarly, the data point in the green square on the “by student count (all schools)” curve means that 50% of all students across all schools are enrolled in schools that each have 1,434 or fewer students.

These data highlight a significant divergence between school count and student count. While the majority of schools are small, the majority of students are enrolled in larger schools. This affects how we select the sampled schools. Initially, we selected the sampled schools in a way that ensured, in Figure 1, the curve of “by school count (sampled schools)” closely matched the curve of “by school count (all schools).” However, this led to a lack of sufficient data points for large schools because small schools make up the majority. Specifically, schools with 1,000 or fewer students account for 75% of all schools while hosting only 34% of all students. To address this issue, we intentionally sampled more schools of larger sizes. This leads to the divergence of the “all schools” curves and the “sampled schools” curves in the figure. We will delve further into this topic during the data analysis.

2.2 Collecting Club Data

The NCES data contain no information about extracurricular clubs offered by schools. To collect such data, traditional approaches would typically resort to surveys to contact individuals at specific schools, with approvals from numerous ethics committees. However, this approach is challenging to scale to hundreds of schools nationwide, which might be a key reason why, after more than half a century of research on extracurricular activities, there are still no publicly available large-scale datasets providing complete lists of clubs offered by many schools.

In the Internet era, we instead take an innovative approach to collect club data from schools’ public websites. Interestingly, the shift to online learning during the COVID-19 pandemic facilitated this data collection process, as an increasing number of schools have adopted websites as a primary channel for communication.

Specifically, we sampled thousands of schools, searched their respective websites, and identified 229 schools with complete lists of clubs they offer. It is important to note that while many school websites include information about a limited subset of clubs as examples, they do not present a comprehensive list of all available clubs. We have excluded these schools from our study. On average, it takes searching through more than 15 schools to identify one with sufficiently complete club data. Consequently, we conducted searches across thousands of schools to identify the 229 sampled schools. In total, the 229 schools hosted 221,300 students and offered 5,983 clubs.

Among a school’s clubs, this study excludes clubs for varsity sports such as swimming and baseball for two reasons: (1) They form a substantial category on their own and merit a separate, dedicated study, and (2) the treatment of varsity sports on school websites is highly inconsistent. While some schools classify them as clubs, others categorize them separately under athletics. Furthermore, common varsity sports are frequently not reported on websites despite the likelihood that most schools offer them. For the sake of maintaining consistency in this study, we have excluded varsity sports and intend to explore them in future work. However, we do include intramural sports in this study, such as trapshooting and mountain biking, as they tend to be consistently reported as clubs on different school websites. Similarly, we excluded student councils and class-specific clubs such as Class 2025 because almost every school has these clubs, but the practice of whether to include them as clubs on school websites is inconsistent.

As several co-authors extracted club data from school websites, one challenge was ensuring inter-rater consistency. To address this, we established clear rules to exclude certain schools (e.g., private schools) and adopted a two-step process. First, a software tool randomly sampled thousands of schools from the NCES dataset. Co-authors then searched for each school’s website, identified the specific pages listing club data, and excluded schools that listed only a subset of clubs as examples rather than presenting a comprehensive list of all available clubs. This step filtered out approximately 95% of the sampled schools, leaving around 250. In the second step, a single evaluator reviewed the club webpages of the approximately 250 remaining schools and further excluded any that appeared to have incomplete club data, resulting in a final set of 229 schools. Having a single evaluator performing the second step ensured consistency and avoided inter-rater discrepancies in the final list. The process is scalable, as over 95% of the effort occurred in the first step, which was parallelized across multiple co-authors. A limitation, however, is that some valid schools may have been prematurely excluded in the first step by different data collectors, without a chance for the final evaluator to review them – potentially introducing bias and reducing the dataset size.

3.1 Effect of School Size

Our analysis starts with the school size factor. For an initial intuitive grasp of the effect of school size, we plot in Figure 2 the relationship between school size ( X enroll ) and the number of clubs ( Y club ) for each of the 229 sampled schools. Each dot represents one school, and the “moving average” curve shows the overall trend. To compute the moving-average curve, we first sort the schools according to their sizes. Then, for each group of 10 consecutive schools in the sorted list, we calculate their average school size ( x ) and average number of clubs ( y ). These mean values are then plotted as one data point ( x , y ) on the moving-average curve. By iterating this procedure for all groups of 10 consecutive schools in the sorted list, the complete moving-average curve can be plotted.

Figure 2

The number of clubs follows the trend of square root of school size. Without thorough modeling and exploration of different curve-fitting functions, one might oversimplify the relationship between school size and number of clubs, mistakenly assuming a linear correlation that could lead to suboptimal results.

The moving-average curve reveals a strong positive correlation between X enroll and Y club . The Pearson correlation coefficient is 0.63 , with p < 0.001 . To further understand their relationship, we fit the data to the following model: Y club = β 1 f ( X enroll ) + ε . Note that this model does not have a bias term because, intuitively, when X enroll approaches zero, Y club should also approach zero. This is confirmed by the trend in Figure 2.

The moving-average curve in Figure 2 shows that Y club grows more slowly than X enroll . This pattern suggests modeling the relationship with a slow-growing function. After exploring various forms of function f ( ⋅ ) , including f ( x ) = x , f ( x ) = x , f ( x ) = x 3 , f ( x ) = log ( x ) , and a logistic function, we find that f ( x ) = x yields the best fit, specifically, with the following coefficient:

This function is also plotted in Figure 2, showing good alignment with the moving-average curve. Note that we fit the function to the entire dataset instead of the moving-average curve, which merely aids visualization.

We derive several observations from Figure 2. First, the simplicity of the function is appealing, indicating a square root relationship between Y club and X enroll . Second, it reveals the nonlinear nature of their relationship, highlighting the necessity for a square root transformation of X enroll in the final multiple regression model. Finally, as the school size increases, the prediction error also increases, indicating the presence of heteroscedasticity. This issue requires attention in the final multiple regression model.

In Figure 2, one might note that a big fraction of the sampled schools are relatively small, i.e., with 500 or fewer students. This is because most schools are smaller schools, as shown by the curve of “by school count (all schools)” in Figure 1. To have sufficient samples for larger schools, we already intentionally sampled a bigger fraction of larger schools, as shown by the curve of “by school count (sampled schools)” in Figure 1.

Finally, from Figure 2, one might be inclined to assume that students at smaller schools are at a disadvantage, given their schools offer fewer clubs. However, Figure 3 presents a contrasting perspective, unveiling that the number of clubs per student is, in fact, higher in smaller schools. Because the number of clubs follows the trend of X enroll , the number of clubs per student follows the trend of X enroll X enroll = 1 X enroll . On the one hand, students at smaller schools enjoy better opportunities for leadership roles and active club engagement. On the other hand, smaller schools have a reduced variety of club choices.

Figure 3

The number of clubs per student decreases as the school size increases. On the one hand, smaller schools tend to have fewer clubs, limiting students’ choices. On the other hand, students at smaller schools enjoy better opportunities for leadership roles and active club engagement.

3.2 Effect of Pupil-to-Teacher Ratio

In addition to school size, the pupil-to-teacher ratio ( X teacher ) also affects the number of clubs. For an initial intuitive grasp, we graph in Figure 4 the 229 sampled schools based on their X teacher and X enroll values. Surprisingly, the figure appears to imply a positive correlation between X teacher and Y club , which is confirmed by a correlation coefficient of 0.20 and p < 0.01 . It seems counterintuitive for schools with a higher X teacher (meaning fewer teachers per student) to have more clubs, as schools rely on teachers to assume advisory roles for each club.

Figure 4

Schools with higher X teacher values tend to have more clubs, which is counterintuitive at first glance. This is because larger schools tend to be associated with both higher X teacher and higher Y club values, leading to the side effect of a positive correlation between X teacher and Y club .

Further investigation reveals that this counterintuitive result is merely a side effect of the correlation between X teacher and X enroll , as depicted in Figure 5. The seemingly paradoxical relationship between X teacher and Y club is primarily driven by the fact that larger schools, with higher X teacher , tend to have more clubs, rather than the direct effect of X teacher on Y club .

Figure 5

X teacher is correlated with X enroll . Larger schools tend to have higher X teacher values.

To assess the true effect of X teacher , we design a t -test between two groups of schools categorized by higher and lower X teacher values, respectively, while ensuring the two groups have comparable characteristics in school size. As portrayed in Figure 5, smaller schools generally exhibit lower X teacher values. Consequently, if we were to partition the groups solely based on raw X teacher values, the lower- X teacher group would largely consist of smaller schools, and the higher- X teacher group would mainly include larger schools. This would result in a comparison predominantly between schools of differing sizes, which is undesirable.

To prevent this, we follow the division lines in Figure 6 to create school groups for comparison. Specifically, we first sort the schools by their sizes. For every 30 adjacent schools in the sorted list with similar school sizes, we identify division lines that divide the 30 schools into three groups of equal sizes with higher, medium, and lower X teacher values, respectively. As observed in the figure, these division lines ascend as the school size increases, accommodating the tendency for larger schools to have higher X teacher values. Through this approach, schools assigned to the higher- and lower- X teacher groups possess comparable X enroll values, averaging 790 and 879, respectively. In contrast, their X teacher values differ significantly, averaging 19 and 11, respectively. These characteristics align well with the experiment’s objectives.

Figure 6

Partition the 229 sampled schools into three groups based on their X teacher values.

We compare the higher- and lower- X teacher groups in Figure 7. The figure shows that, after excluding the influence of school size, schools in the lower- X teacher group tend to have more clubs. A t -test between the two groups validates that the difference is statistically significant, with t -statistic = − 2.42  and p <  0.02. On average, schools in the lower- X teacher group have 9.3 more clubs than those in the higher- X teacher group. This contrasts with the seemingly positive correlation between X teacher and Y club .

Figure 7

Schools in the lower- X teacher group tend to have more clubs, meaning that among schools of similar sizes, those with a lower pupil-to-teacher ratio tend to have more clubs.

Furthermore, Figure 7 illustrates that the gap in Y club values between the two groups tends to widen as X enroll increases. This suggests that X enroll magnifies the effect of X teacher , which is why we will employ X enroll as a multiplicative factor for other independent variables such as X teacher in the final multiple regression model.

3.3 Effect of Household Income

Next, we analyze the effect of household income by using X lunch as an indicator of household income. X lunch is calculated as l ∕ X enroll , where l is the number of students eligible for free or reduced-price lunch. A lower X lunch indicates a higher income.

To assess whether we need to exclude the influence of school size before analyzing the effect of X lunch , we calculate the correlation between X lunch and school size. Unlike school size’s strong correlation with X teacher , it has little correlation with X lunch , reflected in a correlation coefficient of − 0.08  and p = 0.22 . This suggests that we can directly observe the effect of X lunch on Y club without needing to exclude the influence of school size.

Figure 8 illustrates the relationship between X lunch and Y club . When X lunch is less than 40%, an increase in X lunch tends to cause a decrease in Y club . However, when X lunch exceeds 40%, a further increase in X lunch no longer affects Y club . In general, the logistic function – visualized in Figure 9 and defined mathematically in equation (9) – is commonly used to represent this flattening-out effect. Specifically, the logistic function depicted in Figure 8 aligns well with the moving-average curve. Its parameters, K , X 0 , L , and C , are inferred using the Levenberg–Marquardt algorithm to fit the data. Our empirical evaluation further demonstrates that among various functions we explored, the logistic function indeed offers the best fit.

Figure 8

Only when X lunch < 40 % , higher-income schools tend to have more clubs than lower-income schools.

Figure 9

This logistic function g ( X lunch ) is used to transform X lunch in the final multiple regression model.

Because of the nonlinear relationship between X lunch and Y club , we will apply the following nonlinear transformation to X lunch in the final multiple regression model.

g ( X lunch ) is graphed in Figure 9. Within the domain of [0, 1] for X lunch , the range of g ( X lunch ) starts at 0.337 and levels off at 1 when X lunch approaches approximately 0.4.

Overall, because the curve in Figure 8 flattens out when X lunch exceeds 40%, it suggests that in terms of the effect on club count, there is a significant difference between the higher- and middle-income groups, but there is little difference between the middle- and lower-income groups. For context, the average value of X lunch for all schools is 50%.

To further quantify the statistical significance of the effect of household income, we follow the method shown in Figure 6 to partition the 229 sampled schools into three groups based on the values of X lunch : one-third higher-income schools (i.e., lower X lunch ), one-third middle-income schools, and one-third lower-income schools. We compare Y club of the higher- and lower-income groups, as depicted in Figure 10. The figure shows that the higher-income group has more clubs than the lower-income group, specifically, with an average difference of 16.7 clubs. For context, on average, a school has 26.1 clubs. A t -test between these groups confirms that their difference is statistically significant, with t -statistics = − 4.77  and p < 0.0001 .

Figure 10

Higher-income schools tend to have more clubs.

Similar to Figures 7 and 10 also shows that the disparity in Y club values between the two groups tends to expand as X enroll increases. This implies that X enroll amplifies the influence of X lunch , which is why we will include X enroll as a multiplicative factor for other independent variables in the final multiple regression model.

3.4 Effect of Racial Demographics

In this section, we analyze the effect of racial demographics. We begin by examining the largest racial group, the White students. For each school, we compute the fraction of White students as X white = w ∕ X enroll , where w is the number of White students in the school. X white is correlated with school size, with a coefficient of − 0.34  and p < 0.0001 . This indicates that larger schools tend to have smaller proportions of White students. Moreover, X white is even more strongly correlated with X lunch , with a coefficient of − 0.54  and p < 0.0001 . This indicates that lower-income schools tend to have smaller proportions of White students.

Since both X enroll and X lunch influence how X white affects Y club , to assess the true effect of X white in a controlled environment, we design an experiment to compare two groups of schools with comparable values in both X lunch and X enroll , but differ significantly in their X white values. The subsequent procedure outlines our approach. First, we sort the 229 sampled schools by their X lunch values and partition the sorted list into N lunch buckets. Each bucket comprises the same number of schools, and the schools within the same bucket have similar X lunch values as they are adjacent on the list sorted by X lunch . For schools in each bucket, we resort them by their X enroll values and further partition them into N enroll bins. Each bin comprises the same number of schools, and the schools within the same bin have similar X enroll values as they are adjacent on the list sorted by X enroll . In total, there are N lunch × N Enroll bins, where each bin has the same number of schools, and within each bin, schools exhibit similar values in both X lunch and X enroll . Finally, for the schools in each bin, we partition them into three subgroups based on their X white values, forming the higher-, medium-, and lower- X white subgroups. The aggregation of all schools in the higher- X white subgroups across all N lunch × N Enroll bins forms the overall higher- X white group. Similarly, the assembly of all schools in the N lunch × N Enroll lower- X white subgroups forms the overall lower- X white group. In our experiment, both N Lunch and N Enroll are set to 5.

The resulting higher- and lower- X white groups exhibit desired characteristics. Schools in these two groups have comparable X enroll values, averaging 969 and 994 students, respectively, along with comparable X lunch values, averaging 46 and 49%, respectively. However, their X white values differ significantly at 70 and 28%, respectively – meeting our intended criteria. Remarkably, despite substantial X white discrepancies, they have a comparable number of clubs ( Y club ), averaging 25.4 and 27.5, respectively. Interestingly, the higher- X white group shows a slightly lower number of clubs. The t -test between the Y club values of these two groups yields a p -value of 0.19, insufficient to establish a significant difference between them.

To visualize the pattern, we plot the schools in the higher- and lower- X white groups in Figure 11. The figure shows that despite the notable difference in their X white values, there is no distinct disparity in Y club between the two groups.

Figure 11

For two school groups with comparable X enroll and X lunch values but significant differences in their X white values, there is no evidence that X white significantly affects Y club .

Similar t -test results for other racial groups are summarized in Table 3. Overall, when comparing two groups of schools with higher and lower fractions of students from a specific racial group, similar to what is depicted in Figure 11, the marginal difference in Y club and the large p value in the t -test indicate a lack of statistical significance in the difference. An intuitive interpretation of these results is that schools with comparable student enrollments and household incomes do not display significant disparities in the number of clubs due to variations in racial demographics. The lack of direct affect on Y club is a key reason for us to consider excluding the factors representing racial demographics from the final multiple regression model, as described in the next section.

3.5 Multiple Regression Model

After understanding the characteristics of individual independent variables, in this section, we design a multiple regression model that incorporates all the factors. Specifically, we model the number of clubs as

where β i are coefficients, ε is the model residual, X race such as X white is the fraction of students of a specific race in a school, f ( X Enroll ) is the nonlinear transformation of X enroll derived from equation (1), and g ( X lunch ) is the nonlinear transformation of X lunch obtained from equation (2).

A distinctive aspect of this model is the role of f ( X enroll ) as a multiplicative factor for all other independent variables. This choice is closely tied to the significant role of school size. When f ( X enroll ) approaches 0, indicating a small student population, Y club should converge to zero, regardless of X teacher , X lunch , and X race . This behavior is effectively captured by utilizing f ( X enroll ) as a multiplicative term. Furthermore, the effect of other independent variables is amplified by school size. For instance, if higher income levels are assumed to lead to more clubs, this effect should result in a more substantial absolute increase in the number of clubs in larger schools. This amplification finds empirical support in the data, as depicted in Figures 7 and 10. Also note that this model does not have a bias term in order to achieve the effect that as f ( X enroll ) approaches 0, Y club converges to zero.

In addition to this model, we have explored numerous others, including those incorporating additional terms related to X lunch and X teacher but without f ( X enroll ) as their multiplicative factor, or adding a bias term, etc. Moreover, we have also explored introducing independent variables besides school demographics, such as those representing how schools present their club information on their websites, with or without teacher contact information, and with or without detailed club descriptions. Nevertheless, these models added complexity without meaningfully improving data fitting. Consequently, we have excluded their use.

3.6 Further Simplifying the Model

Although the model in equation (3) seems intuitive and comprehensive, we find that the terms representing racial demographics, ∑ race β race X race , do not contribute significantly to enhancing data fitting. Therefore, we eliminate those terms and adopt the simplified model below:

We use the linear regression implementation in the statsmodels package (Statsmodels, 2023) to compute the coefficients. Information concerning these coefficients is summarized in Table 4. Due to the presence of heteroscedasticity in the residual, as evident in Figures 7 and 10, we employ the HC3 covariance matrix estimator in statsmodels to calculate heteroskedasticity-robust standard errors. This estimator implements the algorithm proposed by MacKinnon and White (MacKinnon & White, 1985).

The adjusted R 2 for the simplified model is 0.60, while the p -value for the F-statistic is less than 1 0 − 60 , indicating that the model is a reasonable fit for the data. Note that when dealing with models lacking a bias term, the statsmodels package defaults to presenting uncentered and adjusted R 2 values, which is 0.84 for our model. When computing R 2 , it calculates the sum of squares total as ∑ Y club 2 , instead of the conventional ∑ ( Y club − Y ¯ club ) 2 due to the absence of the bias term to help center the residual. To avoid undue optimism, we have opted to focus on the lower, centered R 2 value of 0.60, instead of the uncentered R 2 value of 0.84.

The mean of the model residual, ε , is only − 0.44 , in comparison to the mean of Y club at 26.1. This suggests that the decision to exclude a bias term is justified. Furthermore, empirically, we note that the absence of the bias term enhances the model’s stability by preventing the regression from adjusting the bias term to compensate for errors elsewhere in the coefficients.

The signs of the coefficients in Table 4 align with the previous analysis results for individual variables. Specifically, X enroll shows a positive correlation with Y club , while X lunch and X teacher display negative correlations with Y club after excluding the influence of X enroll . Furthermore, the partial R 2 values for X enroll , X lunch , and X teacher are 0.50, 0.30, and 0.08, respectively. These values reflect the relative magnitude of their affect on Y club , aligning with our prior evaluation of these variables individually.

After establishing the simplified model in equation (5) as a solid baseline, we delve into the reasoning behind excluding the ∑ race β race X race terms. When solely introducing the X white variable, without incorporating other X race variables, to the simplified model, the adjusted R 2 remains 0.60. This indicates that X white does not contribute to enhancing data fit. On the other hand, if all X race variables are added to the model, they cause multicollinearity as the X race variables always sum up to one and have strong correlations. Hence, at least one X race variable should be dropped. Eliminating an X race variable representing a racial group with a small school population inadequately mitigates multicollinearity, as other high-value X race variables still have strong correlations. The most effective resolution is removing the X race variable with the largest value, namely, X white . With X white eliminated and multicollinearity resolved, incorporating additional X race variables, regardless of their combination, never raises the adjusted R 2 beyond 0.61, which is hardly an improvement over the simplified model.

These observations indicate that the inclusion of X race variables complicates the model without improving its data fit. Furthermore, the t -test results in Table 3 explicitly confirm that these variables do not have a significant effect on Y club for schools with comparable X enroll and X lunch values despite the significant variation in the X race values. Guided by these insights, we choose to exclude the X race terms from the simplified model.

3.7 School Initiatives in Improving Club Offerings

On the one hand, the multiple regression model’s reasonable accuracy, with an adjusted R 2 of 0.60, represents a noteworthy achievement, considering the intricate interplay of numerous big and small factors affecting operations of schools ranging from 50 to 4,300 students in this study. On the other hand, it underscores that schools still have ample opportunities to take initiatives to enhance club offerings, since 40% of the variance remains unattributed to hard-to-change factors like school size.

To illustrate the substantial effect of school initiatives, we compare top- and bottom-performing schools in terms of club offerings in Figure 12. For each school, we calculate the ratio between its actual club count and the club count prediction derived from equation (5). A school with a higher ratio outperforms a school with a lower ratio. In Figure 12, we plot the top and bottom 1/4 schools with the highest and lowest values in this ratio, respectively. This figure reveals a significant disparity between the top and bottom schools. Specifically, the average of the ratios for the top 1/4 schools is 163%, whereas that for the bottom 1/4 schools is only 51%. In other words, the top schools have 3.2 times more clubs than the bottom schools after school demographics are taken into account, emphasizing the pivotal role of individual school initiatives.

Figure 12

Out of the 229 sampled schools, this figure compares the relative club counts between the top 58 schools and the bottom 58 schools. After accounting for school demographics, the top and bottom schools still show a significant difference in club count, emphasizing the pivotal role of individual school initiatives.

3.8 Analysis Summary and Lessons Learned

Overall, the simplified model demonstrates a reasonable level of accuracy, indicated by its adjusted R 2 value of 0.60. Notable features of this model include: (1) addressing the nonlinearity of X enroll and X lunch through the transformations f ( X enroll ) and g ( X lunch ) as defined in equation (4); (2) utilizing X enroll exclusively as a multiplicative factor for X lunch and X teacher ; (3) refraining from incorporating a bias term to force Y club to approach zero as X enroll approaches zero; and (4) excluding the X race variables from the model despite the common perception that racial demographics might correlate with X enroll and X lunch , potentially influencing Y club .

A valuable lesson we have learned is that the design of this model greatly benefited from a thorough initial analysis of the characteristics of individual variables, rather than blindly including them in a linear regression model. This aided us in identifying and designing the aforementioned features of the model.

In contrast, if the traditional linear regression model shown below were used, it would produce misleading and unreasonable results despite its seemingly acceptable R 2 value of 0.51.

Specifically, the regression result for the aforementioned model misleadingly suggests that X teacher may not affect Y club due to β 3 ’s large p -value of 0.12. More importantly, the bias term, β 0 , has a value of 30.0, with p < 0.001 , indicating high confidence. However, this large bias is highly unreasonable since the mean of Y club is only 29.1, and when X enroll approaches zero, Y club should also approach zero, rather than taking on the large bias value of 30.0. These results highlight the importance of introducing the aforementioned notable features in our model.

Another valuable lesson we have learned is that when validating observations, we have frequently strived to design experiments that directly compare two population groups differing significantly in one independent variable while maintaining similar characteristics across other independent variables. Concrete examples include the comparisons shown in Figures 7, 10, and 11, which are further summarized in Table 3. This direct validation approach has enhanced our confidence in and improved our interpretation of the summaries generated by statistical software.