live, firsts, others = first.MakeFrames()
formula = 'totalwgt_lb ~ agepreg'
model = smf.ols(formula, data=live)
results = model.fit()
regression.SummarizeResults(results)

##Results are also available as parameters:
inter = results.params['Intercept']
slope = results.params['agepreg']
slope_pvalue = results.pvalues['agepreg']
results.rsquared

##this gives the p-value associated with the model as a whole
results.f_pvalue

residuals = results.resid

##this returns a sequence of values corresponding
##to agepreg.
fitted_values = results.fittedvalues

# results.summary() provides a lot of info
#the following is easier:
regression.SummarizeResults(results)

"""
Explanation: regression - fitting a model to data.  The goal is to describe the relationship between the dependent variables and the explanatory variables.
multiple regression multiple explanatory variables.
linear regression wen the relationship is linear.
$$
y = \beta_0 + \beta_1x_1 + \beta_2x_2 + \epsilon
$$
where $\beta_0$ is the intercept, $\beta_1$ is the parameter association with $x_1$, $\beta_2$ is the parameter associated with $x_2$ and $\epsilon$ is the residual due to random variation or unknown factors.
ordinary least squares.  Given a sequence of values for y and sequences for $x_1$ and $x_2$, we can find the beta parameters that minimize $epsilon^2$
End of explanation
"""


diff_weight = firsts.totalwgt_lb.mean() - others.totalwgt_lb.mean()
diff_weight

diff_age = firsts.agepreg.mean() - others.agepreg.mean()
diff_age

results = smf.ols('totalwgt_lb ~ agepreg', data=live).fit()
slope = results.params['agepreg']
slope

##expected difference in birthweight for
##first babies and others, due to mother's age:
slope * diff_age

##this only accounts for half of the observed difference.
live['isfirst'] = live.birthord == 1
formula = 'totalwgt_lb ~ isfirst'
results = smf.ols(formula, data=live).fit()

regression.SummarizeResults(results)
print 

formula = 'totalwgt_lb ~ agepreg'
results = smf.ols(formula, data=live).fit()

regression.SummarizeResults(results)

"""
Explanation: spurious - result for which there is no obvious mechanism that would explain it. e.g. why would first babies be lighter than others?  Perhaps because mothers of first babies are younger...
End of explanation
"""


formula = 'totalwgt_lb ~ isfirst + agepreg'
results = smf.ols(formula, data=live).fit()
regression.SummarizeResults(results)

"""
Explanation: isfirst is a categorical variable, i.e. it falls into categories like True or False.  The estimated parameter is the effect on birth weight when isfirst is true.  So the result, -0.125 lbs is the difference in birth weight between first babies and others.
End of explanation
"""


live['agepreg2'] = live.agepreg**2
formula = 'totalwgt_lb ~ isfirst + agepreg + agepreg2'
results = smf.ols(formula, data=live).fit()
regression.SummarizeResults(results)

"""
Explanation: In the combined model, the parameter isfirst is smaller by about half, which means that part of the apparent effict is accounted for by agepreg.  Also, notice the p-value is now on the border of statistical significance.
could the contribution of agepreg be nonlinear? (Note, this is still linear regression, because the explanatory variable is now just the result of a nonlinear function)
End of explanation
"""


live = live[live.prglngth>30]
resp = chap01soln.ReadFemResp()
resp.index = resp.caseid
join = live.join(resp, on='caseid', rsuffix='_r')

t = []
for name in join.columns:

    try:
        if join[name].var() < 1e-7:
            continue
            
        formula = 'totalwgt_lb ~ agepreg + ' + str(name)
        model = smf.ols(formula, data=join)
        if model.nobs < len(join)/2: #num of observations
            continue
        
        results = model.fit()
    except(ValueError, TypeError):
        continue
        
    t.append((results.rsquared, name))
t

t.sort(reverse=True)
for mse, name in t[:30]:
    print name, mse

"""
Explanation: Negative coefficient in agepreg2 corroborates with downward curving of figure 10.2
in this example, the mother's age acts as a control variable for the diference in age between first-time mothers and others, making it possible to isolate the effect of isfirst
End of explanation
"""


formula = ('totalwgt_lb ~ agepreg + C(race) + babysex==1 + '
           'nbrnaliv > 1 +paydu==1 + totincr')
results = smf.ols(formula, data=join).fit()

regression.SummarizeResults(results)

"""
Explanation: race is a proxy variable; that is apparent correlations are often caused, at least in part, by other factors.
data mining
    Sometimes you start with a theory and use data to test it.  Other times you start with data and go looking for possible theories.
End of explanation
"""


y = np.array([0,1,0,1])
x1 = np.array([0,0,0,1])
x2 = np.array([0,1,1,1])

beta = [-1.5,2.8,1.1]

log_o = beta[0] + beta[1]*x1 + beta[2]*x2
print "log_o", log_o
o = np.exp(log_o)
print "o", o
p = o / (o + 1)
print "p", p

likes = y * p + (1-y) * (1-p)
print "likes", likes
like = np.prod(likes)
print "overall likelihood", like

"""
Explanation: C(race) tells Patsy to treat race as a categorical variable
Poisson regression if the dependent variable is an integer count.
logisitc regression when the dependet variable is boolean.  Expresses predictions in odds rather than probabilities--i.e. the ratio of the probability that it will occur to the probability that it will not.
$$
\text{odds} = p / (1 - p)
$$
$$
\text{probability} =  o / (o + 1)
$$
Logistic regression is based on the following model:
$$
\log o = \beta_0 + \beta_1x_1 + \beta_2x_2 + \epsilon
$$
logistic regression does not have a closed form solution, so it is solved by guessing an initial solution and improving it iteratively--looking for MLE, or the set of parameters that maximize the likelyhood of the data.
End of explanation
"""


live, first, others = first.MakeFrames()
df = live[live.prglngth>30]

df['boy'] = (df.babysex==1).astype(int)
model = smf.logit('boy ~ agepreg', data=df)
results = model.fit()
regression.SummarizeResults(results)

formula = 'boy ~ agepreg + hpagelb + birthord + C(race)'
model = smf.logit(formula, data=df)
results = model.fit()
regression.SummarizeResults(results)

##these dataframes contain the variables
endog = pandas.DataFrame(model.endog, columns=[model.endog_names])
exog = pandas.DataFrame(model.exog, columns=model.exog_names)
exog.head(5)

"""
Explanation: Probability is p when y == 1, 1-p when y == 0.  The goal is to maximize the likelihood of the data.  
endogenous variable - dependent variable
exogenous - explanatory
End of explanation
"""


actual = endog['boy']
baseline = actual.mean()
#this is the baseline prediction if you were to grab
#a random record from the sample
baseline

predict = (results.predict() >= 0.5)

true_pos = predict * actual
true_neg = (1 - predict) * (1 - actual)
##accuracy is a fraction of correct guesses
acc =  (sum(true_pos) + sum(true_neg)) /  len(actual)
print acc

columns = ['agepreg', 'hpagelb','birthord', 'race']
new = pandas.DataFrame([[35, 39, 3, 2]], columns=columns)
y = results.predict(new)
y

"""
Explanation: Measuring Accuracy
End of explanation
"""


df = join[join.prglngth > 30]
def datamine(df, dep_variable):
    """
    expl_variables:  a list of explanatory variables to
        be considered in the mining.
    """
    t = []
    for name in df.columns:
        try:
            if df[name].var() < 1e-7:
                continue
            
            formula = str(dep_variable) + " ~ " + str(name)
            model = smf.ols(formula, data=df)
            if model.nobs < len(df) / 2:
                continue
                
            results = model.fit()
        except (ValueError, TypeError):
            continue
        
        t.append((name,results))
    return t
   
t = datamine(df, 'prglngth')

def SortAndPrint(t):
    """
    t is a list of tuples with regression results 
       in the form (name of parameter, results object)
    """
    t.sort(key=lambda x: x[1].rsquared, reverse=True)
    for name, res in t:
        try:
            if res.pvalues[name] < 0.1:
                print name, res.params[name], res.rsquared
        except KeyError, e:
            print "EXCEPTION, ",e
SortAndPrint(t)

model = smf.ols('prglngth ~ birthord + C(race) + C(paydu) + + C(diabetes)', data=df)
results = model.fit()
regression.SummarizeResults(results)

"""
Explanation: Exercises
exercise 11.1
End of explanation
"""


def GoMiningBabySex(df):
    df.dropna(subset=['babysex'])
    df['boy'] = (df.babysex==1).astype(int)
    
    t = []
    
    for name in df.columns:
        try:
            if df[name].var() < 1e-7:
                continue
            
            formula = "boy ~ " + str(name)
            model = smf.logit(formula, data=df)
            nobs = len(model.endog)
            if nobs < len(df) / 2:
                continue
                
            results = model.fit()
        except:
            continue
        
        t.append((results.prsquared, name))
    return t

bs = GoMiningBabySex(df)


regression.MiningReport(bs)

model = smf.logit('boy ~ fmarout5==5 + totalwgt_lb + nplaced + infever==1', data=df)
results = model.fit()
regression.SummarizeResults(results)

"""
Explanation: Exercise 11.2
End of explanation
"""


def makeJoin():
    live, firsts, others = first.MakeFrames()
    live = live[live.prglngth>30]
    resp = chap01soln.ReadFemResp()
    resp.index = resp.caseid
    join = live.join(resp, on='caseid', rsuffix='_r')
    return join

join = makeJoin()

join = join.dropna(subset=['agepreg'])
join.numbabes.replace([97], np.nan, inplace=True)
join['is35'] = ((join.agepreg >= 35) * (join.agepreg < 36)).astype(int)
formula = 'numbabes ~ C(race) + age_r + C(havedeg) + C(totinc)'
model = smf.poisson(formula, data=join)
results = model.fit()
results.summary()

columns = ['race', 'age_r', 'havedeg', 'totinc']
new = pandas.DataFrame([[1, 35, 1, 14]], columns=columns)
y = results.predict(new)
y


"""
Explanation: Exercise 11.3
Use a Poisson Regression to predict how many children a woman has born if she is 35 years old, black, college graduate, with household income greater than 75000
End of explanation
"""


resp = chap01soln.ReadFemResp()

formula = 'rmarital ~ age_r + C(race) + C(educat) + C(totincr)'
model = smf.mnlogit(formula, data=resp)
results = model.fit()
results.summary()

columns = ['age_r','race','educat','totincr']
new = pandas.DataFrame([[25, 2, 12, 11]], columns=columns)
y = results.predict(new)
y

"""
Explanation: Exercise 11.4
use multinomial logistic regression to determine the categorical marital status of a 25 year old white high school graduate whose annual income is about $45,000
End of explanation
"""


columns = ['age_r','race','educat','totincr']
new = pandas.DataFrame([[43, 2, 16, 14]], columns=columns)
y = results.predict(new)
y

"""
Explanation: 43% are currently married
14% are not married but currently living with opp sex partner
An extremely small amount are widowed
2.6% are divorced
1.7% are separated due to marital discord
38% have never been maried.
End of explanation
"""


