Source code for statsmodels.sandbox.regression.try_ols_anova

''' convenience functions for ANOVA type analysis with OLS

Note: statistical results of ANOVA are not checked, OLS is
checked but not whether the reported results are the ones used
in ANOVA

includes form2design for creating dummy variables

TODO:
 * ...
 *

'''

from __future__ import print_function
from statsmodels.compat.python import lmap
import numpy as np
#from scipy import stats
import statsmodels.api as sm

[docs]def data2dummy(x, returnall=False):
    '''convert array of categories to dummy variables
    by default drops dummy variable for last category
    uses ravel, 1d only'''
    x = x.ravel()
    groups = np.unique(x)
    if returnall:
        return (x[:, None] == groups).astype(int)
    else:
        return (x[:, None] == groups).astype(int)[:,:-1]

[docs]def data2proddummy(x):
    '''creates product dummy variables from 2 columns of 2d array

    drops last dummy variable, but not from each category
    singular with simple dummy variable but not with constant

    quickly written, no safeguards

    '''
    #brute force, assumes x is 2d
    #replace with encoding if possible
    groups = np.unique(lmap(tuple, x.tolist()))
    #includes singularity with additive factors
    return (x==groups[:,None,:]).all(-1).T.astype(int)[:,:-1]

[docs]def data2groupcont(x1,x2):
    '''create dummy continuous variable

    Parameters
    ----------
    x1 : 1d array
        label or group array
    x2 : 1d array (float)
        continuous variable

    Notes
    -----
    useful for group specific slope coefficients in regression
    '''
    if x2.ndim == 1:
        x2 = x2[:,None]
    dummy = data2dummy(x1, returnall=True)
    return dummy * x2

# Result strings
#the second leaves the constant in, not with NIST regression
#but something fishy with res.ess negative in examples ?
#not checked if these are all the right ones

anova_str0 = '''
ANOVA statistics (model sum of squares excludes constant)
Source    DF  Sum Squares   Mean Square    F Value    Pr > F
Model     %(df_model)i        %(ess)f       %(mse_model)f   %(fvalue)f %(f_pvalue)f
Error     %(df_resid)i     %(ssr)f       %(mse_resid)f
CTotal    %(nobs)i    %(uncentered_tss)f     %(mse_total)f

R squared  %(rsquared)f
'''

anova_str = '''
ANOVA statistics (model sum of squares includes constant)
Source    DF  Sum Squares   Mean Square    F Value    Pr > F
Model     %(df_model)i      %(ssmwithmean)f       %(mse_model)f   %(fvalue)f %(f_pvalue)f
Error     %(df_resid)i     %(ssr)f       %(mse_resid)f
CTotal    %(nobs)i    %(uncentered_tss)f     %(mse_total)f

R squared  %(rsquared)f
'''


[docs]def anovadict(res):
    '''update regression results dictionary with ANOVA specific statistics

    not checked for completeness
    '''
    ad = {}
    ad.update(res.__dict__)  #dict doesn't work with cached attributes
    anova_attr = ['df_model', 'df_resid', 'ess', 'ssr','uncentered_tss',
                 'mse_model', 'mse_resid', 'mse_total', 'fvalue', 'f_pvalue',
                  'rsquared']
    for key in anova_attr:
        ad[key] = getattr(res, key)
    ad['nobs'] = res.model.nobs
    ad['ssmwithmean'] = res.uncentered_tss - res.ssr
    return ad


[docs]def form2design(ss, data):
    '''convert string formula to data dictionary

    ss : string
     * I : add constant
     * varname : for simple varnames data is used as is
     * F:varname : create dummy variables for factor varname
     * P:varname1*varname2 : create product dummy variables for
       varnames
     * G:varname1*varname2 : create product between factor and
       continuous variable
    data : dict or structured array
       data set, access of variables by name as in dictionaries

    Returns
    -------
    vars : dictionary
        dictionary of variables with converted dummy variables
    names : list
        list of names, product (P:) and grouped continuous
        variables (G:) have name by joining individual names
        sorted according to input

    Examples
    --------
    >>> xx, n = form2design('I a F:b P:c*d G:c*f', testdata)
    >>> xx.keys()
    ['a', 'b', 'const', 'cf', 'cd']
    >>> n
    ['const', 'a', 'b', 'cd', 'cf']

    Notes
    -----

    with sorted dict, separate name list wouldn't be necessary
    '''
    vars = {}
    names = []
    for item in ss.split():
        if item == 'I':
            vars['const'] = np.ones(data.shape[0])
            names.append('const')
        elif not ':' in item:
            vars[item] = data[item]
            names.append(item)
        elif item[:2] == 'F:':
            v = item.split(':')[1]
            vars[v] = data2dummy(data[v])
            names.append(v)
        elif item[:2] == 'P:':
            v = item.split(':')[1].split('*')
            vars[''.join(v)] = data2proddummy(np.c_[data[v[0]],data[v[1]]])
            names.append(''.join(v))
        elif item[:2] == 'G:':
            v = item.split(':')[1].split('*')
            vars[''.join(v)] = data2groupcont(data[v[0]], data[v[1]])
            names.append(''.join(v))
        else:
            raise ValueError('unknown expression in formula')
    return vars, names

[docs]def dropname(ss, li):
    '''drop names from a list of strings,
    names to drop are in space delimeted list
    does not change original list
    '''
    newli = li[:]
    for item in ss.split():
        newli.remove(item)
    return newli

if __name__ == '__main__':

    # Test Example with created data
    # ------------------------------

    nobs = 1000
    testdataint = np.random.randint(3, size=(nobs,4)).view([('a',int),('b',int),('c',int),('d',int)])
    testdatacont = np.random.normal( size=(nobs,2)).view([('e',float), ('f',float)])
    import numpy.lib.recfunctions
    dt2 = numpy.lib.recfunctions.zip_descr((testdataint, testdatacont),flatten=True)
    # concatenate structured arrays
    testdata = np.empty((nobs,1), dt2)
    for name in testdataint.dtype.names:
        testdata[name] = testdataint[name]
    for name in testdatacont.dtype.names:
        testdata[name] = testdatacont[name]


    #print(form2design('a',testdata)

    if 0: # print(only when nobs is small, e.g. nobs=10
        xx, n = form2design('F:a',testdata)
        print(xx)
        print(form2design('P:a*b',testdata))
        print(data2proddummy((np.c_[testdata['a'],testdata['b']])))

        xx, names = form2design('a F:b P:c*d',testdata)

    #xx, names = form2design('I a F:b F:c F:d P:c*d',testdata)
    xx, names = form2design('I a F:b P:c*d', testdata)
    xx, names = form2design('I a F:b P:c*d G:a*e f', testdata)


    X = np.column_stack([xx[nn] for nn in names])
    # simple test version: all coefficients equal to one
    y = X.sum(1) + 0.01*np.random.normal(size=(nobs))
    rest1 = sm.OLS(y,X).fit() #results
    print(rest1.params)
    print(anova_str % anovadict(rest1))


    X = np.column_stack([xx[nn] for nn in dropname('ae f', names)])
    # simple test version: all coefficients equal to one
    y = X.sum(1) + 0.01*np.random.normal(size=(nobs))
    rest1 = sm.OLS(y,X).fit()
    print(rest1.params)
    print(anova_str % anovadict(rest1))


    # Example: from Bruce
    # -------------------

    #get data and clean it
    #^^^^^^^^^^^^^^^^^^^^^

    # requires file 'dftest3.data' posted by Bruce

    # read data set and drop rows with missing data
    dt_b = np.dtype([('breed', int), ('sex', int), ('litter', int),
                   ('pen', int), ('pig', int), ('age', float),
                   ('bage', float), ('y', float)])
    dta = np.genfromtxt('dftest3.data', dt_b,missing='.', usemask=True)
    print('missing', [dta.mask[k].sum() for k in dta.dtype.names])
    m = dta.mask.view(bool)
    droprows = m.reshape(-1,len(dta.dtype.names)).any(1)
    # get complete data as plain structured array
    # maybe doesn't work with masked arrays
    dta_use_b1 = dta[~droprows,:].data
    print(dta_use_b1.shape)
    print(dta_use_b1.dtype)

    #Example b1: variables from Bruce's glm
    #^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

    # prepare data and dummy variables
    xx_b1, names_b1 = form2design('I F:sex age', dta_use_b1)
    # create design matrix
    X_b1 = np.column_stack([xx_b1[nn] for nn in dropname('', names_b1)])
    y_b1 = dta_use_b1['y']
    # estimate using OLS
    rest_b1 = sm.OLS(y_b1, X_b1).fit()
    # print(results)
    print(rest_b1.params)
    print(anova_str % anovadict(rest_b1))
    #compare with original version only in original version
    #print(anova_str % anovadict(res_b0))

    # Example: use all variables except pig identifier
    #^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

    allexog = ' '.join(dta.dtype.names[:-1])
    #'breed sex litter pen pig age bage'

    xx_b1a, names_b1a = form2design('I F:breed F:sex F:litter F:pen age bage', dta_use_b1)
    X_b1a = np.column_stack([xx_b1a[nn] for nn in dropname('', names_b1a)])
    y_b1a = dta_use_b1['y']
    rest_b1a = sm.OLS(y_b1a, X_b1a).fit()
    print(rest_b1a.params)
    print(anova_str % anovadict(rest_b1a))

    for dropn in names_b1a:
        print(('\nResults dropping', dropn))
        X_b1a_ = np.column_stack([xx_b1a[nn] for nn in dropname(dropn, names_b1a)])
        y_b1a_ = dta_use_b1['y']
        rest_b1a_ = sm.OLS(y_b1a_, X_b1a_).fit()
        #print(rest_b1a_.params)
        print(anova_str % anovadict(rest_b1a_))