Objekt Orient

Association between prior smoking and depression

2016-02-22T10:05:00.002+01:00

In this analysis, I examine if the association between nicotine dependence and depression holds true also when there was a major smoking prior to the last 12 months reported.

From the post-hoc tests we see that the X2 value increases and the p-value shrinks when we compare infrequent smokers to frequent smokers (See tables in the Coding section). The significance is lesser than if we only compare the last 12 months, maybe because both are indicator for actual smoking. Future tests should investigate actual smoking and depression.

Appendix Code

%matplotlib inline
import pandas
import numpy
import scipy.stats
import seaborn
import matplotlib.pyplot as plt

In [33]:

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

In [41]:

data['TAB12MDX'] = data['TAB12MDX'].convert_objects(convert_numeric=True)
data['TABP12MDX'] = data['TABP12MDX'].convert_objects(convert_numeric=True)
data['CHECK321'] = data['CHECK321'].convert_objects(convert_numeric=True)
data['S3AQ3B1'] = data['S3AQ3B1'].convert_objects(convert_numeric=True)
data['S3AQ3C1'] = data['S3AQ3C1'].convert_objects(convert_numeric=True)
data['AGE'] = data['AGE'].convert_objects(convert_numeric=True)

#subset data to young adults age 18 to 25 who have smoked in the past 12 months
sub1=data[(data['AGE']>=18) & (data['AGE']<=25) & (data['CHECK321']==1)]

#make a copy of my new subsetted data
sub2 = sub1.copy()

# recode missing values to python missing (NaN)
sub2['S3AQ3B1']=sub2['S3AQ3B1'].replace(9, numpy.nan)
sub2['S3AQ3C1']=sub2['S3AQ3C1'].replace(99, numpy.nan)

#recoding values for S3AQ3B1 into a new variable, USFREQMO
recode1 = {1: 30, 2: 22, 3: 14, 4: 6, 5: 2.5, 6: 1}
sub2['USFREQMO']= sub2['S3AQ3B1'].map(recode1)

/usr/local/lib/python3.4/dist-packages/ipykernel/__main__.py:1: FutureWarning: convert_objects is deprecated.  Use the data-type specific converters pd.to_datetime, pd.to_timedelta and pd.to_numeric.
  if __name__ == '__main__':
/usr/local/lib/python3.4/dist-packages/ipykernel/__main__.py:2: FutureWarning: convert_objects is deprecated.  Use the data-type specific converters pd.to_datetime, pd.to_timedelta and pd.to_numeric.
  from ipykernel import kernelapp as app
/usr/local/lib/python3.4/dist-packages/ipykernel/__main__.py:3: FutureWarning: convert_objects is deprecated.  Use the data-type specific converters pd.to_datetime, pd.to_timedelta and pd.to_numeric.
  app.launch_new_instance()
/usr/local/lib/python3.4/dist-packages/ipykernel/__main__.py:4: FutureWarning: convert_objects is deprecated.  Use the data-type specific converters pd.to_datetime, pd.to_timedelta and pd.to_numeric.
/usr/local/lib/python3.4/dist-packages/ipykernel/__main__.py:5: FutureWarning: convert_objects is deprecated.  Use the data-type specific converters pd.to_datetime, pd.to_timedelta and pd.to_numeric.

In [43]:

# contingency table of observed counts
ct1=pandas.crosstab(sub2['TABP12MDX'], sub2['USFREQMO'])
print (ct1)

# column percentages
colsum=ct1.sum(axis=0)
colpct=ct1/colsum
print(colpct)

# chi-square
print ('chi-square value, p value, expected counts')
cs1= scipy.stats.chi2_contingency(ct1)
print (cs1)

USFREQMO  1.0   2.5   6.0   14.0  22.0  30.0
TAB12MDX                                    
0           64    53    69    59    41   521
1            7    12    19    32    27   799
USFREQMO      1.0       2.5       6.0       14.0      22.0      30.0
TAB12MDX                                                            
0         0.901408  0.815385  0.784091  0.648352  0.602941  0.394697
1         0.098592  0.184615  0.215909  0.351648  0.397059  0.605303
chi-square value, p value, expected counts
(165.27320708055845, 7.4363642083905987e-34, 5, array([[  33.64474457,   30.80152672,   41.70052848,   43.1221374 ,
          32.22313564,  625.50792719],
       [  37.35525543,   34.19847328,   46.29947152,   47.8778626 ,
          35.77686436,  694.49207281]]))

In [45]:

# set variable types 
sub2["USFREQMO"] = sub2["USFREQMO"].astype('category')
# new code for setting variables to numeric:
sub2['TAB12MDX'] = pandas.to_numeric(sub2['TABP12MDX'], errors='coerce')

# old code for setting variables to numeric:
#sub2['TAB12MDX'] = sub2['TAB12MDX'].convert_objects(convert_numeric=True)

# graph percent with nicotine dependence within each smoking frequency group 
seaborn.factorplot(x="USFREQMO", y="TAB12MDX", data=sub2, kind="bar", ci=None)
plt.xlabel('Days smoked per month')
plt.ylabel('Proportion Nicotine Dependent')

recode2 = {1: 1, 2.5: 2.5}
sub2['COMP1v2']= sub2['USFREQMO'].map(recode2)

# contingency table of observed counts
ct2=pandas.crosstab(sub2['TAB12MDX'], sub2['COMP1v2'])
print (ct2)

# column percentages
colsum=ct2.sum(axis=0)
colpct=ct2/colsum
print(colpct)

print ('chi-square value, p value, expected counts')
cs2= scipy.stats.chi2_contingency(ct2)
print (cs2)

recode3 = {1: 1, 6: 6}
sub2['COMP1v6']= sub2['USFREQMO'].map(recode3)

# contingency table of observed counts
ct3=pandas.crosstab(sub2['TAB12MDX'], sub2['COMP1v6'])
print (ct3)

# column percentages
colsum=ct3.sum(axis=0)
colpct=ct3/colsum
print(colpct)

print ('chi-square value, p value, expected counts')
cs3= scipy.stats.chi2_contingency(ct3)
print (cs3)

recode4 = {1: 1, 14: 14}
sub2['COMP1v14']= sub2['USFREQMO'].map(recode4)

# contingency table of observed counts
ct4=pandas.crosstab(sub2['TAB12MDX'], sub2['COMP1v14'])
print (ct4)

# column percentages
colsum=ct4.sum(axis=0)
colpct=ct4/colsum
print(colpct)

print ('chi-square value, p value, expected counts')
cs4= scipy.stats.chi2_contingency(ct4)
print (cs4)

recode5 = {1: 1, 22: 22}
sub2['COMP1v22']= sub2['USFREQMO'].map(recode5)

# contingency table of observed counts
ct5=pandas.crosstab(sub2['TAB12MDX'], sub2['COMP1v22'])
print (ct5)

# column percentages
colsum=ct5.sum(axis=0)
colpct=ct5/colsum
print(colpct)

print ('chi-square value, p value, expected counts')
cs5= scipy.stats.chi2_contingency(ct5)
print (cs5)

recode6 = {1: 1, 30: 30}
sub2['COMP1v30']= sub2['USFREQMO'].map(recode6)

# contingency table of observed counts
ct6=pandas.crosstab(sub2['TAB12MDX'], sub2['COMP1v30'])
print (ct6)

# column percentages
colsum=ct6.sum(axis=0)
colpct=ct6/colsum
print(colpct)

print ('chi-square value, p value, expected counts')
cs6= scipy.stats.chi2_contingency(ct6)
print (cs6)

recode7 = {2.5: 2.5, 6: 6}
sub2['COMP2v6']= sub2['USFREQMO'].map(recode7)

# contingency table of observed counts
ct7=pandas.crosstab(sub2['TAB12MDX'], sub2['COMP2v6'])
print (ct7)

# column percentages
colsum=ct7.sum(axis=0)
colpct=ct7/colsum
print(colpct)

print ('chi-square value, p value, expected counts')
cs7=scipy.stats.chi2_contingency(ct7)
print (cs7)

COMP1v2   1.0  2.5
TAB12MDX          
0          55   49
1          16   16
COMP1v2        1.0       2.5
TAB12MDX                    
0         0.774648  0.753846
1         0.225352  0.246154
chi-square value, p value, expected counts
(0.0069422451870988916, 0.93359698915826184, 1, array([[ 54.29411765,  49.70588235],
       [ 16.70588235,  15.29411765]]))
COMP1v6    1   6
TAB12MDX        
0         55  71
1         16  17
COMP1v6          1         6
TAB12MDX                    
0         0.774648  0.806818
1         0.225352  0.193182
chi-square value, p value, expected counts
(0.090349171475131884, 0.76373372487916535, 1, array([[ 56.26415094,  69.73584906],
       [ 14.73584906,  18.26415094]]))
COMP1v14  1   14
TAB12MDX        
0         55  64
1         16  27
COMP1v14        1         14
TAB12MDX                    
0         0.774648  0.703297
1         0.225352  0.296703
chi-square value, p value, expected counts
(0.70756615276814527, 0.40025297614002775, 1, array([[ 52.15432099,  66.84567901],
       [ 18.84567901,  24.15432099]]))
COMP1v22  1   22
TAB12MDX        
0         55  41
1         16  27
COMP1v22        1         22
TAB12MDX                    
0         0.774648  0.602941
1         0.225352  0.397059
chi-square value, p value, expected counts
(4.0231238896127559, 0.044880510850748734, 1, array([[ 49.03597122,  46.96402878],
       [ 21.96402878,  21.03597122]]))
COMP1v30  1    30
TAB12MDX         
0         55  624
1         16  696
COMP1v30        1         30
TAB12MDX                    
0         0.774648  0.472727
1         0.225352  0.527273
chi-square value, p value, expected counts
(23.387240486657827, 1.3245419870769008e-06, 1, array([[  34.65780014,  644.34219986],
       [  36.34219986,  675.65780014]]))
COMP2v6   2.5  6.0
TAB12MDX          
0          49   71
1          16   17
COMP2v6        2.5       6.0
TAB12MDX                    
0         0.753846  0.806818
1         0.246154  0.193182
chi-square value, p value, expected counts
(0.34652664097266372, 0.55608593058763378, 1, array([[ 50.98039216,  69.01960784],
       [ 14.01960784,  18.98039216]]))

A good political score explains variation in alcohol consumption?

2016-02-14T14:09:00.002+01:00

Does political score explains variation in alcohol consumption?

My model with data from Gapminder is:

model2 = smf.ols(formula='alcconsumption ~ C(polityscore)', data=sub).fit()

print (model2.summary())

At first look the result looks like there polityscore can explain variations in alcohol consumption.

Output:

OLS Regression Results
==============================================================================
Dep. Variable: alcconsumption R-squared: 0.326
Model: OLS Adj. R-squared: 0.227
Method: Least Squares F-statistic: 3.309
Date: Sun, 14 Feb 2016 Prob (F-statistic): 1.79e-05
Time: 13:54:24 Log-Likelihood: -451.31
No. Observations: 158 AIC: 944.6
Df Residuals: 137 BIC: 1009.
Df Model: 20
Covariance Type: nonrobust
==========================================================================================
coef std err t P>|t| [95.0% Conf. Int.]
------------------------------------------------------------------------------------------
Intercept 0.8150 3.197 0.255 0.799 -5.507 7.137
C(polityscore)[T.-9.0] 3.7383 4.127 0.906 0.367 -4.423 11.899
C(polityscore)[T.-8.0] -0.0950 4.521 -0.021 0.983 -9.035 8.845
C(polityscore)[T.-7.0] 4.3783 3.453 1.268 0.207 -2.450 11.206
C(polityscore)[T.-6.0] 3.4917 4.127 0.846 0.399 -4.669 11.653
C(polityscore)[T.-5.0] 4.0350 4.521 0.892 0.374 -4.905 12.975
C(polityscore)[T.-4.0] 2.6267 3.691 0.712 0.478 -4.673 9.926
C(polityscore)[T.-3.0] 2.9000 3.691 0.786 0.433 -4.400 10.200
C(polityscore)[T.-2.0] 1.5470 3.783 0.409 0.683 -5.933 9.027
C(polityscore)[T.-1.0] 6.2400 3.915 1.594 0.113 -1.502 13.982
C(polityscore)[T.0.0] 1.8700 3.691 0.507 0.613 -5.430 9.170
C(polityscore)[T.1.0] 3.7783 4.127 0.915 0.362 -4.383 11.939
C(polityscore)[T.2.0] 1.6083 4.127 0.390 0.697 -6.553 9.769
C(polityscore)[T.3.0] 4.1850 4.521 0.926 0.356 -4.755 13.125
C(polityscore)[T.4.0] 9.1025 3.915 2.325 0.022 1.360 16.845
C(polityscore)[T.5.0] 4.1693 3.625 1.150 0.252 -2.999 11.337
C(polityscore)[T.6.0] 4.3460 3.502 1.241 0.217 -2.579 11.271
C(polityscore)[T.7.0] 3.8642 3.434 1.125 0.262 -2.926 10.655
C(polityscore)[T.8.0] 8.2171 3.361 2.445 0.016 1.571 14.863
C(polityscore)[T.9.0] 8.1736 3.418 2.392 0.018 1.416 14.932
C(polityscore)[T.10.0] 9.6331 3.295 2.923 0.004 3.117 16.149
==============================================================================
Omnibus: 17.250 Durbin-Watson: 1.909
Prob(Omnibus): 0.000 Jarque-Bera (JB): 19.895
Skew: 0.750 Prob(JB): 4.79e-05
Kurtosis: 3.880 Cond. No. 43.1
==============================================================================

Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.

Then we have a look at the relationship between mean polityscore and alcconsumption

m2= sub.groupby('polityscore').mean()

means for alcconsumption by polityscore
             alcconsumption
polityscore                
-10.00                 0.82
-9.00                  4.55
-8.00                  0.72
-7.00                  5.19
-6.00                  4.31
-5.00                  4.85
-4.00                  3.44
-3.00                  3.71
-2.00                  2.36
-1.00                  7.05
0.00                   2.69
1.00                   4.59
2.00                   2.42
3.00                   5.00
4.00                   9.92
5.00                   4.98
6.00                   5.16
7.00                   4.68
8.00                   9.03
9.00                   8.99
10.00                 10.45

The null hypoteses is that the alcohol variations is explained by politics. A tukey test, that accounts for what levels of politicscore are affecting the output shows that there is seldom a level that rejects the null hypotesis. We should be doubting that there is a influence.

import statsmodels.stats.multicomp as multi 
mc1 = multi.MultiComparison(sub['alcconsumption'], sub['polityscore'])
res1 = mc1.tukeyhsd()
print(res1.summary())

Multiple Comparison of Means - Tukey HSD,FWER=0.05
==============================================
group1 group2 meandiff  lower    upper  reject
----------------------------------------------
-10.0   -9.0   3.7383  -11.2896 18.7663 False 
-10.0   -8.0   -0.095  -16.5573 16.3673 False 
-10.0   -7.0   4.3783  -8.1949  16.9516 False 
-10.0   -6.0   3.4917  -11.5363 18.5196 False 
-10.0   -5.0   4.035   -12.4273 20.4973 False 
-10.0   -4.0   2.6267  -10.8147 16.0681 False 
-10.0   -3.0    2.9    -10.5414 16.3414 False 
-10.0   -2.0   1.547   -12.2263 15.3203 False 
-10.0   -1.0    6.24   -8.0168  20.4968 False 
-10.0   0.0     1.87   -11.5714 15.3114 False 
-10.0   1.0    3.7783  -11.2496 18.8063 False 
-10.0   2.0    1.6083  -13.4196 16.6363 False 
-10.0   3.0    4.185   -12.2773 20.6473 False 
-10.0   4.0    9.1025  -5.1543  23.3593 False 
-10.0   5.0    4.1693  -9.0299  17.3685 False 
-10.0   6.0    4.346   -8.4056  17.0976 False 
-10.0   7.0    3.8642  -8.6398  16.3682 False 
-10.0   8.0    8.2171  -4.0208   20.455 False 
-10.0   9.0    8.1736  -4.2707  20.6179 False 
-10.0   10.0   9.6331  -2.3657   21.632 False 
 -9.0   -8.0  -3.8333  -18.8613 11.1946 False 
 -9.0   -7.0    0.64   -9.9864  11.2664 False 
 -9.0   -6.0  -0.2467  -13.6881 13.1947 False 
 -9.0   -5.0   0.2967  -14.7313 15.3246 False 
 -9.0   -4.0  -1.1117  -12.7523 10.5289 False 
 -9.0   -3.0  -0.8383  -12.4789 10.8023 False 
 -9.0   -2.0  -2.1913  -14.2137  9.831  False 
 -9.0   -1.0   2.5017  -10.0716 15.0749 False 
 -9.0   0.0   -1.8683  -13.5089  9.7723 False 
 -9.0   1.0     0.04   -13.4014 13.4814 False 
 -9.0   2.0    -2.13   -15.5714 11.3114 False 
 -9.0   3.0    0.4467  -14.5813 15.4746 False 
 -9.0   4.0    5.3642  -7.2091  17.9374 False 
 -9.0   5.0    0.431   -10.9291  11.791 False 
 -9.0   6.0    0.6077  -10.2291 11.4445 False 
 -9.0   7.0    0.1259  -10.4184 10.6702 False 
 -9.0   8.0    4.4788  -5.7486  14.7061 False 
 -9.0   9.0    4.4352  -6.0382  14.9087 False 
 -9.0   10.0   5.8948  -4.0453  15.8348 False 
 -8.0   -7.0   4.4733  -8.0999  17.0466 False 
 -8.0   -6.0   3.5867  -11.4413 18.6146 False 
 -8.0   -5.0    4.13   -12.3323 20.5923 False 
 -8.0   -4.0   2.7217  -10.7197 16.1631 False 
 -8.0   -3.0   2.995   -10.4464 16.4364 False 
 -8.0   -2.0   1.642   -12.1313 15.4153 False 
 -8.0   -1.0   6.335   -7.9218  20.5918 False 
 -8.0   0.0    1.965   -11.4764 15.4064 False 
 -8.0   1.0    3.8733  -11.1546 18.9013 False 
 -8.0   2.0    1.7033  -13.3246 16.7313 False 
 -8.0   3.0     4.28   -12.1823 20.7423 False 
 -8.0   4.0    9.1975  -5.0593  23.4543 False 
 -8.0   5.0    4.2643  -8.9349  17.4635 False 
 -8.0   6.0    4.441   -8.3106  17.1926 False 
 -8.0   7.0    3.9592  -8.5448  16.4632 False 
 -8.0   8.0    8.3121  -3.9258   20.55  False 
 -8.0   9.0    8.2686  -4.1757  20.7129 False 
 -8.0   10.0   9.7281  -2.2707   21.727 False 
 -7.0   -6.0  -0.8867  -11.513   9.7397 False 
 -7.0   -5.0  -0.3433  -12.9166 12.2299 False 
 -7.0   -4.0  -1.7517  -9.9828   6.4795 False 
 -7.0   -3.0  -1.4783  -9.7095   6.7528 False 
 -7.0   -2.0  -2.8313  -11.5941  5.9314 False 
 -7.0   -1.0   1.8617  -7.6428  11.3662 False 
 -7.0   0.0   -2.5083  -10.7395  5.7228 False 
 -7.0   1.0     -0.6   -11.2264 10.0264 False 
 -7.0   2.0    -2.77   -13.3964  7.8564 False 
 -7.0   3.0   -0.1933  -12.7666 12.3799 False 
 -7.0   4.0    4.7242  -4.7803  14.2287 False 
 -7.0   5.0    -0.209  -8.0384   7.6203 False 
 -7.0   6.0   -0.0323  -7.0811   7.0164 False 
 -7.0   7.0   -0.5141  -7.1043   6.0761 False 
 -7.0   8.0    3.8388  -2.2314   9.909  False 
 -7.0   9.0    3.7952   -2.681  10.2715 False 
 -7.0   10.0   5.2548  -0.3177  10.8273 False 
 -6.0   -5.0   0.5433  -14.4846 15.5713 False 
 -6.0   -4.0   -0.865  -12.5056 10.7756 False 
 -6.0   -3.0  -0.5917  -12.2323 11.0489 False 
 -6.0   -2.0  -1.9447  -13.967  10.0777 False 
 -6.0   -1.0   2.7483  -9.8249  15.3216 False 
 -6.0   0.0   -1.6217  -13.2623 10.0189 False 
 -6.0   1.0    0.2867  -13.1547 13.7281 False 
 -6.0   2.0   -1.8833  -15.3247 11.5581 False 
 -6.0   3.0    0.6933  -14.3346 15.7213 False 
 -6.0   4.0    5.6108  -6.9624  18.1841 False 
 -6.0   5.0    0.6776  -10.6824 12.0377 False 
 -6.0   6.0    0.8543  -9.9825  11.6911 False 
 -6.0   7.0    0.3726  -10.1717 10.9169 False 
 -6.0   8.0    4.7254  -5.5019  14.9528 False 
 -6.0   9.0    4.6819  -5.7916  15.1554 False 
 -6.0   10.0   6.1415  -3.7986  16.0815 False 
 -5.0   -4.0  -1.4083  -14.8497 12.0331 False 
 -5.0   -3.0   -1.135  -14.5764 12.3064 False 
 -5.0   -2.0   -2.488  -16.2613 11.2853 False 
 -5.0   -1.0   2.205   -12.0518 16.4618 False 
 -5.0   0.0    -2.165  -15.6064 11.2764 False 
 -5.0   1.0   -0.2567  -15.2846 14.7713 False 
 -5.0   2.0   -2.4267  -17.4546 12.6013 False 
 -5.0   3.0     0.15   -16.3123 16.6123 False 
 -5.0   4.0    5.0675  -9.1893  19.3243 False 
 -5.0   5.0    0.1343  -13.0649 13.3335 False 
 -5.0   6.0    0.311   -12.4406 13.0626 False 
 -5.0   7.0   -0.1708  -12.6748 12.3332 False 
 -5.0   8.0    4.1821  -8.0558   16.42  False 
 -5.0   9.0    4.1386  -8.3057  16.5829 False 
 -5.0   10.0   5.5981  -6.4007   17.597 False 
 -4.0   -3.0   0.2733  -9.2312   9.7778 False 
 -4.0   -2.0  -1.0797  -11.0481  8.8887 False 
 -4.0   -1.0   3.6133   -7.013  14.2397 False 
 -4.0   0.0   -0.7567  -10.2612  8.7478 False 
 -4.0   1.0    1.1517  -10.4889 12.7923 False 
 -4.0   2.0   -1.0183  -12.6589 10.6223 False 
 -4.0   3.0    1.5583  -11.8831 14.9997 False 
 -4.0   4.0    6.4758  -4.1505  17.1022 False 
 -4.0   5.0    1.5426  -7.6162  10.7014 False 
 -4.0   6.0    1.7193  -6.7818  10.2204 False 
 -4.0   7.0    1.2376  -6.8874   9.3625 False 
 -4.0   8.0    5.5904  -2.1187  13.2996 False 
 -4.0   9.0    5.5469  -2.4859  13.5797 False 
 -4.0   10.0   7.0065  -0.3173  14.3302 False 
 -3.0   -2.0   -1.353  -11.3214  8.6154 False 
 -3.0   -1.0    3.34   -7.2864  13.9664 False 
 -3.0   0.0    -1.03   -10.5345  8.4745 False 
 -3.0   1.0    0.8783  -10.7623 12.5189 False 
 -3.0   2.0   -1.2917  -12.9323 10.3489 False 
 -3.0   3.0    1.285   -12.1564 14.7264 False 
 -3.0   4.0    6.2025  -4.4239  16.8289 False 
 -3.0   5.0    1.2693  -7.8895  10.4281 False 
 -3.0   6.0    1.446   -7.0551   9.9471 False 
 -3.0   7.0    0.9642  -7.1607   9.0892 False 
 -3.0   8.0    5.3171  -2.3921  13.0263 False 
 -3.0   9.0    5.2736  -2.7592  13.3063 False 
 -3.0   10.0   6.7331  -0.5906  14.0568 False 
 -2.0   -1.0   4.693   -6.3502  15.7362 False 
 -2.0   0.0    0.323   -9.6454  10.2914 False 
 -2.0   1.0    2.2313   -9.791  14.2537 False 
 -2.0   2.0    0.0613  -11.961  12.0837 False 
 -2.0   3.0    2.638   -11.1353 16.4113 False 
 -2.0   4.0    7.5555  -3.4877  18.5987 False 
 -2.0   5.0    2.6223   -7.017  12.2616 False 
 -2.0   6.0    2.799   -6.2178  11.8158 False 
 -2.0   7.0    2.3172  -6.3458  10.9803 False 
 -2.0   8.0    6.6701  -1.6042  14.9445 False 
 -2.0   9.0    6.6266  -1.9501  15.2032 False 
 -2.0   10.0   8.0861   0.1697  16.0026  True 
 -1.0   0.0    -4.37   -14.9964  6.2564 False 
 -1.0   1.0   -2.4617  -15.0349 10.1116 False 
 -1.0   2.0   -4.6317  -17.2049  7.9416 False 
 -1.0   3.0    -2.055  -16.3118 12.2018 False 
 -1.0   4.0    2.8625  -8.7781  14.5031 False 
 -1.0   5.0   -2.0707  -12.389   8.2476 False 
 -1.0   6.0    -1.894  -11.6332  7.8452 False 
 -1.0   7.0   -2.3758  -11.7884  7.0369 False 
 -1.0   8.0    1.9771  -7.0791  11.0333 False 
 -1.0   9.0    1.9336  -7.3997  11.2668 False 
 -1.0   10.0   3.3931  -5.3373  12.1236 False 
 0.0    1.0    1.9083  -9.7323  13.5489 False 
 0.0    2.0   -0.2617  -11.9023 11.3789 False 
 0.0    3.0    2.315   -11.1264 15.7564 False 
 0.0    4.0    7.2325  -3.3939  17.8589 False 
 0.0    5.0    2.2993  -6.8595  11.4581 False 
 0.0    6.0    2.476   -6.0251  10.9771 False 
 0.0    7.0    1.9942  -6.1307  10.1192 False 
 0.0    8.0    6.3471  -1.3621  14.0563 False 
 0.0    9.0    6.3036  -1.7292  14.3363 False 
 0.0    10.0   7.7631   0.4394  15.0868  True 
 1.0    2.0    -2.17   -15.6114 11.2714 False 
 1.0    3.0    0.4067  -14.6213 15.4346 False 
 1.0    4.0    5.3242  -7.2491  17.8974 False 
 1.0    5.0    0.391   -10.9691  11.751 False 
 1.0    6.0    0.5677  -10.2691 11.4045 False 
 1.0    7.0    0.0859  -10.4584 10.6302 False 
 1.0    8.0    4.4388  -5.7886  14.6661 False 
 1.0    9.0    4.3952  -6.0782  14.8687 False 
 1.0    10.0   5.8548  -4.0853  15.7948 False 
 2.0    3.0    2.5767  -12.4513 17.6046 False 
 2.0    4.0    7.4942  -5.0791  20.0674 False 
 2.0    5.0    2.561   -8.7991   13.921 False 
 2.0    6.0    2.7377  -8.0991  13.5745 False 
 2.0    7.0    2.2559  -8.2884  12.8002 False 
 2.0    8.0    6.6088  -3.6186  16.8361 False 
 2.0    9.0    6.5652  -3.9082  17.0387 False 
 2.0    10.0   8.0248  -1.9153  17.9648 False 
 3.0    4.0    4.9175  -9.3393  19.1743 False 
 3.0    5.0   -0.0157  -13.2149 13.1835 False 
 3.0    6.0    0.161   -12.5906 12.9126 False 
 3.0    7.0   -0.3208  -12.8248 12.1832 False 
 3.0    8.0    4.0321  -8.2058   16.27  False 
 3.0    9.0    3.9886  -8.4557  16.4329 False 
 3.0    10.0   5.4481  -6.5507   17.447 False 
 4.0    5.0   -4.9332  -15.2515  5.3851 False 
 4.0    6.0   -4.7565  -14.4957  4.9827 False 
 4.0    7.0   -5.2383  -14.6509  4.1744 False 
 4.0    8.0   -0.8854  -9.9416   8.1708 False 
 4.0    9.0   -0.9289  -10.2622  8.4043 False 
 4.0    10.0   0.5306  -8.1998   9.2611 False 
 5.0    6.0    0.1767   -7.936   8.2894 False 
 5.0    7.0   -0.3051  -8.0227   7.4126 False 
 5.0    8.0    4.0478  -3.2308  11.3265 False 
 5.0    9.0    4.0043  -3.6163  11.6248 False 
 5.0    10.0   5.4638  -1.4052  12.3329 False 
 6.0    7.0   -0.4818  -7.4062   6.4426 False 
 6.0    8.0    3.8711  -2.5604  10.3026 False 
 6.0    9.0    3.8276  -2.9885  10.6436 False 
 6.0    10.0   5.2871  -0.6769  11.2512 False 
 7.0    8.0    4.3529  -1.5725  10.2783 False 
 7.0    9.0    4.3093  -2.0313   10.65  False 
 7.0    10.0   5.7689   0.3545  11.1833  True 
 8.0    9.0   -0.0435  -5.8419   5.7548 False 
 8.0    10.0   1.416   -3.3518   6.1839 False 
 9.0    10.0   1.4596  -3.8155   6.7346 False 
----------------------------------------------

Males are more likely to have a major life depression or panic disorders

2016-02-08T10:11:00.000+01:00

1. Summary

1) what you found, making sure you discuss the results for the associations between all of your explanatory variables and your response variable. Make sure to include statistical results (odds ratios, p-values, and 95% confidence intervals for the odds ratios) in your summary.

After adjusting for potential confounding factors (Nicotine dependency, social phobia), the odds of being a male were more than two times higher for participants with major depression than for participants without major depression (OR=1.69, 95% CI = 1.31-2.20, p<.0001).
Panic was also significantly associated with being male, such that participants with panic problems were significantly more likely to be males (OR= 1.51, 95% CI=1.1-2.09, p=.011).
Luckily the Model

2) Report whether or not your results supported your hypothesis for the association between your primary explanatory variable and your response variable.

My primary hypothesis that Males are more likely to have major life depressions is supported by the model
3) Discuss whether or not there was evidence of confounding for the association between your primary explanatory and the response variable (Hint: adding additional explanatory variables to your model one at a time will make it easier to identify which of the variables are confounding variables).

There was no evidence of confounders

2. the output from your logistic regression model.

                           Logit Regression Results                           
==============================================================================
Dep. Variable:                   SEXO   No. Observations:                 1320
Model:                          Logit   Df Residuals:                     1317
Method:                           MLE   Df Model:                            2
Date:                Mon, 08 Feb 2016   Pseudo R-squ.:                 0.01697
Time:                        09:49:06   Log-Likelihood:                -899.14
converged:                       True   LL-Null:                       -914.66
                                        LLR p-value:                 1.817e-07
================================================================================
                   coef    std err          z      P>|z|      [95.0% Conf. Int.]
--------------------------------------------------------------------------------
Intercept       -0.1599      0.066     -2.413      0.016        -0.290    -0.030
PANIC            0.4152      0.164      2.529      0.011         0.093     0.737
MAJORDEPLIFE     0.5272      0.132      3.980      0.000         0.268     0.787
================================================================================
Odds Ratios
              Lower CI  Upper CI   OR
Intercept         0.75      0.97 0.85
PANIC             1.10      2.09 1.51
MAJORDEPLIFE      1.31      2.20 1.69

2016-02-07T14:37:00.001+01:00

Countries with high Internet use have a high Alcohol Consumption

2016-02-07T14:37:00.000+01:00

The Study uses the Gapminder Dataset shows a correlation between Alcohol and Internet. After controlling for urbanization rate, the a doubling in internet usage accounts for beta=11.5% of Alcohol consumption (p< 0,005). A second order Regression alcconsumption ~ urbanrate_c + I(urbanrate_c**2) + internetuserate_c shows a that 34% of the variability can be explained by the model, R^2=0,343, F statistic= 2.56e-11.

Intercept sayst that Alcohol consumption with average internet use and urbanisation rate (variables are centered) is 8.63 litres (estimated average alcohol consumption, adult (15+) per capita consumption in litres pure alcohol)

urban rate (centered) has a small negative effect ~ -6% on alcohol consumption, but increases the explanatory power about 10% (compare the R^2 of the models with/without control for urbanization)

I(urbanrate_c ** 2) is significant with littel effect (- 0,31%)

                            OLS Regression Results                            
==============================================================================
Dep. Variable:         alcconsumption   R-squared:                       0.343
Model:                            OLS   Adj. R-squared:                  0.327
Method:                 Least Squares   F-statistic:                     21.59
Date:                Sun, 07 Feb 2016   Prob (F-statistic):           2.56e-11
Time:                        12:40:28   Log-Likelihood:                -365.20
No. Observations:                 128   AIC:                             738.4
Df Residuals:                     124   BIC:                             749.8
Df Model:                           3                                         
Covariance Type:            nonrobust                                         
=======================================================================================
                          coef    std err          t      P>|t|      [95.0% Conf. Int.]
---------------------------------------------------------------------------------------
Intercept               8.6291      0.501     17.230      0.000         7.638     9.620
urbanrate_c            -0.0596      0.024     -2.466      0.015        -0.107    -0.012
I(urbanrate_c ** 2)    -0.0031      0.001     -3.980      0.000        -0.005    -0.002
internetuserate_c       0.1149      0.018      6.558      0.000         0.080     0.150
==============================================================================
Omnibus:                        5.219   Durbin-Watson:                   1.972
Prob(Omnibus):                  0.074   Jarque-Bera (JB):                5.213
Skew:                           0.316   Prob(JB):                       0.0738
Kurtosis:                       3.760   Cond. No.                         891.
==============================================================================

The result strongly suggests that there is a correlation between alcohol and internet use. It is both highly statistically significant (p>0,0005) and has an effect of beta = 11.5%

There is no evidence of confounding since the internet use rate did not lose statistical significance after the addition of urbanization rate as an explanatory variable.

q-q plot: Residuals are pretty normaly distributed with two outliers on the upper quantiles. This indicates that there is very little systematic error, so probably the model can not be optimized further by adding other variables.

standardized residuals for all observations: We see that most of the residuals are in between two standard deviations. There are few major outliers, so the model captures well the variability in the dataset.

Also, these plots show that the data is pretty normal and explain the response variable well. "Urbanrate" has normally distributed residuals and a well fitting regression line with small negative effect on the response.

leverage plot: The data that has the most influence on the regression outcome has little residuals, which supports our claim that the model fits well the data and has predictive power.

Coursera Regression Course - Assignment 2

2016-01-29T17:25:00.000+01:00

A major life depression is associated with
a higher daily consumption of alcohol

All data from the NESARC Study

Categorical Explanatory Variable

We have a categorical variable encoded with No=0, implying that subjects with no major depression in life are at the baseline 0.

3549-3549
MAJORDEPLIFE
MAJOR DEPRESSION - LIFETIME (NON-HIERARCHICAL)
----------------------------------------------
35254 0. No
7839 1. Yes
----------------------------------------------

Quantitative Response Variable

The amount of ethanol (alcohol) consumed averagely in the past year is our response variable.

3675-3682
ETOTLCA2
AVERAGE DAILY VOLUME OF ETHANOL CONSUMED IN PAST YEAR, FROM ALL TYPES OF ALCOHOLIC BEVERAGES COMBINED
(NOTE: Users may wish to exclude outliers)
--------------------------------------------------
0.0003 - 219.9555 Ounces of ethanol
Blank Unknown
--------------------------------------------------

We see that there is a large number of abstinent people in the dataset. We still don't know if they are associated with depression or not. It would be interesting to see if they are more the one or the other.

Linear Regression Test

                            OLS Regression Results                            
==============================================================================
Dep. Variable:               ETOTLCA2   R-squared:                       0.001
Model:                            OLS   Adj. R-squared:                  0.001
Method:                 Least Squares   F-statistic:                     18.86
Date:                Fri, 29 Jan 2016   Prob (F-statistic):           1.42e-05
Time:                        15:57:58   Log-Likelihood:                -59022.
No. Observations:               26655   AIC:                         1.180e+05
Df Residuals:                   26653   BIC:                         1.181e+05
Df Model:                           1                                         
Covariance Type:            nonrobust                                         
================================================================================
                   coef    std err          t      P>|t|      [95.0% Conf. Int.]
--------------------------------------------------------------------------------
Intercept        0.5366      0.015     35.377      0.000         0.507     0.566
MAJORDEPLIFE     0.1474      0.034      4.342      0.000         0.081     0.214
==============================================================================
Omnibus:                    84158.156   Durbin-Watson:                   2.010
Prob(Omnibus):                  0.000   Jarque-Bera (JB):      20279387185.774
Skew:                          50.216   Prob(JB):                         0.00
Kurtosis:                    4274.926   Cond. No.                         2.62
==============================================================================

The R-squared number is rather low, sugesting that only 0,1% of the "Etanolic Intake" can be explained by a "major depression in life". At least the F-statistic is so low that we can assume the small effect is statistically significant.

Now I would not go further but lets assume that R-squared is higher, in that case we believe that a major depression in a lifetime can predict Beta=14.74% increase in alcohol intake, with a high statistical significance of P>0.000

Graph: a major life depression could increase average ethanol intake from 54% to 68%

Program Code and Output

In [1]:

%matplotlib inline

import numpy
import pandas
import statsmodels.api as sm
import seaborn
import matplotlib.pyplot as plt
import statsmodels.formula.api as smf

In [2]:

# bug fix for display formats to avoid run time errors
pandas.set_option('display.float_format', lambda x:'%.2f'%x)

In [3]:

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

In [4]:

#setting variables you will be working with to numeric
data['p'] = pandas.to_numeric(data['MAJORDEPLIFE'], errors='coerce')

In [5]:

data['ETOTLCA2'] = pandas.to_numeric(data['ETOTLCA2'], errors='coerce')

In [6]:

subGroupEthanol = list(filter(lambda x: x > 10.0, data['ETOTLCA2'].dropna()))
seaborn.distplot(subGroupEthanol)

Out[6]:

In [7]:

reg1 = smf.ols('ETOTLCA2 ~ MAJORDEPLIFE', data=data).fit()
print (reg1.summary())

                            OLS Regression Results                            
==============================================================================
Dep. Variable:               ETOTLCA2   R-squared:                       0.001
Model:                            OLS   Adj. R-squared:                  0.001
Method:                 Least Squares   F-statistic:                     18.86
Date:                Fri, 29 Jan 2016   Prob (F-statistic):           1.42e-05
Time:                        17:11:29   Log-Likelihood:                -59022.
No. Observations:               26655   AIC:                         1.180e+05
Df Residuals:                   26653   BIC:                         1.181e+05
Df Model:                           1                                         
Covariance Type:            nonrobust                                         
================================================================================
                   coef    std err          t      P>|t|      [95.0% Conf. Int.]
--------------------------------------------------------------------------------
Intercept        0.5366      0.015     35.377      0.000         0.507     0.566
MAJORDEPLIFE     0.1474      0.034      4.342      0.000         0.081     0.214
==============================================================================
Omnibus:                    84158.156   Durbin-Watson:                   2.010
Prob(Omnibus):                  0.000   Jarque-Bera (JB):      20279387185.774
Skew:                          50.216   Prob(JB):                         0.00
Kurtosis:                    4274.926   Cond. No.                         2.62
==============================================================================

Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.

In [8]:

# listwise deletion for calculating means for regression model observations
sub1 = data[['ETOTLCA2', 'MAJORDEPLIFE']].dropna()

# group means & sd
print ("Mean")
ds1 = sub1.groupby('MAJORDEPLIFE').mean()
print (ds1)
print ("Standard deviation")
ds2 = sub1.groupby('MAJORDEPLIFE').std()
print (ds2)

# bivariate bar graph
seaborn.factorplot(x="MAJORDEPLIFE", y="ETOTLCA2", data=sub1, kind="bar", ci=None)
plt.xlabel('Major Life Depression')
plt.ylabel('Mean Number AVERAGE DAILY VOLUME OF ETHANOL CONSUMED IN PAST YEAR')

Mean
              ETOTLCA2
MAJORDEPLIFE          
0                 0.54
1                 0.68
Standard deviation
              ETOTLCA2
MAJORDEPLIFE          
0                 1.44
1                 4.03

Out[8]:

In [ ]:

2016-01-20T20:15:00.001+01:00

First Assignment Coursera - Regression Models in Practice

I have chosen the National Epidemiologic Survey on Alcohol and Related Conditions (NESARC)

http://pubs.niaaa.nih.gov/publications/AA70/AA70.htm

1. Describe the sample

a) Describe the study population (who or what was studied).

NESARC contains an extensive battery of questions about present and past alcohol consumption, AUDs, and the use of alcohol treatment services. NESARC also included similar questions related to tobacco and illicit drug use (including nicotine dependence and drug use disorders) as well as questions designed to determine a wide variety of psychiatric disorders such as major depression, anxiety disorders, and personality disorders.

Furthermore, NESARC contained questions that operationalized the criteria set forth in the American Psychiatric Association’s Diagnostic and Statistical Manual of Mental Disorders, Fourth Edition (DSM–IV) for the following psychiatric disorders:

Five mood disorders (major depressive disorder, bipolar I and bipolar II disorders, dysthymia, and hypomania)
Four anxiety disorders (panic with and without agoraphobia, social phobia, specific phobia, and generalized anxiety)
Seven personality disorders (avoidant, dependent, obsessive–compulsive, paranoid, schizoid, histrionic, and antisocial disorders).

b) Report the level of analysis studied (individual, group, or aggregate).

Individuals were randomly selected from a systematic sample of group quarters in these PSUs.The Census 2000 Group Quarters Inventory formed the sampling frame for the group quarters portion of the NESARC sample.

c) Report the number of observations in the data set.

43093 Total, 8209 1. Northeast 8991 2. Midwest 16156 3. South 9737 4. West

d) Describe your data analytic sample (the sample you are using for your analyses).

397-398 S2AQ8A HOW OFTEN DRANK ANY ALCOHOL IN LAST 12 MONTHS

---------------------------------------------

1865 1. Every day

1210 2. Nearly every day

2619 3. 3 to 4 times a week

2914 4. 2 times a week

3261 5. Once a week

3557 6. 2 to 3 times a month

2663 7. Once a month 1805 8. 7 to 11 times in the last year

3210 9. 3 to 6 times in the last year

3637 10. 1 or 2 times in the last year

205 99. Unknown

16147 BL. NA, former drinker or lifetime abstainer

---------------------------------------------------------------------------------------

2. Describe the data collection procedure

a) Report the study design that generated that data (for example: data reporting, surveys, observation, experiment).

Surveys, Data were collected in face-to-face, computer-assisted personal interviews conducted in respondents’ homes. The NESARC response rate was 81 percent.

b) Describe the original purpose of the data collection.

The major purposes of the Wave 1 and Wave 2 NESARC are to:

Determine the prevalence, incidence, stability, and recurrence of AUDs and their associated disabilities in the general U.S. population.
Estimate the magnitude of health disparities in AUDs and their associated disabilities among population subgroups defined by gender, race/ethnicity, disability, sexual orientation, age, and socioeconomic status.
Estimate the size, characteristics, and changing nature of populations of special concern, including alcohol abusers and other people in the general population who are impaired or affected by the use of alcohol (e.g., those engaging in binge drinking or impaired driving).
Estimate changes in AUDs and their associated disabilities over time, and identify factors associated with the natural history of AUDs.
Determine the number of people receiving alcohol treatment through various treatment programs and services, including those not otherwise represented in surveys of treatment facilities; measure the unmet need for alcohol treatment services; and identify barriers to seeking treatment.
Determine the associations between AUDs and their major physical and mental disabilities, differentiating drug-induced disorders from those reflecting true, independent mental conditions.
Determine the boundaries between safe and hazardous drinking levels and patterns for various types of AUDs and their associated medical, social, and psychological sequelae.

c) Describe how the data were collected.

NESARC participants are variable in lifestyle and ages.

To ensure that minority and special populations were well represented in the sample, NESARC oversampled Blacks, Hispanics, and young adults ages 18–24. As a result, the survey produced enough minority respondents to answer questions of race/ethnic disparities in comorbidity and access to health care services

d) Report when the data were collected.

In 2001—2002, the National Institute on Alcohol Abuse and Alcoholism (NIAAA) conducted the first wave of the National Epidemiologic Survey on Alcohol and Related Conditions (NESARC),The second wave of the survey took place from 2004 to 2005.

e) Report where the data were collected.

They represented all regions of the United States. In addition to sampling people living in traditional households, NESARC investigators questioned military personnel living off base and people living in a variety of group accommodations such as boarding or rooming houses and college quarters. By including these different types of housing, the investigators were able to obtain data on people not typically captured by household surveys.

3. Measures section describing your variables and how you managed them to address your own research question

a) Describe what your explanatory and response variables measured.

explanatory variable:

397-398 S2AQ8A HOW OFTEN DRANK ANY ALCOHOL IN LAST 12 MONTHS

---------------------------------------------

1865 1. Every day

1210 2. Nearly every day

2619 3. 3 to 4 times a week

2914 4. 2 times a week

3261 5. Once a week

3557 6. 2 to 3 times a month

2663 7. Once a month 1805 8. 7 to 11 times in the last year

3210 9. 3 to 6 times in the last year

3637 10. 1 or 2 times in the last year

205 99. Unknown

16147 BL. NA, former drinker or lifetime abstainer

---------------------------------------------------------------------------------------

b) Describe the response scales for your explanatory and response variables.

c) Describe how you managed your explanatory and response variables.

response variable

SECTION 9: GENERALIZED ANXIETY (GENERAL ANXIETY) -
3020-3020 S9Q1A EVER HAD 6+ MONTH PERIOD FELT TENSE/NERVOUS/WORRIED MOST OF TIME
3128 1. Yes 38225 2. No 1740 9. Unknown

3021-3021 S9Q1B EVER HAD 6+ MONTH PERIOD FELT VERY TENSE/NERVOUS/WORRIED MOST OF TIME ABOUT EVERYDAY PROBLEMS
358 1. Yes 37867 2. No 1740 9. Unknown 3128 BL. NA, had 6+ month period feeling tense/nervous/worried most of the time

3022-3022 S9Q31 IN WORST PERIOD, EVER WORRY A LOT ABOUT THINGS YOU USUALLY DIDN'T WORRY ABOUT
2336 1. Yes 1115 2. No 35 9. Unknown 39607 BL. NA, never or unkn. if ever had 6+ month pd. being tense/nervous/worried

3023-3023 S9Q32 IN WORST PERIOD, EVER WORRY ABOUT MORE THAN ONE THING
2644 1. Yes 808 2. No 34 9. Unknown 39607 BL. NA, never or unkn. if ever had 6+ month pd. being tense/nervous/worried

3024-3024 S9Q33 IN WORST PERIOD, EVER FIND IT DIFFICULT TO STOP BEING TENSE/NERVOUS/WORRIED
2715 1. Yes 742 2. No 29 9. Unknown 39607 BL. NA, never or unkn. if ever had 6+ month pd. being tense/nervous/worried

3025-3025 S9Q34 IN WORST PERIOD, EVER WORRY ABOUT THINGS THAT WERE VERY UNLIKELY TO HAPPEN
1418 1. Yes 2019 2. No 49 9. Unknown 39607 BL. NA, never or unkn. if ever had 6+ month pd. being tense/nervous/worried

c) Describe how you managed your explanatory and response variables.

confounders:

I will control for different confounders like:

Page 302: SECTION 4A: MAJOR DEPRESSION (LOW MOOD I)
Page 316: SECTION 4B: FAMILY HISTORY (III) OF MAJOR DEPRESSION
Page 319: SECTION 4C: DYSTHYMIA (LOW MOOD II)
Page 329: SECTION 5: MANIA OR HYPOMANIA (HIGH MOOD)
Page 341: SECTION 6: PANIC DISORDERS AND AGORAPHOBIA (ANXIETY)
Page 355: SECTION 7: SOCIAL PHOBIA (SOCIAL SITUATIONS)
Page 371: SECTION 8: SPECIFIC PHOBIA (SPECIFIC SITUATIONS)
Page 386: SECTION 9: GENERALIZED ANXIETY (GENERAL ANXIETY)
Page 403: SECTION 10: PERSONALITY DISORDERS (USUAL FEELINGS/ACTIONS)
Page 419: SECTION 11A: ANTISOCIAL PERSONALITY DISORDER (BEHAVIOR)

Where do the richest Bakers/ Butchers live?

2014-07-07T20:21:00.000+02:00

Which Travel Guide for Tyrol? To Germans does not matter

2013-09-01T13:49:00.003+02:00

My endevor to build an tourism indicator with online booksales led to a first interesting result.
In the category of worldwide travelguides, you see that the rank number (more is worse) increases comparably fast.
Whereas within the Tyrol category there are alot of books doing comparably well. The rank number increases slow.

World travel guides

Italy travel guides

South Tyrol travel guides

Tyrol travel guides

See you at G+

2012-01-09T21:42:00.001+01:00

As long as I have not a longer article to publish you find daily rants and ramblings on certainty, uncertainty and ignorance under https://plus.google.com/u/0/103983654382933904290/posts

Easy way to stop smoking, this time it works!

2011-08-27T14:54:00.014+02:00

It was Supposed to be Easy
Like millions I read the famous book by Allen Carr "Easy way to stop smoking". Like millions, I tought: he has a point. Many quited, but I went on. Now I'v found a secret weapon of mass responsibility, where I can change my behavior for good. The weapon is called Commitment Contract.

UPDATE: it's now over a month since I am without smoking. It was tough, but you just stick with it.

Why Quit?
The consequences of my behavior are well studied:

Men born in 1900-1930 who smoked only cigarettes and continued smoking died on average about 10 years younger than lifelong non-smokers. Cessation at age 60, 50, 40, or 30 years gained, respectively, about 3, 6, 9, or 10 years of life expectancy.

[ study prolonged over a century]

This is the current best knowledge:

Quit before 40 and chances are that your smoking has no consequences.
Quit before 50 and you shorten your life for -5 years.
Quit before 60 and you shorten your life for -7 years.
Don't quit and you shorten your life for -10 years.

Why Smoking?
Addiction is all about dopamine, a hormone that makes you feel happy. A cigarette gives you a small dose of Nicotine which in turn kickstarts the increase of adrenaline, serotonin and dopamine. You are more aware, more reactive and more happy after smoking. Just a little, not much.
To learn what is good for one-self is crucial for the sophisticated survival machines we are. Therefore a mamals brain learns to seek dopamine rewards. With every cigarette you learn that smoking is good, and the absence of nicotine worries you. The phenomena of withdrawal.
Although the damages it does are way superior to the possible benefits, not everything is bad with smoking: current medicine knowledge is that nicotine helps against alzheimer's and parkinson's. But I also love the nicotine rich tomatoes and potatoes and I started to complete my nutrition with Omega 3 DHA and EPA (makes you smart and protects against the major civilatory deseases like Alzheimer and cardio-vascular deseases).

How to Quit?
Good news: smoking does not provide increasing dopamine returns.
For a binge drinker, every drink augments the desire of another drink.
Even for a chain smoker this is not the case. The craving for a cigarette might not disappear but it certainly does not increase.

Decreasing Returns
The first cigarette is the most wanted.
The addiction develops more as a decrease in the inteval between two smokes.

This is important. We now know the maximum craving is the longing for the first cigarette. We can measure it! What is the unit of measurement? An economist would say: money. I'd say: money is not enough, but a good start. How much would you give me for a cigarette if you really need one?

Would you give me 100 Euros?
Would you give me 50 Euros?
Would you give me 10 Euros?
...

There is an amount which you would pay when craving hard. I decided mine is somewhere around 5 to 10 Euros. This is a measure of your addiction.

Hedging against Smoking
If you want to stop smoking you have to make sure that for each cigarette you smoke, you will lose more value than your addiction treshold. For example you can formulate a contract: "I pay 10 Euros for every cigarette I smoke." I did this, but I did more:

I determined a period for the contract: until 2012-01-01 (Short! Somehow I was not ready for more)
I convinced a friend to be my referee, I gave him 100 Euros to take under ward.
I declared in my contract that I have to give 10 Euros to a detestable political organisation for every cigarette I smoke. This is anti-charity. This is admitting we are not homo economicus - it is is for the social animal in us. Donating the money to an anti-charity pushes the commitment from being only a higher price (I pay for it, so I have the right to smoke it) to something that would weaken your identity if you fail.

My incentives to resist the Nicotine-crave are now bigger than the temptations.
During the first week I felt the urge to smoke strongly. Sometimes I couldn't sleep well, I was nervous. But it really never came to my mind to smoke a cigarette. Its clear that I do not want to donate to the detestable political organisation. The idea of doing that kept me away from the cigarettes. I've externalized my right to smoke by installing a penalty. If it works for the traffic police it's certainly good for keeping me away from bad habits.

The ideas presented here come from Ian Ayres' Carrots and Stick. The idea is genius but the book is a little bit boring, I'd recommend his other book "Supercrunchers".

Which day I lose weight

2011-07-28T21:52:00.003+02:00

I could swear it is weekend's.
But making a 'seasonal' plot, thus aggregating all data for each weekday, showed I do not so good on weekends. Especially Sundays. I am supposed to do sport on sundays. But, hey we had really bad weather the whole July. And some social obligations too...

Is it my day?

Boxes and Wiskers Plot of weekdays: the horizontal line in the middle is the median. Upper box border is the 75% quartile, lower border the 25% quartile,so that 50% of values are between 25-75%. Lines outside are the the next 2-98 percentiles. Dots are outliers. Forget the legend, 0.8 is only the opacity of the boxes.

So in hindsight my worst day is Sunday and my best is Wednesday. How does it relate to kcal left to base rate? How would a weekday regression look like?

A weak start
The regression tells that sunday and monday are fatty
The bigger the dots, the more calories left to base rate.

I even cheated one sunday because of overnight stay in the mountains, including booze, I didn't record.
Disclaimer: You don't know what a mountain is, if you where not there. And I know the evidence is 'thin', five data-points per weekday is not what you should base a scientific law upon

More Quantitative Diet Analysis

2011-07-24T01:22:00.003+02:00

Two Questions:

What is the correlation between calorie consumption and weight loss?
Is it happening the same day?

This Saturday, I skipped hiking up to the mountains to answer these naging questions my last diet analysis left over. To answer them, let me first present the data, carefully recorded over 35 days of dieting.

Weight-loss per day

Measured each morning before breakfast and after taking a crap, pardon my french.

How much did I lose on a regular day?
Most of the time lost between -0.07 and 0.5kg

The Percentiles are:

Min. 1st Q. Median Mean 3rd Q. Max.

-0.66kg -0.07 0.16 0.2 0.5 1.5kg

The 1.5 kg outlier draws the attention. Is it a measurement error?

Net caloric intake per day

This is the amount of calories I have eaten, minus the calories that I have burned doing sport.

When I say calories I mean the american metric which corresponds to the european kilo-calories (kcal).

Caloric Balance Sheet
How much calories did I take-in and consume in a day?
Most often 500 kcal to 960 kcal

The Percentiles are:

Min. 1st Q Median Mean 3rd Q Max.

-1009kcal 490 660 630 960 1480 kcal

My friend, the outlier, is here again. Seems reasonable: A day of heavy exercise and little eating... but 1.5 kg? That's a lot.

Caloric intake left to cover the daily base burn rate

Now its getting nasty. Our body consumes calories even if we lay in bed all day. This is the base rate. Whitout getting into details, it depends roughly on your body weight. When you take-in lesser calories than the base-rate (plus exercising) you lose fat. My formula for this spread goes like this:

kcal_left = base-rate - net_kcal

And the distribution of the not consumed calories is:

Distribution of calories-left
I normally underperform the base-rate by 250kcal - 730kcal.

   Min.   1st Q Median  Mean 3rd Q  Max. 

 -240kcal 250   520     570  730    2230kcal 

Median underperformance is 520kcal, roughly a spurned Pizza Magherita. Two times I've consumed slighty more than the base-rate. And the outlier from my sport weekend is here too!

Correlation between calories left and weight-loss

Back to question 1: What is the correlation between calorie consumption and weight loss?

Now I know that I have to maximize calories-left (to the sedentary base-rate) in order to lose weight.

What are the results of my effort?

Figuring a trend

The Linear regression line in blue

shows there is a positive correlation
between calories-left and weight loss.
The border of the darker region
is the mean deviation oft the data

from the regression line,
a.k.a. the standard error.

If you are a genius you certainly have noticed that the outlier is away. It would have greatly influenced the regression and I felt that I want to answer what happened in normal days. Also the correlation is not high (R^2 = 0.3541857), there are other effects that influence the variables a lot and it's clear that grand part of them are measurement errors.

Coefficients of the regression line are:

(Intercept) kCal_Left

0.1174047 0.0001555

Basically this is the linear function:

weightLoss = 0.1174047 + 0.0001555 * kCal_Left

For the median calorie underperformance of 520 kcal this means:

520 kcal/day * 0.0001555 kg/kcal + 0.1174 kg/day ~= 0.2 kg/day

Not bad at all, and its supported by the fact that really I lost 7kg in 34days. 7/34 = 0.20!

Do I lose weight the same day that I 'starve'?

Question 2 is remaining: Is it happening the same day?

Short answer: Yes. Long answer: Yes, but there might be a lag and prolonged effect for sport activities. The evidence in a picture:

Direct impact of 'fasting'
Size of points representing calories left to base-rate
The higher the line, the greater the weight-loss.
A calorie point sits on the weight-loss recorded the morning after.
That is: the loss-spike is caused by the calorie point it carries.

Often a bigger point sits on a higher spike, which leads to the satisfactory conclusion that indeed you lose the kilos the same day as you restrict your calorie intake. The big outlier we saw already before has a different behavior. Here, the big balls hang relatively low. Does sport have a weight loss lag? The mythical regeneration phase?

The days my sport activity burned more than 1000kcal:

Day 12: 1587 kcal

Day 13: 1946 kcal

Day 29: 1099 kcal

This big outlier on day 12 and 13 could cause a after-burner effect, but the evidence is thin.

More data needed!

Exponentially Smoothed Diet

2011-07-21T00:31:00.001+02:00

Body weight prediction based on irregular measurements over five weeks.
Green is the predicted path in the next ten days.

Exponential smoothing is a technique that can be applied to time series data, either to produce smoothed data for presentation, or to make forecasts. The time series data themselves are a sequence of observations. The observed phenomenon may be an essentially random process, or it may be an orderly, but noisy, process. Whereas in the simple moving average the past observations are weighted equally, exponential smoothing assigns exponentially decreasing weights over time.
[http://en.wikipedia.org/wiki/Exponential_smoothing]

I applied Exponential smoothing to the time series I obtained from my livestrong.com recordings. Of course I have the highest motivation to beware this noble technique from being wrong.
Find all the R code and data below.

Lies, Damn Lies and Diets

2011-07-16T16:31:00.008+02:00

The simplest diet that works is to eat not more calories than you spend. Each day. I experimented with livestrong's myplate, a slick calorie counter.

You enter the food that you ate - "penne alla arrabiata - 100g". You enter the exercise or tasks - "computer work - 480 min" and you can instantly figure out your remaining calorie budget.
It's addicting!
Furthermore it's an effective example how a quantitative approach resolves a problem.

My strongest believe is that one should pay more attention to the genuine scientific question: "How much it is?" than to the genuine philosophical question: "What is it?". Size does matter. All diets are deeply 'philosophical' in this sense, they try to answer the question of what to eat instead of how much. Contrary to this approach you have to start with the amount of calories first, and then - in order to eat enough without spending too much of your calorie budget - you learn to choose food that is calorie-cheap. You eat the vegetable dish and not the pizza. When having a calorie rich plate, you simply eat less. Don't be afraid of throwing half of a Magnum Dark Chocolate Icecream away.

Do you see how big the correlation coefficient R must be here? Linear weight loss, sans famine.

Italy's self-financed debt?

2011-07-12T18:57:00.002+02:00

If national debt is financed by its own households, according to common wisdom the exposure is not problematic. Italy has traditionally a high saving rate, and the national banks advertise government bonds appealingly.
Probably a good share of household savings go into BOT, BCT etc.
But the saving rate has diminished from >16% to 7.4% in the last twenty years.
When experts bring up the high saving rates argument, I can't trust them.

http://www.gfmag.com/tools/global-database/economic-data/10396-household-saving-rates.html#axzz1RpkjkzTb

A Foundation of Measurement

2011-06-10T22:14:00.002+02:00

Do you think Math is the science of measurement?
I am not sure about it, but I know a mathematician who is a radical philosopher of measurement and next week I am going to meet him. Thomas Saaty is an Iranian-American professor of mathematics who goes to the core of it. The comparision is both predecessor and generalization of measurement. The pairwise comparison theory of Prof. Saaty is a reconciliation of mathematics and psychology. With the precise understanding of pairwise comparison you can make an abundance of rational decisions. From picking the right software to
Palestine-Israeli peace process. Ok,... at least we try!
See you at http://www.isahp.org/italy2010

The Kondratjev Epedemy

2011-05-17T22:30:00.009+02:00

Somehow this thing called "Kondratjeff Wave" is stubborn. I encounter it frequently when the occassionaly business related speaker crosses my path. The idea of Nicolai Kondratjev was that the regular business cycles of boom and bust are overlayed by sixthy year super cycles. The problem with business cycles is that they can not be predicted, we know that there is a boom and then eventualy a bust. We don't know how long a boom stays and when and how long the bust takes. Thats a fact, otherwise the market actors would antecipate each phase change and therefore modify the phases. You see, its unpredictability is native to the markets. Therefore "cycle" might be a term taken with a grain of salt: its periodicity is not constant nor predictable. Now Kontradieff made his proposal of Supercycle 1929 and he figured out three cycles in history. The observation of an invariance of three cycles really means nothing statistically speaking. Take the US GDP per capita since 1961, and tell me where the cycles are?

The six Kondratieff Cycles
we should have witnessed in the past century

Kondratieff Cycles - hard to spot: remember, they are sixty year periods. It just doesn't fit

Ok, If this is convincing, why the Kondratieff spreads like a virus?

A regular cycle of boom and bust is a persuasive idea.
We are acquainted with the notion of "epoch" and Kondratieff is a visual representation of the "epoch". Therefore by analogy we understand the concept, it is sound. But conclusions drawn from analogy are not waterproof logical operations. They are merely hints.
The host of the meme has the opportunity to show technical skills (it looks kind of scientific) and demonstrate deeper insight into long term mechanics of human condition. What an opportunity for the host.
The witnesses of the meme have the succumb to 1 and 2 and humans are such pattern addicts, they can not resist the temptation to see a regularity.

We need to implant a new meme that represents the truth:

Kontradieff is bogus, they want to fool you with business esoterics.

I have hope for this brand new meme, humans are ashamed of being tricked.

How to decide when more than one criteria is important: ELECTRE

2011-05-15T10:44:00.000+02:00

The French School of decision making has a short answer: Apply an ELECTRE algorithm.
Howto:

Decide what to compare, e.g. products, actions
Define criteria and their scales, importance and veto tresholds (I'd rather die, than choosing a hotel without Internet connection)
Measure them
Construct outranking relations, in the style of a better b if. These are called Recommandations.
Devise an exploitation procedure that uses the Recommandations to e.g. rank the products.

Sounds promising, after a day of research I must give up because the material found on internet is hard to grasp. But the americans, much more elegalitarian when offering ideas, are doing good job with Analytical Hierarchical Process. Subject under investigation and my next post.

Expert Bias and Making a Screencast 'better'

2011-04-14T23:16:00.001+02:00

I've updated the screencast tutorial for EPR. Even if I augmented the quality of the recording, I have the bad feeling that I've complicated the explaination. Dan and Chip Heaths "Made to Stick" argues that experts fail to teach, because they focus too much on 'important' details. So eventually I have to do it again.

How to Solve It! A book review

2011-04-10T14:46:00.011+02:00

"Schule des Denkens" means school of thought. The title of the original manuscript reveals the motivation for this book: to teach better. George Polya wrote it as a guide on how to lecture mathematics. But only the publicized edition, after Polya's migration to the United States, tells in plain English why its relevant for a broader audience: How to Solve It. I am not a teacher, the didactic musings of a lecturer would produce nothing more than a spark of curiosity in me. But "How to Solve It" describes a general procedure for problem solving. This freaks me out! Polya hides this endeavor in favor of the didactic justification, and only later in this book he will reveal that it is about a Modern Heuristic, to "understand the process of solving problems, especially the mental operations typically useful in this process" [p. 129]. The book is divided into four parts, I will discuss them briefly:

Part I: In the Classroom

This 32 pages are all you need to grasp his algorithm of problem solving. Good ideas are simple and the procedure proposed is not counterintuitive. You could easily come a similiar conclusion by your own:

What is the first step of problem solving?

(1) Understanding the problem.

What is the next step?

(2) Devise a plan.

Then?

(3) Execute the plan.

And finally

(4) Look back

Mightily impressed? Then you are a lobotomized PowerPoint disciple! But follow Polya in a Socratic dialog with the classroom and look into the train of thoughts of an educated problem solver. There are many subtilities to discover.

By reading this chapter, more than once I had moments of Heureka!, when Polya guides you to ask the so called right questions and instructs you how to take a different point of view of the problem.

Part II: How to Solve it - a Dialog

The second part compresses the problem solving procedure, the ars inveniendi, in a summary of two pages. I did not gain from this, but it might be helpful as a short rehearsal when time passes by.

Part III: Dictionary of Heuristic

This is a 200 pages collection of heuristics to use as a pattern language for problem solving. The autor advices to take your time read this piece by piece when you are struggling with problems. Which I do.

Part IV: Problems, Hints, Solutions

These 8 pages are filled with exercises and smart hints how to approach the individual problems.

Conclusions

Polya opens your mind for solutions. I will tackle future hard problems only with Polya's algorithm and benefit from the careful order he imposes to the confused mind.

Albeit written for teaching mathematics, I suspect that Polyas work is useful not only for quantitativ problem solving, but for qualitative problems too. Here, I have no proof and only the application of it will tell.

The book is easy to read, and you might master the first fundamental part "In the Classroom" in 3-5 hours. It is also a really cheap book, 13-something Euros, and if you are a problem solver you will need and enjoy it.

Butter in Bangladesh as a S&P 500 predictor

2011-03-31T10:10:00.002+02:00

http://nerdsonwallstreet.typepad.com/my_weblog/files/dataminejune_2000.pdf

Fully Automated Prediction with Random Forests

2011-03-30T14:28:00.002+02:00

Put the input vector down each of the trees in the forest. Each tree gives a classification, and we say the tree "votes" for that class. The forest chooses the classification having the most votes (over all the trees in the forest).

Basically Random Forests automatically generate many decision trees with mostly weak predictive goodness, and gain high predictive power by averaging them out. The algo can be sketched like this:

randomly sample variables and predictors, repeat:

identifying a predictor
repeat down the tree

seeking the most correlated variable
make a binary decision out of it

combine all predictions and average them out, voila!

Result: high predictive goodness sans parameters!

http://stat-www.berkeley.edu/users/breiman/RandomForests/

Switching of the Nuclear Powerplants in Germany now?

2011-03-15T18:03:00.009+01:00

Today chancellor Merkel announced to switch of 7 nuclear power plants as a consequence of the fukushima catastrophe.
The decision means cutting 10% to 15% of Germanys power consumption, using the net power production median of the weakest nuclear powerplants (890 MW) and the median of the overall powerplants (1288 MW)

Update: they named the plants in question

de.wikipedia.org/wiki/Liste_der_Kernreaktoren_in_Deutschland

Unterweser	KKU	NI	E.ON	1.410	1.345

Biblis B	KWB B	HE	RWE	1.300	1.240
Biblis A	KWB A	HE	RWE	1.225	1.167
Philippsburg 1	KKP 1	BW	EnBW	926	890
Isar/Ohu 1	KKI 1	BY	E.ON	912	878
Neckarwestheim 1	GKN 1	BW	EnBW	840	785
Brunsbüttel	KKB	SH	Vattenfall	806	771

These have a net power of 771+785+878+890+1167+1240+1345= 7076 MW or 11.3% of annual consumption.

Update: the last data on power trade balance shows a surplus of 22 TWh

Hot Bozen

2011-03-14T21:54:00.008+01:00

Since I study time series analysis, here a little exercise with data from Hydrographisches Amt South Tyrol:

Its getting hotter:
decomposing mean temperature in Bozen 1981-2010,
the observed data is splitted into a seasonal component and a level

The code for R:

bz<-read.csv("data/bozen.csv.txt", header=FALSE)

bz<-ts(bz[,1], start=c(1981,1), frequency = 12)

plot(decompose(bz))

+2°C is the mean temperature increase since 1971