Introduction to Data Science <br> Winter 2018

# Homework 6
Due Saturday, 2/19 , 3:00am

Please complete this notebook by filling in the cells provided. Once you are finished, "Save and Checkpoint" the notebook, download the ipynb and html versions to your computer, and submit the html version. 

**Each question is worth 10 points** <br>
**Total points: 60**

In this homework you will use a permutation tests to analyze data from a clinical trial. 

A summary of the methods to Table is [here](http://data8.org/datascience/tables.html) <br>
A tutorial for the datascience module is [here](http://data8.org/datascience/tutorial.html) <br>
A cheatsheet for the datascience module is [here](https://github.com/wstuetzle/STAT180/blob/master/Computing/data8_sp17_midterm_ref_sheet.pdf)


```python
from datascience import *
import numpy as np
import seaborn as sns
sns.set_style("whitegrid")
%matplotlib inline
```

### Does Rogain promote hair growth ?

You will analyze the results of a randomized, double blind clinical trial to evaluate the effectiveness of Rogaine in promoting hair growth. There were 1431 subjects, of whom 714 were treated with Rogaine and 717 received a placebo. Here is a table summarizing the outcome of the experiment:

under the assumptio nthat treatment has no effect.
Is this unusual table or not!?
if unusual, threatment has effect! (contradict)

to measure how unusual my table is, if no effct, count is
same for all!!
distance between tables?
sum all the element of the table, 
(1-3) ^2 / 3

last step, is the difference large or small?d


```python
observed_counts = Table.read_table("rogaine_observed_counts.csv")
observed_counts
```




<table border="1" class="dataframe">
    <thead>
        <tr>
            <th>treatment</th> <th>dense_growth</th> <th>minimal_growth</th> <th>moderate_growth</th> <th>new_vellus</th> <th>no_growth</th>
        </tr>
    </thead>
    <tbody>
        <tr>
            <td>Placebo  </td> <td>1           </td> <td>114           </td> <td>29             </td> <td>150       </td> <td>423      </td>
        </tr>
    </tbody>
        <tr>
            <td>Rogaine  </td> <td>5           </td> <td>178           </td> <td>58             </td> <td>172       </td> <td>301      </td>
        </tr>
    </tbody>
</table>



It does look like subjects who received Rogaine fared better than subjects who received a placebo. You want to find out whether there really is an association between treatment and outcome or whether the difference could just have been due to luck of the draw in the random assignment of subjects to Rogaine and Placebo. 

** Here is the basic idea: **

There were 1431 subjects, of whom 714 (let's call that half) were treated with Rogaine.

Consider the "no_growth" column. There were 724 subjects with outcome "no_growth". If there was no treatment effect you would expect them to be split 50/50 between "Placebo" and "Rogaine", just like the overall pool of subjects. So under the null hypotheses of no treatment effect you would expect the "no_growth" column to be (362, 362).

The same argument can be made for the other four columns.

** Question 1. ** Write Python code to compute the table of expected counts (under the hypotheses of no treatment effect). Call this table "expected_counts". Print the table.

** Note: ** If this takes you too long and you run out of patience, you can compute the table of expected counts "by hand".


```python
n_placebo = np.sum(make_array(observed_counts.drop("treatment").row(0)))
n_treatment = np.sum(make_array(observed_counts.drop("treatment").row(1)))
n_total = n_placebo + n_treatment
expected_counts = Table()
for label in observed_counts.labels:
    if label == "treatment":
        expected_counts = expected_counts.with_column(label, observed_counts.column(label))
    else:
        col_sum = observed_counts.column(label).sum()
        expected_counts = expected_counts.with_column(label, (col_sum/n_total) * make_array(n_placebo, n_treatment))
expected_counts
```




<table border="1" class="dataframe">
    <thead>
        <tr>
            <th>treatment</th> <th>dense_growth</th> <th>minimal_growth</th> <th>moderate_growth</th> <th>new_vellus</th> <th>no_growth</th>
        </tr>
    </thead>
    <tbody>
        <tr>
            <td>Placebo  </td> <td>3.00629     </td> <td>146.306       </td> <td>43.5912        </td> <td>161.338   </td> <td>362.759  </td>
        </tr>
    </tbody>
        <tr>
            <td>Rogaine  </td> <td>2.99371     </td> <td>145.694       </td> <td>43.4088        </td> <td>160.662   </td> <td>361.241  </td>
        </tr>
    </tbody>
</table>



Now we need a measure of distance between the observed and expected tables. A commonly used measure is the Chi-square distance, defined as 
\begin{equation}
Chi2 = \sum_{rows} \sum_{columns} \frac{(observed\_count - expected\_count)^2}{expected\_count}
\end{equation}

** Question 2. ** Write a Python function  Chi_square_distance(observed_counts, expected_counts)  that computes the Chi-square distance between a table of observed counts and a table of expected counts. Compute and print the distance.


```python
def Chi_square_distance(observed_counts, expected_counts):
    chi2 = 0
    for i in np.arange(1,observed_counts.num_columns):
        diff_square = (observed_counts[i]- expected_counts[i])**2
        x = sum(diff_square / expected_counts[i])
        chi2 =  x + chi2
    return chi2
    
    
observed_chi_square = Chi_square_distance(observed_counts, expected_counts)
observed_chi_square
```




    48.41577071466682



Next you need to randomly permute the treatment labels, thereby destroying any association between treatment and outcome, and find out what values of the Chi-square distance you would obtain under the null hypotheses of "no association" (or equivalently "no treatment effect")

To do that you first need to make up the raw_data table that gave rise to the observed_counts. The raw_data table will have 1,431 rows (one for each subject), and two columns, "treatment" and "outcome". The treatment column has entries "Rogaine" and "Placebo". The outcome column has entries "no_growth", "new_vellus", ... , "dense_growth"

** Question 3. ** Write Python code to create a raw_data table that such that raw_data.pivot("outcome", "treatment") 
produces the observed counts.

** Hints: **

np.full(714, "Rogaine") produces an array of length 714 filled with the string "Rogaine" <br>
np.concatenate([array1, array2, array3]) glues any number of arrays together into a single array


```python
z = []
for j in [0,1]:
    table = observed_counts.take(j)
    for i in np.arange(1,observed_counts.num_columns):
        z = np.concatenate([z, np.full(table.column(i)[0], table.labels[i])])
z
```




    array(['dense_growth', 'minimal_growth', 'minimal_growth', ...,
           'no_growth', 'no_growth', 'no_growth'],
          dtype='<U32')




```python
raw_data = Table().with_columns(
    "treatment", np.concatenate([np.full(717, "Placebo"), np.full(714, "Rogaine")]),
    "outcome", z
)

raw_data.pivot("outcome", "treatment")
```




<table border="1" class="dataframe">
    <thead>
        <tr>
            <th>treatment</th> <th>dense_growth</th> <th>minimal_growth</th> <th>moderate_growth</th> <th>new_vellus</th> <th>no_growth</th>
        </tr>
    </thead>
    <tbody>
        <tr>
            <td>Placebo  </td> <td>1           </td> <td>114           </td> <td>29             </td> <td>150       </td> <td>423      </td>
        </tr>
    </tbody>
        <tr>
            <td>Rogaine  </td> <td>5           </td> <td>178           </td> <td>58             </td> <td>172       </td> <td>301      </td>
        </tr>
    </tbody>
</table>



Now you are ready for the permutation test:

* generate 1,000 permuted_raw_data tables by randomly permuting the "treatment" column of raw_data
* for each permuted_raw_data table compute the table of counts - let's call it permuted_counts.
* for each permuted_counts table compute permuted_chi_square = Chi_square_distance(permuted_counts, expected_counts)

You end up with 1,000 values for permuted_chi_square. Put them in a table permuted_chi_square_table with one column named "permuted_chi_square".

** Question 4. ** Draw a histogram of the 1,000 values of permuted_chi_square


```python
num_sample = 1000


permuted_chi_square_array = np.zeros(num_sample)

for i in range(num_sample):
    shuffled_treatment = raw_data.sample(n_total, with_replacement=False).column("treatment")
    permuted_raw_data = raw_data.with_column("treatment", shuffled_treatment)
    permuted_count = permuted_raw_data.pivot("outcome", "treatment")
    permuted_chi_square = Chi_square_distance(permuted_count, expected_counts)
    permuted_chi_square_array[i] = permuted_chi_square
    
permuted_chi_square_table = Table().with_columns("permuted_chi_square", permuted_chi_square_array)
```


```python
max_pcs = np.max(permuted_chi_square_table.column("permuted_chi_square"))
permuted_chi_square_table.hist("permuted_chi_square", bins = np.arange(0, max_pcs, 0.5))
```

** Question 5. ** Compute the p_value of the test  (See comment below)


```python
p_value = permuted_chi_square_table.where("permuted_chi_square", are.above(observed_chi_square)).num_rows / num_sample
p_value
# It should be 1/1001
```




    0.0



** Question 6. ** Here are two alternative explanation for the distance between observed and expected counts:

(a) Rogaine has no effect. The manufacturers of Rogaine were lucky: Rogaine just happened to be assigned to subjects whose hair would have improved no matter what.

(b) There is an association between treatment and outcome. Rogaine does promote hair growth.

Choose (a) or (b)


```python
b
```

** Comment on p-value ** 

It obviously can't be zero. It really is 1 / (nperm + 1)

** Reason: ** Suppose there was no treatment effect. Then we would have observed the value chi-square for nperm + 1 random permutations of the treatment labels: the "original" one from the random assignment of treatments to subjects, and nperm "reference" ones for the permutation test. The probability that the chi_square value of the original one is among the k largest of the (nperm + 1) permutations is k / (nperm + 1). In our case k = 1, so the p-value is 1 / 1001.

02eeaa25d8d62ff481a06282e54b5c3434496955
ipynb
Jupyter Notebook
H6.ipynb
ByungSuJung/Stat180
9bdcfea8b60f33f240b2e93ee407c39507359adb
MIT
H6.ipynb
ByungSuJung/Stat180
9bdcfea8b60f33f240b2e93ee407c39507359adb
MIT
H6.ipynb
ByungSuJung/Stat180
9bdcfea8b60f33f240b2e93ee407c39507359adb
MIT
Qwen/Qwen-72B
1. YES
2. YES