Module stikpetP.other.poho_residual_gof_bin

Functions

def ph_residual_gof_bin(data, test='std-residual', expCount=None, mtc='bonferroni', **kwargs)

Post-Hoc Residuals Goodness-of-Fit Test

This function will perform a residuals post-hoc test for each of the categories in a nominal field. This could either be a z-test using the standardized residuals, the adjusted residuals, or any of the one-sample binary tests.

The unadjusted p-values and Bonferroni adjusted p-values are both determined.

This function is shown in this YouTube video and the test is also described at PeterStatistics.com

Parameters

data : list or pandas series

 

test : {"adj-residual", "std-residual", "binomial", "wald", "score"}, optional

test to use

expCount : pandas dataframe, optional

categories and expected counts

mtc : string, optional

any of the methods available in p_adjust() to correct for multiple tests

**kwargs : optional

additional arguments for the specific test that are passed along.

Returns

pandas.DataFrame

A dataframe with the following columns:

• category, the label of the category
• obs. count, the observed count
• exp. count, the expected count
• z-statistic, the standardized residuals
• p-value, the unadjusted significance
• adj. p-value, the adjusted significance
• test, description of the test used.

Notes

The formula used is for the adjusted residual test: $$z = \frac{F_i - E_i}{\sqrt{E_i\times\left(1 - \frac{E_i}{n}\right)}}$$ $$sig = 2\times\left(1 - \Phi\left(\left|z\right|\right)\right)$$

The formula used is for the standardized residual test: $$z = \frac{F_i - E_i}{\sqrt{E_i}}$$ $$sig = 2\times\left(1 - \Phi\left(\left|z\right|\right)\right)$$

With: * $F_i$, the observed count for category $i$ * $E_i$, the expected count for category $i$ * $\Phi\left(\dots\right)$, the cumulative distribution function of the standard normal distribution

If no expected counts are provide it is assumed they are all equal for each category, i.e. $E_i = \frac{n}{k}$

The Bonferroni adjustment is calculated using: $$p_{adj} = \min \left(p \times k, 1\right)$$

The other tests use the formula from the one-sample test variant, using the expected count/n as the expected proportion.

The adjusted residuals will gave the same result as using a one-sample score test. Some sources will also call these adjusted residuals as standardized residuals (Agresti, 2007, p. 38), and the standardized residuals used in this function as Pearson residuals (R, n.d.). Haberman (1973, p. 205) and Sharpe (2015, p. 3) are sources for the terminology used in this function.

Before, After and Alternatives

Before this an omnibus test might be helpful: * ts_pearson_gof for Pearson Chi-Square Goodness-of-Fit Test * ts_freeman_tukey_gof for Freeman-Tukey Test of Goodness-of-Fit * ts_freeman_tukey_read for Freeman-Tukey-Read Test of Goodness-of-Fit * ts_g_gof for G (Likelihood Ratio) Goodness-of-Fit Test * ts_mod_log_likelihood_gof for Mod-Log Likelihood Test of Goodness-of-Fit * ts_multinomial_gof for Multinomial Goodness-of-Fit Test * ts_neyman_gof for Neyman Test of Goodness-of-Fit * ts_powerdivergence_gof for Power Divergence GoF Test

After this you might want to add an effect size measure: * es_post_hoc_gof for various effect sizes

Alternative post-hoc tests: * ph_pairwise_bin for Pairwise Binary Test * ph_pairwise_gof for Pairwise Goodness-of-Fit Tests * ph_residual_gof_gof for Residuals Using Goodness-of-Fit Tests

The binary test that is performed on each category: * ts_binomial_os for One-Sample Binomial Test * ts_score_os for One-Sample Score Test * ts_wald_os for One-Sample Wald Test

More info on the adjustment for multiple testing: * p_adjust

References

Agresti, A. (2007). An introduction to categorical data analysis (2nd ed.). Wiley-Interscience.

Haberman, S. J. (1973). The analysis of residuals in cross-classified tables. Biometrics, 29(1), 205–220. doi:10.2307/2529686

Mangiafico, S. S. (2016). Summary and analysis of extension program evaluation in R (1.20.01). Rutger Cooperative Extension.

R. (n.d.). Chisq.test [Computer software]. https://stat.ethz.ch/R-manual/R-devel/library/stats/html/chisq.test.html

Sharpe, D. (2015). Your chi-square test is statistically significant: Now what? Practical Assessment, Research & Evaluation, 20(8), 1–10.

Author

Made by P. Stikker

Companion website: https://PeterStatistics.com
YouTube channel: https://www.youtube.com/stikpet
Donations: https://www.patreon.com/bePatron?u=19398076

Examples

Examples: get data

>>> import pandas as pd
>>> pd.set_option('display.width',1000)
>>> pd.set_option('display.max_columns', 1000)    
>>> gss_df = pd.read_csv('https://peterstatistics.com/Packages/ExampleData/GSS2012a.csv', sep=',', low_memory=False, storage_options={'User-Agent': 'Mozilla/5.0'})
>>> ex1 = gss_df['mar1'];

Example 1 using default settings:

>>> ph_residual_gof_bin(ex1)
        category  obs. count  exp. count  z-statistic   p-value  adj. p-value                           test
0        MARRIED       972.0       388.2    29.630319  0.000000       0.00000  standardized residuals z-test
1  NEVER MARRIED       395.0       388.2     0.345129  0.729998       1.00000  standardized residuals z-test
2       DIVORCED       314.0       388.2    -3.765964  0.000166       0.00083  standardized residuals z-test
3        WIDOWED       181.0       388.2   -10.516276  0.000000       0.00000  standardized residuals z-test
4      SEPARATED        79.0       388.2   -15.693208  0.000000       0.00000  standardized residuals z-test

Example 2 using a binomial test and Holm correction:

>>> ph_residual_gof_bin(ex1, test="binomial", mtc='holm')
        category  obs. count  exp. count z-statistic        p-value   adj. p-value                                                                         test
0        MARRIED       972.0       388.2        n.a.  2.375122e-191  1.187561e-190        one-sample binomial, with equal-distance method (with p0 for MARRIED)
1  NEVER MARRIED       395.0       388.2        n.a.   3.585178e-01   3.585178e-01  one-sample binomial, with equal-distance method (with p0 for NEVER MARRIED)
2       DIVORCED       314.0       388.2        n.a.   3.215765e-05   6.431530e-05       one-sample binomial, with equal-distance method (with p0 for DIVORCED)
3        WIDOWED       181.0       388.2        n.a.   9.164123e-38   2.749237e-37        one-sample binomial, with equal-distance method (with p0 for WIDOWED)
4      SEPARATED        79.0       388.2        n.a.   3.303218e-94   1.321287e-93      one-sample binomial, with equal-distance method (with p0 for SEPARATED)

Example from Mangiafico (2016, p. 523-524):

>>> observed = ["White"]*20 + ["Black"]*9 + ["AI"]*9 + ["Asian"] + ["PI"] + ["Two+"]
>>> expected = pd.DataFrame({'category' : ["White", "Black", "AI", "Asian", "PI", "Two+"], 'count' : [0.775*41, 0.132*41, 0.012*41, 0.054*41, 0.002*41, 0.025*41]})
>>> ph_residual_gof_bin(observed, expCount=expected, test="adj-residual")
  category  obs. count  exp. count  z-statistic   p-value  adj. p-value                       test
0    White        20.0      31.775    -4.403793  0.000011      0.000064  adjusted residuals z-test
1    Black         9.0       5.412     1.655441  0.097835      0.587011  adjusted residuals z-test
2       AI         9.0       0.492    12.202996  0.000000      0.000000  adjusted residuals z-test
3    Asian         1.0       2.214    -0.838850  0.401553      1.000000  adjusted residuals z-test
4       PI         1.0       0.082     3.209006  0.001332      0.007992  adjusted residuals z-test
5     Two+         1.0       1.025    -0.025008  0.980049      1.000000  adjusted residuals z-test

Expand source code

def ph_residual_gof_bin(data, test="std-residual", expCount=None, mtc='bonferroni', **kwargs):
    '''
    Post-Hoc Residuals Goodness-of-Fit Test
    ---------------------------------------
    
    This function will perform a residuals post-hoc test for each of the categories in a nominal field. This could either be a z-test using the standardized residuals, the adjusted residuals, or any of the one-sample binary tests.

    The unadjusted p-values and Bonferroni adjusted p-values are both determined.

    This function is shown in this [YouTube video](https://youtu.be/YjhGhsphXRQ) and the test is also described at [PeterStatistics.com](https://peterstatistics.com/Terms/Tests/PostHocAfterGoF.html)
    
    Parameters
    ----------
    data : list or pandas series
    test : {"adj-residual", "std-residual", "binomial", "wald", "score"}, optional 
        test to use
    expCount : pandas dataframe, optional 
        categories and expected counts
    mtc : string, optional
        any of the methods available in p_adjust() to correct for multiple tests
    **kwargs : optional
        additional arguments for the specific test that are passed along.
    
    Returns
    -------
    pandas.DataFrame
        A dataframe with the following columns:
    
        - *category*, the label of the category
        - *obs. count*, the observed count
        - *exp. count*, the expected count
        - *z-statistic*, the standardized residuals
        - *p-value*, the unadjusted significance
        - *adj. p-value*, the adjusted significance
        - *test*, description of the test used.

    Notes
    -----
    The formula used is for the adjusted residual test:
    $$z = \\frac{F_i - E_i}{\\sqrt{E_i\\times\\left(1 - \\frac{E_i}{n}\\right)}}$$
    $$sig = 2\\times\\left(1 - \\Phi\\left(\\left|z\\right|\\right)\\right)$$
    
    The formula used is for the standardized residual test:
    $$z = \\frac{F_i - E_i}{\\sqrt{E_i}}$$
    $$sig = 2\\times\\left(1 - \\Phi\\left(\\left|z\\right|\\right)\\right)$$
    
    With:
    * $F_i$, the observed count for category $i$
    * $E_i$, the expected count for category $i$
    * $\\Phi\\left(\\dots\\right)$, the cumulative distribution function of the standard normal distribution

    If no expected counts are provide it is assumed they are all equal for each category, i.e. $E_i = \\frac{n}{k}$
    
    The Bonferroni adjustment is calculated using:
    $$p_{adj} = \\min \\left(p \\times k, 1\\right)$$

    The other tests use the formula from the one-sample test variant, using the expected count/n as the expected proportion.

    The adjusted residuals will gave the same result as using a one-sample score test. Some sources will also call these adjusted residuals as standardized residuals (Agresti, 2007, p. 38), and the standardized residuals used in this function as Pearson residuals (R, n.d.). Haberman (1973, p. 205) and Sharpe (2015, p. 3) are sources for the terminology used in this function. 

    Before, After and Alternatives
    ------------------------------
    Before this an omnibus test might be helpful:
    * [ts_pearson_gof](../tests/test_pearson_gof.html#ts_pearson_gof) for Pearson Chi-Square Goodness-of-Fit Test
    * [ts_freeman_tukey_gof](../tests/test_freeman_tukey_gof.html#ts_freeman_tukey_gof) for Freeman-Tukey Test of Goodness-of-Fit
    * [ts_freeman_tukey_read](../tests/test_freeman_tukey_read.html#ts_freeman_tukey_read) for Freeman-Tukey-Read Test of Goodness-of-Fit
    * [ts_g_gof](../tests/test_g_gof.html#ts_g_gof) for G (Likelihood Ratio) Goodness-of-Fit Test
    * [ts_mod_log_likelihood_gof](../tests/test_mod_log_likelihood_gof.html#ts_mod_log_likelihood_gof) for Mod-Log Likelihood Test of Goodness-of-Fit
    * [ts_multinomial_gof](../tests/test_multinomial_gof.html#ts_multinomial_gof) for Multinomial Goodness-of-Fit Test
    * [ts_neyman_gof](../tests/test_neyman_gof.html#ts_neyman_gof) for Neyman Test of Goodness-of-Fit
    * [ts_powerdivergence_gof](../tests/test_powerdivergence_gof.html#ts_powerdivergence_gof) for Power Divergence GoF Test
    
    After this you might want to add an effect size measure:
    * [es_post_hoc_gof](../effect_sizes/eff_size_post_hoc_gof.html#es_post_hoc_gof) for various effect sizes
    
    Alternative post-hoc tests:
    * [ph_pairwise_bin](../other/poho_pairwise_bin.html#ph_pairwise_bin) for Pairwise Binary Test
    * [ph_pairwise_gof](../other/poho_pairwise_gof.html#ph_pairwise_gof) for Pairwise Goodness-of-Fit Tests
    * [ph_residual_gof_gof](../other/poho_residual_gof_gof.html#ph_residual_gof_gof) for Residuals Using Goodness-of-Fit Tests

    The binary test that is performed on each category:
    * [ts_binomial_os](../tests/test_binomial_os.html#ts_binomial_os) for One-Sample Binomial Test
    * [ts_score_os](../tests/test_score_os.html#ts_score_os) for One-Sample Score Test
    * [ts_wald_os](../tests/test_wald_os.html#ts_wald_os) for One-Sample Wald Test

    More info on the adjustment for multiple testing:
    * [p_adjust](../other/p_adjustments.html#p_adjust)
    
    References
    ----------
    Agresti, A. (2007). *An introduction to categorical data analysis* (2nd ed.). Wiley-Interscience.
    
    Haberman, S. J. (1973). The analysis of residuals in cross-classified tables. *Biometrics, 29*(1), 205–220. doi:10.2307/2529686
    
    Mangiafico, S. S. (2016). Summary and analysis of extension program evaluation in R (1.20.01). Rutger Cooperative Extension.
    
    R. (n.d.). Chisq.test [Computer software]. https://stat.ethz.ch/R-manual/R-devel/library/stats/html/chisq.test.html
    
    Sharpe, D. (2015). Your chi-square test is statistically significant: Now what? Practical Assessment, *Research & Evaluation, 20*(8), 1–10.
    
    Author
    ------
    Made by P. Stikker
    
    Companion website: https://PeterStatistics.com  
    YouTube channel: https://www.youtube.com/stikpet  
    Donations: https://www.patreon.com/bePatron?u=19398076

    Examples
    --------
    Examples: get data
    >>> import pandas as pd
    >>> pd.set_option('display.width',1000)
    >>> pd.set_option('display.max_columns', 1000)    
    >>> gss_df = pd.read_csv('https://peterstatistics.com/Packages/ExampleData/GSS2012a.csv', sep=',', low_memory=False, storage_options={'User-Agent': 'Mozilla/5.0'})
    >>> ex1 = gss_df['mar1'];

    Example 1 using default settings:
    >>> ph_residual_gof_bin(ex1)
            category  obs. count  exp. count  z-statistic   p-value  adj. p-value                           test
    0        MARRIED       972.0       388.2    29.630319  0.000000       0.00000  standardized residuals z-test
    1  NEVER MARRIED       395.0       388.2     0.345129  0.729998       1.00000  standardized residuals z-test
    2       DIVORCED       314.0       388.2    -3.765964  0.000166       0.00083  standardized residuals z-test
    3        WIDOWED       181.0       388.2   -10.516276  0.000000       0.00000  standardized residuals z-test
    4      SEPARATED        79.0       388.2   -15.693208  0.000000       0.00000  standardized residuals z-test

    Example 2 using a binomial test and Holm correction:
    >>> ph_residual_gof_bin(ex1, test="binomial", mtc='holm')
            category  obs. count  exp. count z-statistic        p-value   adj. p-value                                                                         test
    0        MARRIED       972.0       388.2        n.a.  2.375122e-191  1.187561e-190        one-sample binomial, with equal-distance method (with p0 for MARRIED)
    1  NEVER MARRIED       395.0       388.2        n.a.   3.585178e-01   3.585178e-01  one-sample binomial, with equal-distance method (with p0 for NEVER MARRIED)
    2       DIVORCED       314.0       388.2        n.a.   3.215765e-05   6.431530e-05       one-sample binomial, with equal-distance method (with p0 for DIVORCED)
    3        WIDOWED       181.0       388.2        n.a.   9.164123e-38   2.749237e-37        one-sample binomial, with equal-distance method (with p0 for WIDOWED)
    4      SEPARATED        79.0       388.2        n.a.   3.303218e-94   1.321287e-93      one-sample binomial, with equal-distance method (with p0 for SEPARATED)

    Example from Mangiafico (2016, p. 523-524):
    >>> observed = ["White"]*20 + ["Black"]*9 + ["AI"]*9 + ["Asian"] + ["PI"] + ["Two+"]
    >>> expected = pd.DataFrame({'category' : ["White", "Black", "AI", "Asian", "PI", "Two+"], 'count' : [0.775*41, 0.132*41, 0.012*41, 0.054*41, 0.002*41, 0.025*41]})
    >>> ph_residual_gof_bin(observed, expCount=expected, test="adj-residual")
      category  obs. count  exp. count  z-statistic   p-value  adj. p-value                       test
    0    White        20.0      31.775    -4.403793  0.000011      0.000064  adjusted residuals z-test
    1    Black         9.0       5.412     1.655441  0.097835      0.587011  adjusted residuals z-test
    2       AI         9.0       0.492    12.202996  0.000000      0.000000  adjusted residuals z-test
    3    Asian         1.0       2.214    -0.838850  0.401553      1.000000  adjusted residuals z-test
    4       PI         1.0       0.082     3.209006  0.001332      0.007992  adjusted residuals z-test
    5     Two+         1.0       1.025    -0.025008  0.980049      1.000000  adjusted residuals z-test
    
    '''
    if type(data) is list:
        data = pd.Series(data)
            
    freq = data.value_counts()
        
    if expCount is None:
        #assume all to be equal
        n = sum(freq)
        k = len(freq)
        categories = list(freq.index)
        expC = [n/k] * k
        
    else:
        #check if categories match
        nE = 0
        n = 0
        for i in range(0, len(expCount)):
            nE = nE + expCount.iloc[i,1]
            n = n + freq[expCount.iloc[i,0]]
        
        expC = []
        for i in range(0,len(expCount)):
            expC.append(expCount.iloc[i, 1]/nE*n)
            
        k = len(expC)
        categories = list(expCount.iloc[:,0])
    
    results = pd.DataFrame()
    resRow=0
    for i in range(0, k):
        cat = categories[i]
        n1 = freq[categories[i]]
        e1 = expC[i]
        if test=="std-residual":
            statistic = (n1 - e1)/(e1**0.5)
            sig = 2 * (1 - NormalDist().cdf(abs(statistic))) 
            test_description = "standardized residuals z-test"
        elif test=="adj-residual":
            statistic = (n1 - e1)/((e1*(1 - e1/n))**0.5)
            sig = 2 * (1 - NormalDist().cdf(abs(statistic))) 
            test_description = "adjusted residuals z-test"
        else:
            tempA = [categories[i]]*n1
            tempB = ['all other']*(n - n1)
            temp_data = tempA + tempB
            if test=="binomial":
                test_result = ts_binomial_os(temp_data, codes=[categories[i], 'all other'], p0=e1/n, **kwargs)
                statistic = "n.a."
                sig = test_result.iloc[0, 0]
                test_description = test_result.iloc[0, 1]
            else:
                if test=="wald":
                    test_result = ts_wald_os(temp_data, codes=[categories[i], 'all other'], p0=e1/n, **kwargs)
                elif test=="score":
                    test_result = ts_score_os(temp_data, codes=[categories[i], 'all other'], p0=e1/n, **kwargs)
    
                statistic = test_result.iloc[0, 1]
                sig = test_result.iloc[0, 2]
                test_description = test_result.iloc[0, 3]
        
                
        results.at[i, 0] = cat
        results.at[i, 1] = n1
        results.at[i, 2] = e1
        results.at[i, 3] = statistic
        results.at[i, 4] = sig
        results.at[i, 5] = sig
        results.at[i, 6] = test_description

    results.iloc[:,5] = p_adjust(results.iloc[:,4], method=mtc)
    results.columns = ["category", "obs. count", "exp. count", "z-statistic", "p-value", "adj. p-value", "test"]
    
    return results