# Illustration of datascience Tables for multivariate analysis

David E. Culler

This notebook illustrates some of the use of datascience Tables to perform regressions on multiple variables. In doing so, it shows some of the elegant ways that computational concepts and statistical concepts naturally fit together.

``````

In [108]:

# HIDDEN
from datascience import *
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plots
plots.style.use('fivethirtyeight')

``````

## Read data from a local or remote csv file to create a Table

``````

In [109]:

# Source https://onlinecourses.science.psu.edu/stat501/node/380
births=Table.read_table("data/birthsmokers.csv")
births

``````
``````

Out[109]:

Wgt Gest Smoke

2940 38   1

3130 38   0

2420 36   1

2450 34   0

2760 39   1

2440 35   1

3226 40   0

3301 42   1

2729 37   0

3410 40   0

... (22 rows omitted)

``````

## Looking at the raw data

`Table.scatter` produces a scatter plot of columns versus a specific column. Here we look at how birthweight varies with gestation. And we note whether the mother smoked.

``````

In [161]:

# First let's just look at what is here.  This needs to be a scatter, rather
# than a plot because there is no simple ordering - just relationships between
# birthweight and gestation time along with whether the mother smokes.
births.scatter('Gest')

``````
``````

``````
``````

In [162]:

births.num_rows

``````
``````

Out[162]:

32

``````
``````

In [163]:

# How many samples in each category
births.where('Smoke').num_rows

``````
``````

Out[163]:

16

``````

## Fitting a line to the data on a scatter plot

Here we drop the `smoke` column and look at the birthweight for the whole population.

``````

In [114]:

# As there is a trend among birthweight and gestation, we can show a fit line to try to
# capture it
births.drop('Smoke').scatter('Gest', fit_line=True)

``````
``````

``````

## Partitioning data

The question is whether smoking causes the trend to be substantially different. Split the data into two tables using the `Smoke` column. Then we can find the trends for each.

``````

In [164]:

nosmoke=births.where('Smoke',0).drop('Smoke')
smoke=births.where('Smoke',1).drop('Smoke')

``````
``````

In [165]:

# we can attempt to find the trend for each
nosmoke.scatter('Gest',fit_line=True)

``````
``````

``````
``````

In [166]:

smoke.scatter('Gest',fit_line=True, marker='x')

``````
``````

``````

## Build a model by fitting a line to data in a Table

Selecting a column of a Table yields a numpy array, allowing the numpy numerical tools to be used for fitting curves to the data. For documentation on the `polyfit` function, try `help(np.polyfit)`. It returns the polynomial coefficients, highest power first, which is the slope for a line.

A linear model is only meaningful in the range around the normal gestation period, and indeed the intercept is not physcially meaningful

``````

In [123]:

# what is the coeeficients of the line fitted to these data?

np.polyfit(nosmoke['Gest'],nosmoke['Wgt'],1)

``````
``````

Out[123]:

array([  147.20689655, -2546.13793103])

``````
``````

In [124]:

np.polyfit(smoke['Gest'],smoke['Wgt'],1)

``````
``````

Out[124]:

array([  139.02874903, -2474.56410256])

``````

We see that both the constant and the weight increase per week is lower for the smokers.

## Higher order functions as models

At this point, we could do mx+b all over the place, or we could utilize higher order functions to capture the concept of a model. Here is an example, building such a model directly from the data. It takes the set of `x` and `y` values and returns a function that evaluates the model built from that data at a particular `x`.

``````

In [167]:

# Build a linear model from data and return it is a function
def make_lm(x, y):
m,b = np.polyfit(x, y, 1)
def lm(a):
return m*a+b
return lm

``````
``````

In [168]:

# Create model of non-smokers that returns estimated weight a function of weeks gestation
nosmoker_weight = make_lm(nosmoke['Gest'], nosmoke['Wgt'])

``````
``````

In [169]:

# Evaluate it at normal gestations
nosmoker_weight(40)

``````
``````

Out[169]:

3342.1379310344846

``````
``````

In [170]:

smoker_weight = make_lm(smoke['Gest'], smoke['Wgt'])

``````
``````

In [171]:

smoker_weight(40)

``````
``````

Out[171]:

3086.5858585858577

``````

## Drawing a conclusion

Using the models to remove the contribution due to gestation time, we can attempt to draw a conclusion about the effect of smoking at typical gestation age.

``````

In [133]:

# based on this data set, fitting the data to models of weight as a function of gestation
# We might conclude that at 40 weeks the effect of smoking on birthweigth in grams is
smoke_diff = nosmoker_weight(40) - smoker_weight(40)
smoke_diff

``````
``````

Out[133]:

255.55207244862686

``````
``````

In [139]:

# Or in relative terms
"{0:.1f}%".format(100*(nosmoker_weight(40)-smoker_weight(40))/nosmoker_weight(40))

``````
``````

Out[139]:

'7.6%'

``````

## Use the models to build a Table and visualize the effect of the categorical parameter - smoking

``````

In [142]:

# Create a table with a column containing the independent variable
estimated_birthweight = Table().with_column('week', np.arange(32,44))
# Add columns of dependent variables
estimated_birthweight['nosmoke'] = estimated_birthweight.apply(nosmoker_weight,'week')
estimated_birthweight['smoke'] = estimated_birthweight.apply(smoker_weight,'week')
estimated_birthweight

``````
``````

Out[142]:

week nosmoke smoke

32   2164.48 1974.36

33   2311.69 2113.38

34   2458.9  2252.41

35   2606.1  2391.44

36   2753.31 2530.47

37   2900.52 2669.5

38   3047.72 2808.53

39   3194.93 2947.56

40   3342.14 3086.59

41   3489.34 3225.61

... (2 rows omitted)

``````
``````

In [143]:

# plot it to visualize
estimated_birthweight.plot('week',overlay=True)

``````
``````

``````

## Determining if the conclusion is sound

At this point, we might ask how accurately these linear models fit the data. The 'residual' in the fit would give us some idea of this. But the error in summarizing a collection of empirical data with an analytical model is only a part of the question.

The deeper point is that this data is not "truth", it is merely a sample of a population, and a tiny one at that. We should be asking, is an inference drawn from this data valid for the population that the data is intended to represent. Of course, the population is not directly observable, only the sample of it. How can we use the sample we have to get some idea of how representative it is of the larger population. That is what bootstrap seeks to accomplish.

Tables provide a method `sample` for just this purpose. Here we return to looking at the coefficients, rather than build a function for the model.

``````

In [147]:

# Construct a new model by forming a new sample from our existing one and fiting a line to that
# Here we create quite a general function, which takes a table and column names over which
# the model is to be formed.
def rboot(table, x, y):
sample = table.sample(table.num_rows, with_replacement=True)
return np.polyfit(sample[x],sample[y],1)

``````
``````

In [148]:

# Try it out for non-smokers.  Note that every time this cell is evaluated (ctrl-enter)
# the result is a little different, since a new sample is drawn.
rboot(nosmoke, 'Gest', 'Wgt')

``````
``````

Out[148]:

array([  156.52402235, -2870.49608939])

``````
``````

In [149]:

# And for smokers
rboot(smoke,'Gest','Wgt')

``````
``````

Out[149]:

array([  146.4230423 , -2788.49144914])

``````

## Bootstrap

Using this model builder as a function, draw many samples and form a model for each.
Then we can look at the distribution of the model parameters over lots of models.

This illustrates the construction of tables by rows. The Table constructor accepts the column names and `with_rows` fills them in row by row.

`hist` forms and shows a histogram of the result.

``````

In [48]:

# Bootstrap a distribution of models by drawing many random samples, with replacement, from our samples
num_samples = 1000
nosamples = Table(['slope','intercept']).with_rows([rboot(nosmoke,'Gest','Wgt') for i in range(num_samples)])
nosamples.hist(bins=50,normed=True, overlay=False)

``````
``````

``````

And we repeate this for the other category.

``````

In [50]:

smokesamples = Table(['slope','intercept']).with_rows([rboot(smoke,'Gest','Wgt') for i in range(num_samples)])
smokesamples.hist(bins=50,normed=True, overlay=False)

``````
``````

``````

## Summary of sample distributions of the regression

At this point we could compute a statistic over the sample distributions of these parameters, such as the total variational distance, or the mean.

``````

In [150]:

nosamples.stats([np.min,np.mean,np.max])

``````
``````

Out[150]:

statistic slope intercept

amin      115.09  -4062.35

mean      148.435 -2590.19

amax      185.417 -1277.11

``````
``````

In [151]:

smokesamples.stats([np.min,np.mean,np.max])

``````
``````

Out[151]:

statistic slope intercept

amin      99.3578 -4622.06

mean      140.249 -2525.15

amax      191.544 -889.25

``````

## Estimation of birthweights at 40 weeks

Selecting a column of a Table yields a numpy array. Arithmetic operators work elementwise on the entire array.

``````

In [152]:

smokesamples['slope']*40+smokesamples['intercept']

``````
``````

Out[152]:

array([ 3054.33291615,  3071.45874587,  3104.57205882,  3115.62472885,
3095.25548589,  3095.76451819,  3115.8479307 ,  3082.71710526,
3054.5940348 ,  3061.26005362,  3081.67708333,  3118.06451613,
3094.25158151,  3066.21015011,  3089.47261663,  3091.34962406,
3031.62903226,  3086.93859649,  3045.14313725,  3063.68311404,
3072.21691974,  3032.4365942 ,  3045.31818182,  3074.80638723,
3102.66857143,  3070.49838188,  3069.46820405,  3086.96212121,
3063.79456706,  3111.2038835 ,  3116.02455357,  3102.54857143,
3021.59727385,  3148.7797619 ,  3094.43887147,  3057.30898405,
3042.54303279,  3076.22222222,  3071.17910448,  3086.91304348,
3040.26369168,  3099.42772612,  3095.79840637,  3117.00611621,
3036.45175439,  3133.84635417,  3123.0327381 ,  3139.11567164,
3073.41779789,  3066.02380952,  3071.76666667,  3068.26756497,
3071.83107089,  3069.05555556,  3077.48526523,  3086.85201794,
3114.05018182,  3138.35339506,  3098.83704974,  3095.66269841,
3095.42741935,  3029.92070485,  3051.53529412,  3084.58501119,
3105.21212121,  3069.44384058,  3095.28247423,  3060.33688699,
3074.01552511,  3097.47475509,  3099.03537532,  3075.86286286,
3125.93452381,  3116.24717408,  3025.55135773,  3103.32307692,
3082.68850073,  3104.11232877,  3130.34501718,  3111.34080299,
3092.02974628,  3098.58165548,  3078.84292237,  3095.33763838,
3045.89331437,  3036.39846743,  3122.21871345,  3072.37058824,
3101.125     ,  3092.60416667,  3048.49683544,  3107.02461538,
3094.69815195,  3092.38541667,  3089.83333333,  3103.84782609,
3051.56006768,  3066.16339869,  3112.77272727,  3085.66497462,
3078.0518688 ,  3090.59661017,  3085.75801749,  3095.30974978,
3050.34246575,  3092.51082251,  3112.328     ,  3037.5       ,
3084.27427184,  3065.27212242,  3110.75409836,  3121.22179239,
3059.29850746,  3077.42627346,  3076.875     ,  3063.46240602,
3073.16845878,  3097.12      ,  3043.77077364,  3101.98830409,
3067.67986799,  3088.34679335,  3068.67491749,  3115.02428941,
3093.22622623,  3042.5070922 ,  3054.51319648,  3048.78679505,
3083.35326587,  3083.33062551,  3089.75217736,  3067.88410256,
3080.28636364,  3071.77777778,  3148.98433048,  3108.99775785,
3055.06349206,  3093.83806818,  3051.21521522,  3060.04681195,
3067.33990148,  3086.21011947,  3051.31355932,  3077.08185986,
3098.75793651,  3134.30250784,  3049.04105572,  3052.14347357,
3081.22114609,  3084.79069767,  3062.27319588,  3079.78726592,
3080.34330387,  3122.48078267,  3095.35267349,  3082.125     ,
3078.61528822,  3102.85025641,  3095.43654822,  3047.06333333,
3108.50470219,  3056.21225071,  3085.20833333,  3074.15942029,
3123.45277778,  3077.2962963 ,  3121.66754444,  3090.57692308,
3057.88329519,  3048.53697383,  3071.63690476,  3077.04202483,
3093.14109091,  3057.45842217,  3028.25799574,  3094.61610113,
3011.18421053,  3091.65428571,  3074.92250923,  3097.74358974,
3090.03590426,  3043.25416667,  3085.53630363,  3088.11904762,
3088.06329114,  3055.97818599,  3073.8125    ,  3098.79166667,
3080.44949495,  3028.16017316,  3119.99385818,  3150.53030303,
3075.0625    ,  3065.07489879,  3094.25      ,  3080.72916667,
3121.74307116,  3073.94923858,  3104.32118165,  3089.5375    ,
3110.87990196,  3092.38784629,  3121.59288538,  3043.291364  ,
3051.18956044,  3094.0177683 ,  3033.57780196,  3093.14422369,
3120.41433779,  3096.20048309,  3086.65079365,  3092.29909526,
3089.0458613 ,  3046.08614232,  3112.1264214 ,  3130.39835165,
3112.03296703,  3073.71580123,  3064.17857143,  3074.39      ,
3057.14666667,  3104.7220339 ,  3116.10394265,  3073.89277389,
3115.97947214,  3088.34601449,  3055.19301587,  3078.4375    ,
3118.07333333,  3067.59936909,  3138.63809524,  3124.07828283,
3061.31111111,  3047.04209329,  3054.23484267,  3097.25498575,
3062.25396825,  3073.42307692,  3066.89119171,  3050.09861111,
3092.38721461,  3068.2411985 ,  3100.37808642,  3093.7293578 ,
3048.71493902,  3052.39200864,  3087.28382838,  3056.83552056,
3078.58666667,  3085.27493606,  3096.54785479,  3062.99780702,
3066.83093829,  3117.95100612,  3094.20588235,  3079.00367985,
3110.69558824,  3086.37820513,  3098.48133199,  3073.17877095,
3089.13675214,  3092.11607143,  3125.9793578 ,  3092.10566357,
3051.97790869,  3082.23330516,  3118.25252525,  3076.03846154,
3120.13882532,  3057.1300813 ,  3123.75247525,  3108.57260274,
3106.82226212,  3099.0421187 ,  3097.88202692,  3086.67559524,
3080.58006873,  3093.14224446,  3059.4103139 ,  3073.3844086 ,
3100.97260274,  3032.41666667,  3083.65441176,  3098.08196721,
3099.73433584,  3009.21684588,  3123.71342079,  3091.16165414,
3054.2246696 ,  3073.01932367,  3106.31363636,  3052.36959064,
3096.75132275,  3047.39908257,  3097.29035339,  3051.15351812,
3146.65720294,  3124.22532027,  3041.39310345,  3078.78968254,
3095.43389831,  3092.8797619 ,  3086.20906801,  3085.25505051,
3064.61463415,  3104.4       ,  3099.47058824,  3095.46995885,
3051.41229909,  3123.58539149,  3116.09395973,  3058.04476629,
3082.23037101,  3077.54385965,  3088.17383513,  3072.03970588,
3101.90625   ,  3058.95035461,  3097.75      ,  3095.81317204,
3060.45415596,  3069.12527964,  3120.0625    ,  3073.8937409 ,
3122.37862319,  3153.15749235,  3098.3245614 ,  3101.002002  ,
3054.01515152,  3133.125     ,  3086.90066858,  3044.11232877,
3089.375     ,  3131.62931034,  3105.88228005,  3105.86435569,
3148.51645721,  3126.94680851,  3076.2259887 ,  3053.83673469,
3094.43654822,  3079.21003135,  3026.72352941,  3064.92156863,
3075.26315789,  3109.78240741,  3088.12881356,  3073.49434815,
3072.93916024,  3107.92862625,  3129.43977591,  3076.19583333,
3063.75059761,  3028.77199504,  3075.37341772,  3118.36011905,
3112.56377079,  3105.37128713,  3052.72200264,  3033.10861423,
3053.29113924,  3095.55835432,  3075.77548387,  3082.81345566,
3092.62475697,  3050.69827586,  3083.12445278,  3063.21251476,
3097.48396501,  3136.47433116,  3097.06269113,  3064.38627451,
3110.72336066,  3044.97339593,  3085.61842105,  3113.44666667,
3088.80614657,  3090.53942652,  3093.03559322,  3084.        ,
3035.50462963,  3096.69240196,  3088.51331719,  3139.67991521,
3123.00676133,  3113.48214286,  3076.24128686,  3083.92088316,
3072.62724014,  3047.29258903,  3115.99337748,  3090.6285495 ,
3087.27734171,  3064.62881356,  3038.22607261,  3062.72607261,
3112.36943907,  3091.37829691,  3109.78052434,  3090.56798246,
3105.40252828,  3104.004313  ,  3081.70562771,  3113.7403599 ,
3097.35098039,  3142.85984848,  3061.11372549,  3068.96873121,
3116.89473684,  3105.35574668,  3045.6941896 ,  3083.72376543,
3082.2020202 ,  3094.98878205,  3134.9623312 ,  3143.34166667,
3104.61840929,  3099.46530148,  3074.84908537,  3110.75      ,
3079.80079681,  3111.26666667,  3085.54768154,  3059.56402737,
3101.79276637,  3057.2899106 ,  3054.13492063,  3125.73833671,
3107.83179012,  3042.79615048,  3050.25539258,  3093.92171717,
3075.8989899 ,  3063.61066667,  3060.04027421,  3140.8134058 ,
3082.25373134,  3108.36458333,  3118.89640365,  3092.1059761 ,
3055.25793651,  3070.50666667,  3058.83333333,  3041.53976765,
3071.66440678,  3073.41604351,  3114.76130389,  3057.67207573,
3061.87147335,  3121.46778351,  3115.65853659,  3062.25625423,
3076.01344086,  3049.05128205,  3101.9672619 ,  3138.18857143,
3091.80401235,  3077.94690265,  3085.49258649,  3067.085963  ,
3103.23771224,  3159.98699764,  3049.14912281,  3079.94927536,
3083.81761006,  3036.21825397,  3086.09388336,  3169.27653471,
3071.10762332,  3039.66580977,  3055.65287588,  3095.5767098 ,
3094.92114208,  3113.93108974,  3093.99679487,  3111.94206715,
3109.54905336,  3040.30676085,  3091.68275862,  3100.38155803,
3093.64578313,  3104.11711712,  3078.68888889,  3110.10164271,
3115.29108245,  3083.30218538,  3086.50531582,  3080.37878788,
3065.22079772,  3069.10586177,  3083.19520548,  3073.85      ,
3037.87878788,  3090.27705295,  3074.58823529,  3075.14114514,
3081.8046595 ,  3033.46872461,  3073.39933993,  3089.30661809,
3057.29292929,  3093.82352941,  3030.73058252,  3143.28305085,
3058.14779874,  3104.91385768,  3092.1875    ,  3080.53677173,
3059.09251101,  3084.72674419,  3057.19179487,  3043.84686347,
3092.52580195,  3071.45333333,  3091.03921569,  3077.58580858,
3079.64391951,  3128.18678161,  3120.82200647,  3076.30344828,
3101.87339972,  3098.08195849,  3097.74143955,  3096.90849243,
3070.87641296,  3013.20055198,  3112.77333333,  3103.22941176,
3137.65553373,  3092.86850153,  3139.01591814,  3072.90181181,
3074.12334802,  3039.4375    ,  3065.56372549,  3076.84361792,
3098.02666667,  3104.91880342,  3099.5448592 ,  3122.21141975,
3098.36561743,  3099.0330912 ,  3113.05090138,  3107.59311741,
3066.4290429 ,  3088.25696594,  3084.99059266,  3036.5962963 ,
3108.7428088 ,  3130.34090909,  3106.11041009,  3138.79903148,
3131.09574468,  3067.95643939,  3087.80952381,  3144.83911672,
3130.38596491,  3075.46799117,  3075.4       ,  3082.18181818,
3020.54666667,  3139.38967742,  3047.36092715,  3016.75782537,
3104.84518349,  3057.87991927,  3082.88555347,  3079.04701931,
3041.42794118,  3090.2479564 ,  3061.67164179,  3096.24354244,
3029.05555556,  3105.34299517,  3008.93670886,  3095.275     ,
3079.45745746,  3101.84782609,  3117.45138889,  3142.38523852,
3178.70550162,  3061.64280332,  3125.41269841,  3091.7699115 ,
3081.3       ,  3105.03013699,  3100.02366464,  3090.51920694,
3086.24248777,  3065.60206961,  3088.93189964,  3059.33796296,
3052.91928251,  3043.68504772,  3053.72294705,  3001.44235925,
3078.60178571,  3087.36210317,  3034.99820789,  3097.37990196,
3081.38793879,  3109.58928571,  3082.79627715,  3069.52970297,
3045.60252366,  3086.84115226,  3129.69179785,  3047.50825083,
3066.38116592,  3089.80555556,  3092.18128655,  3047.2693299 ,
3105.36531082,  3070.07590759,  3069.10144928,  3061.9650594 ,
3131.1684435 ,  3065.59259259,  3013.9562109 ,  3079.56819694,
3053.80294118,  3103.31958763,  3092.31018519,  3088.88528139,
3079.81944444,  3078.58095238,  3098.04949054,  3083.4040404 ,
3071.61677513,  3116.41917808,  3099.07681366,  3085.06372549,
3117.75179856,  3051.92888889,  3077.4256927 ,  3147.65949821,
3133.85390428,  3105.997151  ,  3111.20524345,  3088.4267844 ,
3105.5327381 ,  3055.27727273,  3103.47619048,  3040.19634703,
3087.88095238,  3166.93721973,  3125.48809524,  3091.07862759,
3015.05946792,  3091.96873156,  3086.91304348,  3082.58333333,
3081.875     ,  3080.18435754,  3131.04371585,  3076.1589437 ,
3088.08725603,  3058.31015143,  3092.44214047,  3084.78840125,
3074.2601476 ,  3107.61730899,  3062.90419636,  3108.1954023 ,
3066.33241758,  3112.45694532,  3085.59649123,  3072.60648148,
3077.3925328 ,  3027.68912416,  3077.12113402,  3110.39491691,
3072.71155683,  3065.62465501,  3045.62587904,  3116.65800866,
3033.25      ,  3069.43529412,  3100.55987055,  3038.86171761,
3111.15209666,  3082.58868502,  3052.96140351,  3054.94642857,
3113.69547872,  3126.25612849,  3052.33333333,  3095.58328096,
3088.85826002,  3082.14370245,  3097.23809524,  3070.66013072,
3087.85416667,  3090.7405303 ,  3072.65307533,  3092.25      ,
3164.04347826,  3096.7086758 ,  3056.67166213,  3115.6116208 ,
3065.20661824,  3073.11396011,  3062.86883273,  3091.79310345,
3098.16666667,  3086.43824701,  3055.32279315,  3124.9375    ,
3072.89027778,  3107.30555556,  3061.86065574,  3049.86184211,
3060.72222222,  3103.19468025,  3123.3131068 ,  3104.62390671,
3087.92473118,  3056.34266667,  3092.9375    ,  3067.62400578,
3038.52439863,  3109.1712963 ,  3061.6082397 ,  3074.15686275,
3097.32169866,  3127.87235841,  3149.3125    ,  3084.4880137 ,
3064.4422043 ,  3087.32379438,  3079.15614237,  3098.25635739,
3105.47732181,  3129.43626287,  3038.32932933,  3014.77083333,
3109.37885463,  3135.58983891,  3074.        ,  3094.59920635,
3096.90334044,  3049.50971251,  3097.33673469,  3113.35550277,
3031.06324111,  3124.025     ,  3028.14545455,  3111.66906044,
3060.93525809,  3076.6875    ,  3055.18664384,  3092.4969697 ,
3053.3172043 ,  3108.93026941,  3086.83236152,  3052.94897959,
3102.62037037,  3125.42622951,  3057.65277778,  3134.95673077,
3092.71308833,  3103.57027226,  3046.67489712,  3038.48101266,
3115.81520967,  3073.25      ,  3108.94927536,  3112.90505051,
3070.3614304 ,  3103.07083333,  3107.52807392,  3089.79772543,
3142.13043478,  3078.87118644,  3089.2689769 ,  3063.35277622,
3075.08380952,  3064.98230088,  3109.13841599,  3100.72570194,
3058.22229776,  3061.0880069 ,  3106.75694444,  3097.47058824,
3063.48110831,  3080.44338876,  3105.88291139,  3056.57058824,
3033.97759104,  3107.35820896,  3114.40551181,  3074.80896027,
3087.73672727,  3092.92982456,  3089.3064877 ,  3135.28402367,
3093.76767677,  3044.42156863,  3074.09744299,  3071.98790323,
3091.20430108,  3105.79733333,  3044.91333333,  3115.92056583,
3097.09722222,  3085.92152466,  3079.04892966,  3041.17838638,
3039.72222222,  3077.79385045,  3095.62660256,  3090.51646284,
3070.0929878 ,  3108.20138889,  3072.81899642,  3095.34099617,
3096.05592347,  3160.05606901,  3049.59951338,  3079.83137255,
3121.675     ,  3125.63833897,  3103.88965517,  3040.01029412,
3092.99305556,  3070.39393939,  3097.14814815,  3110.640553  ,
3099.52009185,  3102.91748527,  3108.27090779,  3085.54297694,
3085.25925926,  3087.53272577,  3087.68509213,  3125.43236715,
3078.96078431,  3090.88538682,  3099.47619048,  3046.9057971 ,
3087.6092233 ,  3115.10514665,  3072.0315534 ,  3078.375     ,
3114.91435768,  3083.37525355,  3047.94871795,  3074.11882229,
3082.78333333,  3121.54104478,  3087.25520833,  3106.64975845,
3099.06837607,  3065.3016055 ,  3114.9235474 ,  3097.09452736,
3076.44736842,  3095.15942029,  3095.01619433,  3075.97055731,
3070.37647059,  3105.6429137 ,  3070.87356322,  3080.24216524,
3063.81481481,  3070.91338583,  3094.43804756,  3055.57230769,
3074.96485788,  3142.40212633,  3099.00414938,  3049.75      ,
3112.27419355,  3057.78794813,  3094.95707763,  3103.63333333,
3065.4024864 ,  3090.4518028 ,  3111.90065829,  3095.64156946,
3113.38325991,  3036.50048123,  3080.83043478,  3098.08012821,
3066.83152521,  3062.90100251,  3095.5212766 ,  3145.64106583,
3096.37004405,  3065.21461187,  3109.56743536,  3073.15154185,
3077.81432896,  3132.3597786 ,  3046.33333333,  3061.67253177,
3057.23809524,  3049.94794521,  3079.78356164,  3070.62271415,
3042.68290909,  3090.24284764,  3091.98684211,  3077.7103321 ,
3047.44658325,  3070.66745843,  3065.09455587,  3070.69303797,
3076.38778878,  3082.43612335,  3097.04545455,  3020.02888087,
3104.19362745,  3070.75632911,  2984.13891835,  3113.49737119,
3074.33855799,  3072.77990868,  3087.82571912,  3107.77578475,
3086.10187266,  3019.45416667,  3084.47619048,  3076.34935897,
3092.38608807,  3089.19978632,  3137.95014245,  3090.74085366,
3070.90279115,  3061.36479401,  3073.21630094,  3064.820059  ,
3082.01650165,  3101.90305206,  3103.16618076,  3067.44736842,
3093.48448687,  3100.0426009 ,  3060.94736842,  3094.4532967 ,
3102.02542373,  3066.55987055,  3083.31439394,  3053.80020182,
3063.48748749,  3107.71707317,  3052.6152019 ,  3091.13064133,
3099.27118644,  3104.07590759,  3084.35887612,  3075.1875    ,
3072.19444444,  3070.51587302,  3117.87916667,  3124.34287867,
3112.02956989,  3134.95989123,  3046.31066667,  3124.24050633,
3114.67745415,  3070.9607698 ,  3063.40233918,  3074.66396761,
3069.69349845,  3142.66666667,  3101.49404762,  3060.69926199,
3069.37239583,  3028.50785669,  3121.02130493,  3051.21509824,
3071.58441558,  3112.01029412,  3099.96448087,  3040.37115839,
3084.01619433,  3076.88405797,  3086.20803783,  3106.12060302,
3083.26649529,  3067.30188679,  3133.35579196,  3069.00682594,
3093.83828611,  3118.20821662,  3058.83817427,  3094.01084337,
3118.41393787,  3077.51979566,  3028.18253968,  3082.07142857,
3148.46401515,  3023.09221311,  3135.7491814 ,  3042.125     ,
3082.97697698,  3058.27013631,  3059.56076759,  3099.89900427])

``````
``````

In [153]:

# So now we have an estimate of the distribution of birthweights at week 40 for
# something closer to the populations that these small samples represent.
weights_40 = Table().with_columns([
('nosmoke', nosamples['slope']*40+nosamples['intercept']),
('smoke', smokesamples['slope']*40+smokesamples['intercept'])])
weights_40

``````
``````

Out[153]:

nosmoke smoke

3368.64 3054.33

3345.99 3071.46

3378.07 3104.57

3324.19 3115.62

3350.16 3095.26

3362.7  3095.76

3330.87 3115.85

3330.15 3082.72

3376.09 3054.59

3354.27 3061.26

... (990 rows omitted)

``````
``````

In [154]:

weights_40['Smoke Wgt Loss'] = weights_40['nosmoke'] - weights_40['smoke']

``````
``````

In [155]:

# what do we expect the distribution of birthweight reduction due to smoking to look like
# for the population represented by the original sample?
weights_40.select('Smoke Wgt Loss').hist(bins=30,normed=True)

``````
``````

``````
``````

In [156]:

smoke_diff

``````
``````

Out[156]:

255.55207244862686

``````
``````

In [157]:

def firstQtile(x) : return np.percentile(x,25)
def thirdQtile(x) : return np.percentile(x,25)
summary_ops = (min, firstQtile, np.median, np.mean, thirdQtile, max)

``````
``````

In [158]:

summary = weights_40.stats(summary_ops)
summary

``````
``````

Out[158]:

statistic nosmoke smoke Smoke Wgt Loss

min        3250.37 2984.14 141.919

firstQtile 3324.78 3066.42 232.572

median     3347.94 3085.71 263.573

mean       3347.23 3084.8  262.435

thirdQtile 3324.78 3066.42 232.572

max        3456.79 3178.71 421.632

``````
``````

In [159]:

summary['diff']=summary['nosmoke']-summary['smoke']

``````
``````

In [160]:

# the bottom line
summary

``````
``````

Out[160]:

statistic nosmoke smoke Smoke Wgt Loss diff

min        3250.37 2984.14 141.919        266.234

firstQtile 3324.78 3066.42 232.572        258.358

median     3347.94 3085.71 263.573        262.231

mean       3347.23 3084.8  262.435        262.435

thirdQtile 3324.78 3066.42 232.572        258.358

max        3456.79 3178.71 421.632        278.083

``````

## Visualizing the separation of these distributions

``````

In [67]:

weights_40.select(['smoke','nosmoke']).hist(overlay=True,bins=30,normed=True)

``````
``````

``````

# Empirical p-values

A more formal approach would be to take as the null hypothesis that smoking does not affect the birthweight. Repeatedly split the data into random halves and model the birthweight difference. What is the chance that the difference we see in summary table is an artifact?

``````

In [78]:

# As an example, split the original data into two random halves
A, B = births.split(births.num_rows//2)

``````
``````

In [76]:

A

``````
``````

Out[76]:

Wgt Gest Smoke

3500 42   1

2729 37   0

2940 38   1

3530 42   0

2520 35   0

2715 36   1

3244 39   0

3459 40   0

3130 38   0

2957 39   1

... (6 rows omitted)

``````
``````

In [77]:

B

``````
``````

Out[77]:

Wgt Gest Smoke

2740 38   1

2920 38   0

3446 42   1

3095 39   0

2841 36   0

3322 39   0

3175 41   1

2619 35   0

3130 39   1

2420 36   1

... (6 rows omitted)

``````
``````

In [84]:

make_lm(A['Gest'], A['Wgt'])(40)

``````
``````

Out[84]:

3193.9845474613712

``````
``````

In [88]:

make_lm(B['Gest'], B['Wgt'])(40)

``````
``````

Out[88]:

3211.4556962025335

``````

## Capturing statistical computations and general tools

Rather than compute the null hypothesis for this particular table, we can build a very general tool, as a function, that will do it in general.

Then we can use it to build a sample distribution under the null hypothesis

``````

In [104]:

def null_diff_at(tbl, x, y, w):
A, B = tbl.split(tbl.num_rows//2)
return make_lm(A[x], A[y])(w) - make_lm(B[x], B[y])(w)

``````
``````

In [105]:

null_diff_at(births, 'Gest', 'Wgt', 40)

``````
``````

Out[105]:

77.092406857131209

``````
``````

In [106]:

null = Table().with_column('Diff', [null_diff_at(births, 'Gest', 'Wgt', 40) for i in range(1000)])
null.hist()

``````
``````

``````

What is the probablility that we got a 260g difference in birthweight at 40 weeks as an artifact of the sample?

Zero

``````

In [107]:

null.stats()

``````
``````

Out[107]:

statistic Diff

min       -205.91

max       221.76

median    -2.18387

sum       -1092.65

``````
``````

In [ ]:

``````