Python Code: ANOVA Test

# -*- coding: utf-8 -*-
"""
Created on Wed May 25 08:52:38 2016

@author: DH KIM
"""

# ANOVA

import numpy as np
import pandas as pd
import statsmodels.formula.api as smf
import seaborn as sns
import matplotlib.pyplot as plt

data = pd.read_csv('nesarc_pds.csv', low_memory=False)

#setting variables you will be working with to numeric
data['S1Q10A'] = pd.to_numeric(data['S1Q10A'], errors='coerce')
data['S3AQ3B1'] = pd.to_numeric(data['S3AQ3B1'], errors='coerce')
data['S3AQ3C1'] = pd.to_numeric(data['S3AQ3C1'], errors='coerce')
data['CHECK321'] = pd.to_numeric(data['CHECK321'], errors='coerce')
data['SEX'] = pd.to_numeric(data['SEX'], errors='coerce')

#subset data to those people who have income and have smoked in the past 12 months
sub=data[ (data['S1Q10B']!=0) & (data['CHECK321']==1)]

sub2 = sub[['SEX', 'S3AQ3B1', 'S3AQ3C1', 'S1Q10A']].copy()

sub2['S3AQ3B1']=sub2['S3AQ3B1'].replace(9, np.nan)
sub2['S3AQ3C1']=sub2['S3AQ3C1'].replace(99, np.nan)

recode1 = {1: 30, 2: 22, 3: 14, 4: 5, 5: 2.5, 6: 1}
sub2['USFREQMO']= sub2['S3AQ3B1'].map(recode1)

# mean value of personal income in last 12 month
avgIncome = sub2['S1Q10A'].mean()
LvIncome = [0, sub2['S1Q10A'].quantile(0.3), sub2['S1Q10A'].quantile(0.7), sub2['S1Q10A'].max()]

# Split into 3 groups: 'low', 'medium', 'high'
splitIntoCat = pd.cut(sub2.S1Q10A, LvIncome, labels=['low', 'medium', 'high'])
sub2['INCOMECAT'] = splitIntoCat

# num of ciga per month
sub2['NUMCIGMO_EST']=sub2['USFREQMO'] * sub2['S3AQ3C1']
sub1 = sub2[['NUMCIGMO_EST', 'INCOMECAT', 'SEX']].dropna()
sub1.SEX = sub1.SEX - 1 # 0: Male, 1: Female

# using ols function for calculating the F-statistic and associated p value
model = smf.ols(formula='NUMCIGMO_EST ~ C(INCOMECAT)', data=sub1).fit()
print (model.summary())

print ('means for smoking amount by personal income')
print(sub1.groupby('INCOMECAT').mean())

print ('standard deviations for smoking amount by personal income')
print(sub1.groupby('INCOMECAT').std())

sns.factorplot(x="INCOMECAT", y="NUMCIGMO_EST", data=sub1, kind="bar", ci=None)
plt.xlabel('Income Level')
plt.ylabel('Smoking Amount')

# Moderator: SEX
print ('means for numcigmo_est by personal income according to SEX')
print(sub1.groupby('SEX').mean())

sub_m=sub1[ sub1.SEX == 0]
sub_f=sub1[ sub1.SEX == 1]

print ('Association between income level and smoking amount for Male')
model_m = smf.ols(formula='NUMCIGMO_EST ~ C(INCOMECAT)', data=sub_m).fit()
print (model_m.summary())

print ('Association between income level and smoking amount for Female')
model_f = smf.ols(formula='NUMCIGMO_EST ~ C(INCOMECAT)', data=sub_f).fit()
print (model_f.summary())

print ('Means for smoking amount by personal income for Male')
print (sub_m.groupby('INCOMECAT').mean())

print ('Means for smoking amount by personal income for Female')
print (sub_f.groupby('INCOMECAT').mean())

sns.factorplot(x="INCOMECAT", y="NUMCIGMO_EST", data=sub_m, kind="bar", ci=None)
plt.title('Smoking Amount according to income level for Male')
plt.xlabel('Income Level')
plt.ylabel('Smoking Amount')

sns.factorplot(x="INCOMECAT", y="NUMCIGMO_EST", data=sub_f, kind="bar", ci=None)
plt.title('Smoking Amount according to income level for Female')
plt.xlabel('Income Level')
plt.ylabel('Smoking Amount')

ANOVA Test in python code and interpretation - 2

ANOVA Test & Post hoc Test - Tukey HSD

# =============================================================

import numpy as np
import pandas as pd
import statsmodels.formula.api as smf
import statsmodels.stats.multicomp as multi

data = pd.read_csv(‘nesarc_pds.csv’, low_memory=False)

#setting variables you will be working with to numeric
data['S1Q10A’] = pd.to_numeric(data['S1Q10A’], errors='coerce’)
data['S3AQ3B1’] = pd.to_numeric(data['S3AQ3B1’], errors='coerce’)
data['S3AQ3C1’] = pd.to_numeric(data['S3AQ3C1’], errors='coerce’)
data['CHECK321’] = pd.to_numeric(data['CHECK321’], errors='coerce’)

#subset data to those people who have income and have smoked in the past 12 months
sub=data[ (data['S1Q10B’]!=0) & (data['CHECK321’]==1)]

sub.loc[:,'S3AQ3B1’]=sub['S3AQ3B1’].replace(9, np.nan)
sub.loc[:,'S3AQ3C1’]=sub['S3AQ3C1’].replace(99, np.nan)

#recoding number of days smoked in the past month
recode1 = {1: 30, 2: 22, 3: 14, 4: 5, 5: 2.5, 6: 1}
sub['USFREQMO’]= sub.loc[:,'S3AQ3B1’].map(recode1)

# mean value of personal income in last 12 month
avgIncome = sub.loc[:,'S1Q10A’].mean(0)
LvIncome = [0, sub.loc[:,'S1Q10A’].quantile(0.3), sub.loc[:,'S1Q10A’].quantile(0.7), sub.loc[:,'S1Q10A’].max()]

# Split into 3 groups: 'low’, 'medium’, 'high’
splitIntoCat = pd.cut(sub.S1Q10A, LvIncome, labels=['low’, 'medium’, 'high’])
sub['INCOMECAT’] = splitIntoCat

# num of ciga per month
sub['NUMCIGMO_EST’]=sub['USFREQMO’] * sub['S3AQ3C1’]
sub1 = sub[['NUMCIGMO_EST’, 'INCOMECAT’]].dropna()

# using ols function for calculating the F-statistic and associated p value
model = smf.ols(formula='NUMCIGMO_EST ~ C(INCOMECAT)’, data=sub1)
results = model.fit()
print (results.summary())

print ('means for numcigmo_est by personal income’)
print(sub1.groupby('INCOMECAT’).mean())

print ('standard deviations for numcigmo_est by personal income’)
print(sub1.groupby('INCOMECAT’).std())

# ============ Post hoc test, Tukey HSD ===========================

mc1 = multi.MultiComparison(sub1['NUMCIGMO_EST’], sub1['INCOMECAT’])
res1 = mc1.tukeyhsd()
print(res1.summary())

ANOVA Test in python code and interpretation - 1

 I’d like to know the association btw personal income level(explanatory) and smoking amount(quantitative response). I thought there’s a certain relationship. From my experience, I’ve quit smoking since new year 2016 for raising cigarette prices surprisingly. However, I wouldn’t quit the smokes if I’ve got more paid. 

 From NESARC(’S1Q10A’), I split into 3 levels as ‘low’, 'medium’, 'high’ according to the level of personal income to make it categorical. The 'low’ level was set up to 30%, 'medium’ level was up to 70% and more than 70%(30% top of personal income) are categorized as 'high’.

 ANOVA results revealed that there’re no significant relationship btw personal income and number of cigarettes smoked, F(2, 9345)=1.130 and P-value=0.323. And the belows are means and std for no of cigarettes smoked by personal income level.

means for numcigmo_est by personal income
INCOMECAT        NUMCIGMO_EST
low                       424.175214
medium               413.702265
high                     423.927979

standard deviations for numcigmo_est by personal income
INCOMECAT       NUMCIGMO_EST           
low                      342.501714
medium              311.897224
high                    335.876155

 There’re no significant difference each others and here are the results of post hoc test. Especially, ‘meandiff’ value btw 'high’ and 'low’ is very small. Of course, there’re all 'False’ results in reject columns when the significant level(alpha) is set to 0.05.

Multiple Comparison of Means - Tukey HSD,FWER=0.05
==============================================
group1 group2 meandiff  lower    upper  reject
———————————————-
high   low    0.2472  -20.5833 21.0778 False
high  medium -10.2257 -29.6021  9.1507 False
low   medium -10.4729 -29.5268  8.5809 False
———————————————-

 If someone ask me whether this result are credible or not, I would say 'No’. There might be lots of uncontrolled fators, for one simple example, the different ciga. cost by State.