Data Regression Tools - Week 2

Running Linear Regression with mean of 0

Summary

The results of the linear regression model show that identifying as a born-again Christian is significantly P<0.00014) and negatively associated (Beta =-4.9709) with the number of times marijuana is smoked in the lifetime of the adolescent respondents.

Prep work

The data I had originally chosen was not suitable for this activity so I chose the categorical explanatory variable "Are you a born-again Christian" (H1RE5) and the quantitative response variable "What is the total number of times you have smoked marijuana in your life?" (H1TO31).

Respondents who identified as a religion other than Christian or those who had not identified themselves as religious were removed from the cleaned dataset. Additionally, respondents who had reported never smoking marijuana (response 997) were changed to 0 and refused and don't know responses (996 and 998) were changed to NA.

The P value (.00014) was very low indicating a significant probability that a relationship exists. The R-square value (.004), however, was also very low which means the explanatory variable is not a good overall predictor of the response variable. The original mean was subtracted from the values in the explanatory variable column to set the mean of the explanatory variable to 0.

The output

Statistics for H1TO31 - Times smoke marijuana (Lifetime)

count    3603.000000

mean        6.764918

std        39.185240

min         0.000000

25%         0.000000

50%         0.000000

75%         0.000000

max       900.000000

Name: H1TO31, dtype: float64

Statistics for religious people who claim to be born again 

(Adjusted Mean)

count    3.603000e+03

mean    -2.434294e-17

std      4.998146e-01

min     -4.840411e-01

25%     -4.840411e-01

50%     -4.840411e-01

75%      5.159589e-01

max      5.159589e-01

Name: H1RE5, dtype: float64

The mean of the explanatory variable is 0.

Count of Christians who identify as 'Born Again'

0- Does not identify, 1-Identifies

-0.484041    1859

 0.515959    1744

Name: H1RE5, dtype: int64

                            OLS Regression Results                            

==============================================================================

Dep. Variable:                 H1TO31   R-squared:                       0.004

Model:                            OLS   Adj. R-squared:                  0.004

Method:                 Least Squares   F-statistic:                     14.54

Date:                Tue, 20 Dec 2016   Prob (F-statistic):           0.000140

Time:                        08:29:54   Log-Likelihood:                -18322.

No. Observations:                3603   AIC:                         3.665e+04

Df Residuals:                    3601   BIC:                         3.666e+04

Df Model:                           1                                         

Covariance Type:            nonrobust                                         

========================================================================

                 coef    std err          t      P>|t|      [95.0% Conf. Int.]

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

Intercept      6.7649      0.652     10.382      0.000         5.487     8.042

H1RE5         -4.9709      1.304     -3.813      0.000        -7.527    -2.415

========================================================================

Omnibus:                     6202.059   Durbin-Watson:                   2.016

Prob(Omnibus):                  0.000   Jarque-Bera (JB):          5264759.670

Skew:                          11.928   Prob(JB):                         0.00

Kurtosis:                     188.742   Cond. No.                         2.00

=======================================================================

Warnings:

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

The code

Prepping the data

1. import pandas as pd  
2. import numpy as np  
3. import matplotlib.pyplot as plt  
4. import statsmodels.formula.api as smf  
5.   
6. filename = ('addhealth_pds.csv')  
7. data = pd.read_csv(filename, low_memory=False)  
8. pd.set_option('display.max_rows', None, 'display.max_columns', None)  
9. #%%  
10. #This section cleans up the data  
11.   
12. data['H1RE5']=data['H1RE5'].replace(8, np.nan)  
13. data['H1RE5']=data['H1RE5'].replace(7, np.nan)  
14. data['H1RE5']=data['H1RE5'].replace(6, np.nan)   
15. data['H1RE5']=data['H1RE5'].replace(9, np.nan)  
16. data['H1TO31']=data['H1TO31'].replace(996, np.nan) #change "refused" to NA  
17. data['H1TO31']=data['H1TO31'].replace(998, np.nan)  
18. data['H1TO31']=data['H1TO31'].replace(997, 0)  
19. data['H1TO31']=data['H1TO31'].replace(999, np.nan)  
20.   
21. clean = data[['H1RE5', 'H1TO31']].dropna()  

Generating statistics

1. clean['H1RE5']=clean['H1RE5']-(clean['H1RE5'].mean())  
2. bornagain = clean['H1RE5'].value_counts()  
3. mjstats =clean['H1TO31'].describe()  
4. print("Statistics for H1TO31 - Times smoke marijuana (Lifetime)\n")  
5. print(mjstats)  
6.   
7. bornstats=clean['H1RE5'].describe()  
8. print("Statistics for religious people who claim to be born again \  
9. \n(Adjusted Mean)")  
10. print(bornstats)  
11. print("\n\nThe mean of the explanatory variable is "\  
12.       + str(int(clean['H1RE5'].mean())) + ".\n\n")  
13.   
14.   
15. print("Count of Christians who identify as 'Born Again'")  
16. print("0- Does not identify, 1-Identifies")  
17. print(bornagain)  
18.   
19. model = smf.ols('H1TO31 ~ H1RE5', data=clean).fit()  
20. print(model.summary())  

Data Regression Tools - Week 1

Describing the data

The sample of 6504 students in 7^th to 12^th grade was taken from 80 high schools and 52 of their feeder schools from across the United States. The ethnic composition was predominantly White (N=4038, 62.3%), Black (N=1483, 22.9%), Asian(N=212, 3.3%), Native American/Indian (N=68, 1.0%), Other(N=362, 5.6%), and Multiple(N=322, 5.0%). The gender composition was Female (N=3356, 51.6%), Male(N=3147, 48.4%) and one respondent was unknown.

The school-based clustered sampling of the 132 school was primarily conducted by Quality Education Data Inc. through the use of in-school questionnaires, in-home interviews, school administrator questionnaires, in-home interviews, and parent questionnaires.

Grade retention was drawn from the question “Were you ever retained or held back a grade?” Students who responded “Yes” were then asked which grade(s) they had been retained. The data was combined into one variable and binned into four categories, “Never held back”, “Elementary School”, “Secondary” and “Multiple”

Violence was drawn from three different questions with the stem “How often the last 12 months have you…”

• ·         Gotten into a fight?
• ·         Pulled a knife or gun on someone?
• ·         Shot or stabbed someone?

Respondents who indicated any level of involvement in any of the three categories were listed as positive for violence while all who had not indicated any involvement were coded as negative for violence.

The age variable was calculated by subtracting the birth year from the year the survey was completed. Since the age variable itself is not used in the calculations, it's exact calculation isn't pertinent. Three respondents had a negative age calculation using this method and that data was disregarded.

Code to extract the above information (except pie chart)

1. import pandas as pd  
2. import numpy as np  
3. import matplotlib.pyplot as plt  
4.   
5.   
6. filename="addhealth_pds.csv"  
7. data = pd.read_csv(filename, low_memory=False)  
8.   
9.   
10. data['H1GI6A']=data['H1GI6A'].replace(6, np.nan)  
11. data['H1GI6A']=data['H1GI6A'].replace(8, np.nan)  
12. data['H1GI6B']=data['H1GI6B'].replace(6, np.nan)  
13. data['H1GI6C']=data['H1GI6C'].replace(6, np.nan)  
14. data['H1GI6D']=data['H1GI6D'].replace(6, np.nan)  
15. data['H1GI6E']=data['H1GI6E'].replace(6, np.nan)  
16.   
17. data['H1GI6B']=data['H1GI6B'].replace(8, np.nan)  
18. data['H1GI6C']=data['H1GI6C'].replace(8, np.nan)  
19. data['H1GI6D']=data['H1GI6D'].replace(8, np.nan)  
20. data['H1GI6E']=data['H1GI6E'].replace(8, np.nan)  
21. data['BIO_SEX']=data['BIO_SEX'].replace(6, np.nan)  
22.   
23. data['NUMETHNIC']=(data['H1GI6A']+ data['H1GI6B'] + data['H1GI6C'] + data['H1GI6D']  
24.  + data['H1GI6E'])  
25.   
26. data["AGE"]=data['IYEAR'] - data['H1GI1Y']  
27. #data["AGE"]=data["AGE"].replace(-1, np.nan)  
28.   
29. age=data['AGE'].value_counts()  
30. print(age)  
31. agedesc=data['AGE'].describe()  
32. print('\nStatistical Description of Age')  
33. print(agedesc)  
34.   
35.   
36. genderlabel={1:"Male", 2:'Female', 6:'unknown'}  
37. data['BIO_SEX']=data['BIO_SEX'].map(genderlabel).astype('category')  
38. gender = data['BIO_SEX'].value_counts()  
39.   
40. genderprc =data['BIO_SEX'].value_counts(normalize=True)  
41. print("\nCount of Gender")  
42. print(gender)  
43. print("\nPercent of each Gender")  
44. print(genderprc)  
45.   
46. def ethnic (row):  
47.     if row["NUMETHNIC"] > 1:  
48.         return 6  
49.     if row['H1GI6A'] == 1:  
50.         return 1  
51.     if row['H1GI6B'] == 1:  
52.         return 2  
53.     if row['H1GI6C'] == 1:  
54.         return 3  
55.     if row['H1GI6D'] == 1:  
56.         return 4  
57.     if row['H1GI6E'] == 1:  
58.         return 5  
59.     else:  
60.         return np.nan  
61.   
62. ethniclabels = {1:'White', 2:"Black", 3:'Native American', 4:"Asian", 5:'Other',  
63.                 6:'Multiple'}  
64. data['ETHNICITY'] = data.apply(lambda row: ethnic (row), axis=1)  
65.   
66. data['ETHNICITY']=data['ETHNICITY'].map(ethniclabels)  
67. subwhite = data[(data['ETHNICITY'] == 'White')]    
68. whiteage = subwhite['AGE'].tolist()   
69.   
70. subblack = data[(data['ETHNICITY'] == 'Black')]    
71. blackage = subblack['AGE'].tolist()   
72.   
73. subnative = data[(data['ETHNICITY'] == 'Native American')]    
74. nativeage = subnative['AGE'].tolist()  
75.   
76. subasian = data[(data['ETHNICITY'] == 'Asian')]    
77. asianage = subasian['AGE'].tolist()  
78.   
79. subother = data[(data['ETHNICITY'] == 'Other')]    
80. otherage = subother['AGE'].tolist()  
81.   
82. submulti = data[(data['ETHNICITY'] == 'Multiple')]    
83. multiage = submulti['AGE'].tolist()  
84.    
85. ethnicities=[whiteage, blackage, nativeage, asianage, otherage, multiage]   
86.   
87. ethnicity = data['ETHNICITY'].value_counts()  
88. ethnicprc =data['ETHNICITY'].value_counts(normalize=True)  
89. print("\nCounts for Ethnic Groups")  
90. print(ethnicity)  
91. print("\nPercentage of each Ethnic Group")  
92. print(ethnicprc)  
93.   
94. #clean=data[['AGE', 'ETHNICITY', 'BIO_SEX']].dropna()   
95. plt.boxplot(ethnicities, labels = ('White', 'Black', 'Nat. Am.',  
96.                                    'Asian', 'Other', 'Multiple'))  
97. plt.xlabel("Ethnicity")  
98. plt.ylabel("Age")  
99. plt.title("Age Range of ADDHEALTH Respondents")  
100. plt.show()  

The output (boxplot is shown above)

Ages (Count)

17    1169

 16    1162

 18    1149

 15    1081

 14     903

 13     578

 19     389

 20      43

 21      18

 12       9

-1        3

Name: AGE, dtype: int64

Statistical Description of Age

count    6504.000000

mean       16.028905

std         1.810004

min        -1.000000

25%        15.000000

50%        16.000000

75%        17.000000

max        21.000000

Name: AGE, dtype: float64

Count of Gender

Female    3356

Male      3147

Name: BIO_SEX, dtype: int64

Percent of each Gender

Female    0.51607

Male      0.48393

Name: BIO_SEX, dtype: float64

Counts for Ethnic Groups

White              4038

Black              1483

Other               362

Multiple            322

Asian               212

Native American      68

Name: ETHNICITY, dtype: int64