GapMinder Blog - Association between Internet Use and Suicide

This is a blog wherein I record and submit my weekly progress towards the Data Analysis and Interpretation Specialization on Coursera

[reetikasarkar-blog]

Regression Modelling in Practice (Week 4 : Testing a Logistic Regression Model)

The aim of this week is to test a logistic regression model. For my anaysis, I used the same response variable from the previous week, i.e., suicideper100th. My primary regressor was alcconsumption. Since the dataset from GapMinder comprises only quantitative variables, data management had to be done before carrying out the analysis. The following were done:

(i) Since the response variable is quantitative, I first binned it into two categories - low (=0) and high (=1), based on mean. 

(ii) The next step was to bin a quantitative predictor variable into categories for which I chose income levels (incomeperperson variable). I binned it into 4 categories based on the quartiles. 

(iii) The remaining predictors were centered based on their respective means and the means were checked to verify that the means had shifted close to 0.

NOTE: I have not chosen to bin the quantitative predictors for testing logistic regression and the interpretations for odds ratios have been done accordingly.

For the analysis, the SAS Code I used is as under:

libname mydata "/courses/d1406ae5ba27fe300/" access=readonly;

data new; set mydata.gapminder;

label suicidegroup = "Level of Suicides";

******************************************************************** DATA MANAGEMENT ********************************************************************; * binning response variable; if suicideper100th eq . then suicidegroup = .; else if suicideper100th le 8 then suicidegroup = '0'; /* low */ else suicidegroup = '1'; /* high */

* binning predictor variable incomeperperson into 4 categories; if incomeperperson eq . then incomegroup = .; else if incomeperperson le 750 then incomegroup = 0; else if incomeperperson le 3000 then incomegroup = 1; else if incomeperperson le 10000 then incomegroup = 2; else incomegroup = 3;

run;

* centering quantitative predictors; data new2; set new; if suicidegroup ne . and internetuserate ne . and alcconsumption ne . and   urbanrate ne . and incomegroup ne .; internetuserate_c=internetuserate - 34.0510892; alcconsumption_c=alcconsumption - 6.8731395; urbanrate_c=urbanrate - 55.2483721; run;

proc means; var alcconsumption internetuserate urbanrate; run;

* check means of centered variables; proc means; var internetuserate_c alcconsumption_c urbanrate_c; run;

*adding 4 category incomegroup variable with reference group as 0; PROC GLM; class incomegroup (ref="0"); model suicidegroup=alcconsumption_c internetuserate_c urbanrate_c incomegroup/solution; run;

* change reference group to 1; PROC GLM; class incomegroup (ref="1"); model suicidegroup=alcconsumption_c internetuserate_c urbanrate_c incomegroup/solution; run;

*********************************************************************** LOGISTIC REGRESSION ***********************************************************************;

PROC LOGISTIC descending; model suicidegroup=alcconsumption; run;

PROC LOGISTIC descending; model suicidegroup=internetuserate; run;

PROC LOGISTIC descending; model suicidegroup=urbanrate;                run;

* adding internetuserate and urbanrate; PROC LOGISTIC descending; model suicidegroup=alcconsumption internetuserate urbanrate; run;

ANALYSIS

1) Centering the explanatory variables

The means procedure returned the means of the predictors before and after centering and it was seen that post centering, the means of the regressors shifted very close to 0.

2) Reference Group Coding/Parameterization for incomegroup variable

First the chosen predictor incomeperperson was binned into 4 categories to form a new variable incomegroup with levels 0 - low income, 1 - moderate income, 2 - high income and 3 - very high income.

The GLM procedure was used to examine the explanatory variable differences on the response variable after adjusting for other explanatory variables in the model. 

The model predicted that alcconsumption_c and urbanrate_c were the predictors that significantly effected suicide levels. This was further verified on fitting a logistic regression model and observing the fit statistics, as shown later.

First each of the income groups was compared to the reference group = 0 to test whether low income has anything to do with high suicide (suicidegroup=1), after controlling for other predictors in the model. The ‘Estimate’ column in the table below depicts the comparison between a) incomegroup0 and incomegroup1, b) incomegroup0 and incomegroup2, and c) incomegroup0 and incomegroup3. All the three comparisons had non-significant p-values at 5% level of significance indicating that none of the income groups have a significant effect on high suicide levels, in comparison to the low income group (incomegroup=0). 

Then income groups were compared to the reference group = 1 to test whether moderate income has anything to do with high suicide (suicidegroup=1), after controlling for other predictors in the model. The ‘Estimate’ column in the table below depicts the comparison between a) incomegroup1 and incomegroup0, b) incomegroup1 and incomegroup2, and c) incomegroup1 and incomegroup3. All the three comparisons had non-significant p-values at 5% level of significance indicating that none of the income groups have a significant effect on high suicide levels, in comparison to the moderate income group (incomegroup=1).

Note:The X'X matrix has been found to be singular, and a generalized inverse was used to solve the normal equations. Terms whose estimates are followed by the letter 'B' are not uniquely estimable.

3) Logistic Regression

A) alcconsumption as the only predictor

Alcohol consumption has a significant positive effect on the probability of having a high suicide level (Coeff. = 0.1542 and p-value < 0.05), when it is the only predictor being fitted.

The odds ratio estimate is 1.167 and there is a 95% chance that the predicted odds ratio will lie between 1.085 and 1.255 in subsequent samples drawn from the population. The coefficient of 1.167 (>1) indicates that for an increase of 1 litre in alcohol consumption, the odds of having a high level of suicide increases by a factor 1.167.

B) internetuserate as the only predictor

Internet use does not have a significant effect on the probability of having a high suicide level (Coeff. = 0.00872 and p-value(=0.1222) > 0.05), when it is the only predictor being fitted.

The odds ratio estimate is 1.009 and there is a 95% chance that the predicted odds ratio will lie between 0.998 and 1.020 in subsequent samples drawn from the population. The coefficient of 1.009 is almost equal to 1 and indicates that for an increase of 1 unit in internet use rate, the odds of having a high level of suicide remains the same. Also, the 95% confidence interval for the odds ratio is very close to 1 so internet use cannot be said to have a significant effect on suicide levels.

C) urbanrate as the only predictor

Urbanisation rate does not have a significant effect on the probability of having a high suicide level (Coeff. = -0.0125 and p-value(=0.0680) > 0.05), when it is the only predictor being fitted.

The odds ratio estimate is 0.998 and there is a 95% chance that the predicted odds ratio will lie between 0.974 and 1.001 in subsequent samples drawn from the population. The coefficient of 0.998 is almost equal to 1 and indicates that for an increase of 1 unit in rate of urbanisation, the odds of having a high level of suicide remains the same (slight decrease as parameter estimate is negative). Also, the 95% confidence interval for the odds ratio is very close to 1 so urbanisation cannot be said to have a significant effect on suicide levels.

D) Model based on alcconsumption, internetuserate and urbanrate

In the model consisting of all the above three predictors, the coefficients and the p-values obtained are: 

> alcconsumption (coefficient = 0.1791, p-value = 0.0002)

> internetuserate (coefficient = 0.0124, p-value = 0.1933)

> urbanrate (coefficient = -0.0356, p-value = 0.0007)

So at 5% level of significance, alcconsumption and urbanrate are the predictors that are significantly associated with suicide level. There is thus a deviation in the estimates corresponding to the predictor urbanrate as can be seen from part C wherein it was observed that it does not have significant association with suicide level. This is indicative of a confounding effect of additional regressors on urbanrate.

The odds ratio estimates and the corresponding 95% confidence intervals are as under:

> alcconsumption (estimate = 1.196, CI = (1.089, 1.314))

The coefficient of 1.196 (>1) indicates that for an increase of 1 litre in alcohol consumption, the odds of having a high level of suicide increases by a factor 1.167. Subsequent samples will include the odds ratio estimate within the above CI in 95% of samples.  

> internetuserate (estimate = 1.012, CI = (0.994, 1.032))

The coefficient of 1.012 is almost equal to 1 and indicates that for an increase of 1 unit in internet use rate, the odds of having a high level of suicide remains the same. Also, the 95% confidence interval for the odds ratio is very close to 1 so internet use cannot be said to have a significant effect on suicide levels.

> urbanrate (estimate = 0.965, CI = (0.945, 0.985))

We cannot not make any inference based on the above CI for urbanrate as its effect is getting confounded by inclusion of additional predictor in the model.

CONCLUSIVE SUMMARY

For my chosen response (suicidegroup) and primary regressor (alcconsumption) there appears to be a significant positive association between and addition of additional predictors, viz., internetuserate and urbanrate does not affect the relationship. So there is no evidence of confounding for the association between my response and primary predictor variable.

Among the three chosen predictors, only alcohol consumption appears to have a significant effect on suicide levels in a country and the effect is of positive nature, i.e., higher levels of alcohol consumption are associated with higher levels of suicide.

[reetikasarkar-blog]

Regression Modelling in Practice (Week 3 : Testing a Multiple Regression Model)

The aim this week is to fit and test a multiple regression model for a given response variable. Moving on from basic linear fitting done on my response variable suicideper100th last week, this week I first tested the fit of a simple linear model based on internetuserate and subsequently on a quadratic model based on internetuserate. I also tested the fit using an additional predictor alcconsumption to explore the association between suicide rate and alcohol consumption, while controlling for other factors. 

For testing the model, we set up

Null hypothesis (H0): there is no association between the response variable and the predictor, against 

Alternate Hypothesis (H1): there is significant association between the response and predictor.

The SAS Code used is as under:

libname mydata "/courses/d1406ae5ba27fe300/" access=readonly;

data new; set mydata.gapminder;

label internetuserate_c = "Centered Internet Use Rate"         alccconsumption_c = "Centered Alcohol Consumption";

run;   ********************************************************************** POLYNOMIAL REGRESSION **********************************************************************;

* scatterplot with linear regression line suicideper100th response variable; proc sgplot;  reg x=internetuserate y=suicideper100th / lineattrs=(color=blue thickness=2) clm;  yaxis label="Suicide Rate (per 100,000 persons)";  xaxis label="Internet Use Rate"; run;

* scatterplot with linear and quadratic regression line; proc sgplot;  reg x=internetuserate y=suicideper100th / lineattrs=(color=blue thickness=2) degree=1 clm;  reg x=internetuserate y=suicideper100th / lineattrs=(color=green thickness=2) degree=2 clm;  yaxis label="Suicide Rate (per 100,000 persons)";  xaxis label="Internet Use Rate"; run;

* centering quantitative explanatory variables; data new2; set new; if suicideper100th ne . and internetuserate ne . and alcconsumption ne .; internetuserate_c=internetuserate-33.7858281; alcconsumption_c=alcconsumption-6.7778409; run;

proc means; var internetuserate alcconsumption; run;

* check coding; proc means; var internetuserate_c alcconsumption_c; run;

* linear regression model; PROC glm;  model suicideper100th=internetuserate_c/solution clparm; run;

PROC glm;  model suicideper100th=alcconsumption_c/solution clparm; run;

* polynomial regression model; PROC glm;  model suicideper100th=internetuserate_c internetuserate_c*internetuserate_c/solution clparm; run;

********************************************************************** EVALUATING MODEL FIT **********************************************************************;

* multiple regression adding employment rate; PROC glm;  model suicideper100th=internetuserate_c internetuserate_c*internetuserate_c alcconsumption_c/solution clparm; run;

* request regression diagnostic plots; PROC glm PLOTS(unpack)=all; model suicideper100th=internetuserate_c internetuserate_c*internetuserate_c  alcconsumption_c/solution clparm; output residual=res student=stdres out=results; run;

* plot of standardized residuals for each observation; proc gplot; label stdres="Standardized Residual"; plot stdres*country/vref=0;  run;

* using proc reg to get a partial regression plot; * calculate quadratic terms; data partial; set new2; internetuserate2=internetuserate_c*internetuserate_c; run;

*partial regression plot;

PROC reg plots=partial;

model suicideper100th=internetuserate internetuserate2 alcconsumption/partial;

run;

ANALYSIS

1) Centering the predictors

The means of the predictors internetuserate and alcconsumption before centering were 33.7641864 and 6.7778409 respectively. The means after centering were reduced to almost 0 for both the predictors.

2) Testing the regressors

A) For testing suicide rate and internet use rate, the estimated parameter coefficient was 0.0162 (positive) and the corresponding p-value (0.3433) > 0.05, so at 5% level of significance we do not have sufficient evidence to conclude that there is a significant association between internet use and suicide. Figure 1 shows that most of the observations were within the 95% prediction limits, barring those above 21 (approx).

Figure 1. Fit Plot for suicide rate based on internet use rate (centered)

B) For testing suicide rate and alcohol consumption, the estimated parameter coefficient was 0.497 (positive) and the corresponding p-value (<0.0001) < 0.05, so at 5% level of significance we can conclude that there is a significant association between alcohol consumption and suicide. Figure 2 shows that most of the observations were within the 95% prediction limits and were crowded around values of alcohol consumption (centered) lower than 10.

Figure 2. Fit Plot for suicide rate based on alcohol consumption (centered)

C) For testing the fit of a quadratic model based on internetuserate, from Figure 3 and Figure 4, it was observed that for both the predictors (internetuserate_c and internetuserate_c*internetuserate_c), most of the observations were outside the 95% confidence limits. Also, the observations were highly scattered around the line of regression, indicating the fit may not be good. This was further substantiated using the results from the GLM procedure according to which the model has a very high p-value (0.2586) indicating that the model coefficients are not significant, which is indeed true as can be inferred from the high p-values of 0.9618 and 0.1792 for internetuserate_c and  internetuserate_c * internetuserate_c, respectively. Thus, quadratic fit based on internetuserate is also not very effective in predicting suicide rates.

Figure 3. Fit Plot for simple linear model based on internetuserate

Figure 4. Fit Plot for quadratic model based on internetuserate

D) For testing the fit of a model based on internetuserate_c,  internetuserate_c * internetuserate_c and alcconsumption_c, the p-value was <0.0001, which indicated that the fit was significant or that the predictors were significantly associated with the response. The final predicted model was 

Y = 8.3536 - 0.0728*X1 + 0.0016*X1*X1 + 0.6649*X2 

where X1 is internetuserate and X2 is alcohol consumption (in litres)

It was thus observed that on addition of a variable alcohol consumption in the model, internet use rate showed a negative association with suicide rate which is contrary to what was observed for a simple linear regression model. Also, earlier the corresponding regression coefficients for internetuserate_c and internetuserate_c*internetuserate_c were not significant and after addition of alcconsumption variable, their association became significant (p-values are 0.0009 and 0.0129, respectively). This indicates that alcohol consumption has a confounding effect on the relationship between internet use rate and suicide rate. Effect of alcohol consumption however, remains significant even on addition of variables to the model indicating a strong association with suicide rates. 

3) Validation of Assumptions

The residuals are symmetrically distributed around the middle of the plot (Figure 5 & 6), i.e., around residual=0 and also no pattern can be observed in the plot indicating that the residuals and predicted values are unrelated or the actual and predicted values of the response variable are somewhat close to each other. The plots of residuals by regressors (Figure 7) however, indicates that although internetuserate_c and alcconsumption_c seem to be uncorrelated with the residuals, for internetuserate_c*internetuserate_c a pattern in the plot was observed indicating differences in the predicted and actual response values.

Figure 5. Plot of Residuals vs. Predicted response

Figure 6. Plot of Studentized Residuals vs. Predicted response

Figure 7. Plots of residuals by regressors

In linear regression, an outlier is an observation with large residual. An observation with an extreme value on a predictor variable is called a point with high leverage. Leverage is a measure of how much an observation has effected its fitted value from the regression model. From the plot of outlier and leverage points (Figure 8), it is seen that points with studentized residuals higher than 2 are considered outliers and points with leverage value higher than 0.045 (approx.) are considered points with high leverage. The data needs to be treated for these before making any conclusive comment on the underlying relationship of the predictors and the response.

Figure 8. Plot of Outliers and Leverage Points

If the residuals can be taken as coming from normal distribution then the assumption of normality of errors is valid. From figure 9, it can be seen that the Q-Q plot of residuals is not linear and a non-linear pattern (curve) exists which indicates deviation from Normality. This can also be ascertained from figure 10 which shows a more peaked curve of the Kernel than that of Normal distribution.

Figure 9. Q-Q plot of residuals

Figure 10. Plot for distribution of Residuals

In the plot of standardized residuals vs. Country, although some of the points lie far from the line of 0 residual, i.e., fitted=predicted, most of the points lie within a distance of  2 units from the line of 0 deviation.

Figure 11. Plot of Standardized Residuals

In order to check whether the variable under consideration affects the response linearly in the presence of other variables, we use partial regression plots (also called added variable plots). If the points are scattered more or less along a straight line in the partial plot then that feature can be taken as linear. From figure 12 it can be seen that internetuserate and alcconsumption have a linear relationship while internetuserate2 (internetuserate^2) displays a non-linear relationship.

Figure 12. Partial Regression Plots by regressors

[reetikasarkar-blog]

Machine Learning for Data Analysis (Week 4 : Running a k-Means Cluster Analysis)

Cluster analysis or clustering is the task of grouping a set of objects in such a way that objects in the same group (called a cluster) are more similar (in some sense or another) to each other than to those in other groups (clusters). 

A k-means cluster analysis was conducted to identify underlying subgroups of countries based on their similarity of responses on 9 variables that represent characteristics that could have an impact on suicide rate. 

ANALYSIS

1. Data Management

The clustering variables namely, incomeperperson, urbanrate, alcconsumption, employrate, co2emissions, oilperperson, internetuserate, polityscore and lifeexpectancy were all quantitative in nature. All clustering variables were standardized to have a mean of 0 and a standard deviation of 1.

The SAS Code for carrying out cluster analysis is as under:

libname mydata "/courses/d1406ae5ba27fe300" access=readonly;

********************************************************************************** DATA MANAGEMENT **********************************************************************************; data clust; set mydata.gapminder;

* create a unique identifier to merge cluster assignment variable with the main data set; country = _n_;

keep country suicideper100th incomeperperson urbanrate alcconsumption employrate co2emissions oilperperson internetuserate polityscore lifeexpectancy;

* delete observations with missing data; if cmiss(of _all_) then delete; run;

ODS GRAPHICS ON;

* Split data randomly into test and training data; proc surveyselect data=clust out=traintest seed = 123 samprate=0.7 method=srs outall; run;

data clus_train; set traintest; if selected=1; run;

data clus_test; set traintest; if selected=0; run;

* standardize the clustering variables to have a mean of 0 and standard deviation of 1; proc standard data=clus_train out=clustvar mean=0 std=1; var incomeperperson urbanrate alcconsumption employrate co2emissions oilperperson internetuserate polityscore lifeexpectancy; run;

%macro kmean(K); proc fastclus data=clustvar out=outdata&K. outstat=cluststat&K. maxclusters= &K. maxiter=300; var incomeperperson urbanrate alcconsumption employrate co2emissions oilperperson internetuserate polityscore lifeexpectancy; run;

%mend;

%kmean(1); %kmean(2); %kmean(3); %kmean(4); %kmean(5); %kmean(6); %kmean(7); %kmean(8); %kmean(9);

* extract r-square values from each cluster solution and then merge them to plot elbow curve; data clus1; set cluststat1; nclust = 1;

if _type_='RSQ';

keep nclust over_all; run;

data clus2; set cluststat2; nclust=2;

if _type_='RSQ';

keep nclust over_all; run;

data clus3; set cluststat3; nclust = 3;

if _type_='RSQ';

keep nclust over_all; run;

data clus4; set cluststat4; nclust = 4;

if _type_='RSQ';

keep nclust over_all; run;

data clus5; set cluststat5; nclust = 5;

if _type_='RSQ';

keep nclust over_all; run;

data clus6; set cluststat6; nclust = 6;

if _type_='RSQ';

keep nclust over_all; run;

data clus7; set cluststat7; nclust = 7;

if _type_='RSQ';

keep nclust over_all; run;

data clus8; set cluststat8; nclust = 8;

if _type_='RSQ';

keep nclust over_all; run;

data clus9; set cluststat9; nclust = 9;

if _type_='RSQ';

keep nclust over_all; run;

data clusrsquare; set clus1 clus2 clus3 clus4 clus5 clus6 clus7 clus8 clus9; run;

* plot elbow curve using r-square values; symbol1 color=blue interpol=join; proc gplot data=clusrsquare; plot over_all*nclust; run;

***************************************************************************************** further examine cluster solution for the number of clusters suggested by the elbow curve *****************************************************************************************

* plot clusters for 4 cluster solution; proc candisc data=outdata4 out=clustcan; class cluster; var incomeperperson urbanrate alcconsumption employrate co2emissions oilperperson internetuserate polityscore lifeexpectancy; run;

proc sgplot data=clustcan; scatter y=can2 x=can1 / group=cluster; run;

* validate clusters on suiciderate;

* first merge clustering variable and assignment data with suiciderate data; data suicide_data; set clus_train; keep country suicideper100th; run;

proc sort data=outdata4; by country; run;

proc sort data=suicide_data; by country; run;

data merged; merge outdata4 suicide_data; by country; run;

proc sort data=merged; by cluster; run;

proc means data=merged; var suicideper100th; by cluster; run;

proc anova data=merged; class cluster; model suicideper100th = cluster; means cluster/tukey; run;

2. Elbow Curve to determine optimum number of clusters

Data were randomly split into a training set that included 70% of the observations and a test set that included 30% of the observations. A series of k-means cluster analyses were conducted on the training data for k varying from 1 to 9 clusters, using Euclidean distance. The variance in the clustering variables that was accounted for by the clusters (r-square) was plotted for each of the nine cluster solutions in an elbow curve (Figure 1) to provide guidance for choosing the number of clusters to interpret.

Figure 1. Elbow curve of r-square values from the clusters

Although there was a sharp change (elbow) for k=2, another elbow for k=7 is also seen. So, the elbow curve was inconclusive, suggesting that the 2 and 7 cluster solutions might be interpreted. The results below are for an interpretation of the 4-cluster solution.

3. Canonical Discriminant Analysis

Canonical discriminant analysis was used to reduce the 11 clustering variable to a few variables that accounted for most of the variance in the clustering variables. A scatterplot of the first two canonical variables by cluster (Figure 2 shown below) indicated that the observations in cluster 1 were somewhat closely packed with moderate within cluster variance. Clusters 2, 3 and 4 were generally distinct, but the observations had greater spread suggesting higher within cluster variance. The clusters did not overlap very much with each other. The results of this plot suggest that the best cluster solution may have fewer than 4 clusters, so it will be especially important to also evaluate the cluster solutions with fewer than 4 clusters.

Figure 2. Plot of the first two canonical variables for the clustering variables by cluster (for 4 clusters)

The scatterplot for first two canonical variables by 3 clusters (Figure 3) indicates that observations in cluster 3 are closely packed and with low variability. Clusters 1 and 2, although distinct, are scattered randomly and there is high variability within these two clusters. Cluster analysis is perhaps not a viable method for grouping the objects in this dataset based on suicide rates.

Figure 3.  Plot of the first two canonical variables for the clustering variables by cluster (for 3 clusters)

The cluster means of the standardized clustering variables (for 4 clusters) show that:

a) cluster 1 had highest magnitude of mean for none of the variables

b) cluster 2 had highest magnitude of mean for income and co2 emissions

c) cluster 3 had highest magnitude of mean for alcohol consumption (negative) and polity score (negative)

d) cluster 4 had highest magnitude of mean for urban rate, employment rate, oil consumption, internet use and life expectancy

Approximated expected overall R square for the model is 0.31531

4. Cluster Means Analysis and ANOVA

The means on the clustering variables showed that on an average, compared to the other clusters, countries in cluster 1 had higher levels of suicides per 100,000 persons (mean=11.64) while cluster 3 had the lowest level (mean=9.18). Cluster 1 had the highest within cluster variability (sd=9.40) and also the highest number of countries from the test data (17) were included in it. Cluster 4 had the lowest within cluster variability (sd=1.44). 

In order to externally validate the clusters, an Analysis of Variance (ANOVA) was conducted to test for significant differences between the clusters on suicide rates (suicideper100th). The null hypothesis was that there is no significant association between the clusters and suicides per 100,000 persons against the alternate hypothesis that the association is significant. Tukey’s test was used for post hoc comparisons between the clusters. 

The R-Square coefficient was 0.066121, i.e., the model based on 4 clusters explained 6.6% of the variability in suicide rate per 100,000 persons. Also, the corresponding p-value based on F-statistic was 0.4760 indicating that there is no association between the clusters and suicide rate. The results of Tukey’s test at 0.05 level of significance also indicated that none of the clusters indicated significant differences with each other.

Boxplots generated under the ANOVA procedure (Figure 4) give a graphical representation of the skewness in the clusters (barring cluster 4, all the remaining clusters 1, 2 and 3 display positive skewness) and also the variability and mean of the number of suicides for countries within each cluster, as already discussed under part 4.

Figure 4. Boxplots indicating distribution of suicides in the 4 clusters

The entire output can be seen here.

[reetikasarkar-blog]

Machine Learning for Data Analysis (Week 3 : Running a Lasso Regression Analysis)

In statistics and machine learning, lasso (least absolute shrinkage and selection operator) (also Lasso or LASSO) is a regression analysis method that performs both variable selection and regularization in order to enhance the prediction accuracy and interpretability of the statistical model it produces.

For my study, I have chosen suicide rate as the response variable and since it is continuous, I have categorised it into two categories, i.e., countries having low suicide rate (<8) coded as 0 and those with higher rates as 1.

A lasso regression analysis was conducted to identify a subset of variables from a pool of 14 quantitative response variables measuring levels of suicide among countries. The predictor variables included internet use rate, income, urbanisation rate, alcohol consumption, CO2 emissions, employment, oil consumption, polity score, life expectancy, armed forces rate, breast cancer prevalence, female employment, HIV prevalence and electricity consumption.

Data were randomly split into a training set that included 70% of the observations and a test set that included 30% of the observations. The least angle regression algorithm with k=10 fold cross validation was used to estimate the lasso regression model in the training set, and the model was validated using the test set. The change in the cross validation average (mean) squared error at each step was used to identify the best subset of predictor variables.

The SAS code for running Lasso regression analysis is as under:

libname mydata "/courses/d1406ae5ba27fe300" access=readonly; data new; set mydata.gapminder;

*********************************************************************************************** DATA MANAGEMENT ***********************************************************************************************; LABEL suicidegroup = "Level of suicide"; IF suicideper100th EQ . THEN suicidegroup = .; ELSE IF suicideper100th LE 8 THEN suicidegroup = '0'; ELSE suicidegroup = '1';

* delete observations with missing data; IF cmiss(of _all_) THEN delete; RUN;

ods graphics on;

* Split data randomly into test and training data; PROC SURVEYSELECT data=new out=traintest seed = 123 samprate=0.7 method=srs outall; run;  

* lasso multiple regression with lars algorithm k=10 fold validation; PROC GLMSELECT data=traintest plots=all seed=123;     partition ROLE=selected(train='1' test='0');     model suicidegroup = internetuserate incomeperperson urbanrate     alcconsumption co2emissions employrate oilperperson polityscore     lifeexpectancy armedforcesrate breastcancerper100th femaleemployrate     hivrate relectricperperson / selection=lar(choose=cv stop=none)    cvmethod=random(10);

RUN;

OUTPUT

The sample for the training data has been drawn using simple random sampling with a selection probability of 0.714286 of each observation.

On applying lasso regression, it was found that all of the variables have been retained in the model. The values of ASE (Mean Square Error) decrease as the number of predictors increase, which is in expected as the model complexity increases. For the test data, a rise in ASE is seen 10th step onward. 

The optimal model has been obtained at the 4th step. The values in the CV Press column first decrease and then sharply rise after the 4th step. The best subset of predictor variables thus includes 4 predictors namely, internetuserate, alcconsumption, urbanrate and oilperperson.

From the graph above it is seen that among the 4 selected predictor variables while oil consumption, internet use and alcohol consumption are positively associated suicide level, urban rate has a negative association.

From the fit criteria based on other methods, it is seen that the best criterion values for AIC, AICC and Adj R-Sq are at the 8th step while for SBC method, it is at the intercept step. These are in contrast with the predicted best subset based on the CV PRESS values.

The ASE for the training data decreases consistently as the number of predictors in the model increases. The ASE for test data shows a comparatively moderate decrease and eventually starts increasing as model complexity increases.

The best predicted model is:

Y = 0.6496 + 0.0036*X1 - 0.0025*X2 + 0.0028*X3 + 0.0043*X4

where Y = suicide level, X1 = internetuserate, X2 = urbanrate, X3 = alcconsumption, X4 = oilperperson

The selected model has 4 degrees of freedom and an F-value of 1.45. Adjusted R square is 0.0444, i.e., the model explains 4.44% of the variability in the response variable.

[reetikasarkar-blog]

Regression Modelling in Practice (Week 2: Basics of Linear Regression)

For conducting linear regression, my response variable is suicideper100th and explanatory variable is internetuserate, which is continuous. I will be conducting the study using SAS in three parts:

1) linear regression without centering the continuous explanatory variable

2)  linear regression after centering the continuous explanatory variable

3) linear regression after recoding continuous explanatory variable to a binary categorical variable (based on mean value)

For both the parts, we set up the hypothesis:

Null Hypothesis, H0: There is no significant association between internet use and suicide rate, against

Alternate Hypothesis, H1: There is a significant association between internet use and suicide rate

My SAS code is as under:

LIBNAME mydata "/courses/d1406ae5ba27fe300" access=readonly; DATA new; set mydata.gapminder;

/* Setting the labels for each of the variables under study */ LABEL internetpenetration="Internet use level" suicideper100th="Suicides per 100,000" internetuserate="Internet users per 100 persons" c_mean_internet="Centered Internet use rate";

/* Centering response variable */ c_mean_internet = internetuserate - 35.6327158;

/* Setting categories for continuous explanatory variable */ IF internetuserate EQ . THEN internetpenetration = .; ELSE IF internetuserate LE 35 THEN internetpenetration = '0'; ELSE internetpenetration = '1';

IF internetpenetration NE .;

/* Obtaining summary of explanatory variable */ PROC MEANS; VAR internetuserate c_mean_internet suicideper100th; RUN;

/*********************************************************************************** PART 1 ************************************************************************************/

/* Scatterplot with regression line for Continuous explanatory variable */ PROC SGPLOT; REG x=internetuserate y=suicideper100th/lineattrs=(color = blue thickness = 2); title "Scatterplot for Association between Internet use rate and Suicide rate"; xaxis label="Internet use rate"; yaxis label="Suicide rate"; RUN; title;

/* Basic Linear Regression */ PROC GLM; MODEL suicideper100th = internetuserate / solution; RUN;

/*********************************************************************************** PART 2 ************************************************************************************/ PROC GLM; MODEL suicideper100th = c_mean_internet / solution; RUN;

/*********************************************************************************** PART 3 ************************************************************************************/

/* Categorical explanatory variable */ PROC GLM; MODEL suicideper100th = internetpenetration /solution; RUN;

/* Means for response variable */ PROC SORT; BY internetpenetration; PROC MEANS; VAR suicideper100th; BY internetpenetration; run;

/* Bar Graph for Response variable by categorical explanatory variable */ PROC GCHART; VBAR internetpenetration / discrete TYPE = mean SUMVAR = suicideper100th; RUN;

QUIT;

OUTPUT

The above table shows that centering the explanatory variable (internetuserate) on the basis of mean, reduces its mean to almost zero.

Part 1

From the scatterplot, it can be seen that the points are very scattered around the line of regression and thus there does not appear to be very strong correlation between internet use rate and number of suicides.

On running the GLM procedure for the model, it was found that the p-value (=0.3634) of association is not significant at 5% level of significance. Hence, we can conclude that there isn’t a significant association between internet use and suicide.

The model predicted is Y = 9.1223 + 0.0157 * X, where Y is number of suicides (response) and X is percentage of internet users (predictor). There is indication of very slight positive association between the response and predictor. This association however, is not significant as can be inferred from the very high p-value, i.e., 0.3634. 

Also, the R square coefficient is very low 0.004722, i.e. the model explains only 0.5% of the variability if the response variable.

Part 2

Centering the explanatory variable does not bring about any changes to the overall nature/fit statistics of the model although there is some change in the intercept coefficient based on the centered explanatory variable. The predicted model is: Y = 9.1223 + 0.01567 * X, where Y is response and X is the usual predictor centered around its mean.

Part 3

Here again, it is seen that the p-value of association (0.7984) is not significant at 5% level of significance. The model predicted is Y = 9.7575 - 0.2444 * X where Y is number of suicides and X indicates the level of internet use, i.e, low or high. The model has an R square coefficient of 0.000374 which is very low.

Since the predictor is a binary categorical variable, the model above basically gives the mean number of suicides for the two categories of the predictor variable, viz, low (= 0) and high (= 1) levels of internet use, as can be verified from the table of means shown below,

[reetikasarkar-blog]

Machine Learning for Data Analysis (Week 2: Running a Random Forest)

Random forests are predictive models that allow for a data driven exploration of many explanatory variables in predicting a response or target variable.

Random forest analysis was performed to evaluate the importance of a series of explanatory variables in predicting a binary, categorical response variable. The following explanatory variables were included as possible contributors to a random forest evaluating suicide levels (response variable) based on internet use, income, employment rate, urbanisation, alcohol consumption, CO2 emissions, oil usage, polity scores and life expectancy. 

The SAS code I used is as under:

LIBNAME mydata "/courses/d1406ae5ba27fe300 " access=readonly; DATA new; set mydata.gapminder;

/* Setting labels for each of the variables */ LABEL suicidegroup = "Level of suicide";  

/* Data management */ IF suicideper100th EQ . THEN suicidegroup = .; ELSE IF suicideper100th LE 8 THEN suicidegroup = '1'; ELSE suicidegroup = '2';

IF suicidegroup NE .; IF internetuserate NE .; IF incomeperperson NE .; IF urbanrate NE .; IF alcconsumption NE .; IF co2emissions NE .; IF employrate NE .; IF oilperperson NE .; IF polityscore NE .; IF lifeexpectancy NE .;

PROC SORT; BY country;

PROC HPFOREST; TARGET suicidegroup/level=nominal; INPUT internetuserate incomeperperson urbanrate alcconsumption co2emissions employrate oilperperson polityscore lifeexpectancy/level=interval;

RUN;

The output generated by the HPFOREST procedure is shown below.

The average square error of the model is 0.222 and the misclassification rate is 0.333, i.e., about 67% of observations are correctly predicted.

The Fit Statistics for the first and the last ten decision trees are shown below.

It can be seen that as number of trees increases, the average square error and misclassification rate of both training observations and the out-of-bag (OOB) observations decreases. The decrease in misclassification rate among the OOB observations however, is not uniform in the sense that lower levels of misclassification are observed with 7 trees (0.333) while for 100 trees it stands at 0.300 which is not a very significant reduction.

From the loss reduction variable importance table above, the explanatory variables with the highest relative importance scores were internetuserate, incomeperperson, urbanrate and alcconsumption, indicating that these factors contribute significantly to good prediction of the model.

[reetikasarkar-blog]

Machine Learning for Data Analysis (Week 1: Making a Decision Tree)

The aim of this study is to draw decision trees for determination of factors affecting suicide rates among countries. The data used for this study has been obtained from the GapMinder dataset. 

Decision trees attempt to find a strong relationship between input values and target values in a group of observations that form a data set. Decision tree analysis was performed to test nonlinear relationships among a series of explanatory variables, namely, internetuserate, incomeperperperson, urbanrate, alcconsumption, co2emissions, employrate, oilperperson, polityscore, lifeexpectancy and a binary, categorical response variable - suicideper100th (number of suicides per 100,000 people). All possible separations (categorical) have been tested. For the present analyses, the entropy “goodness of split” criterion was used to grow the tree and a cost complexity algorithm was used for pruning the full tree into a final subtree.

The following explanatory variables were included as possible contributors to a classification tree model evaluating suicide rate/group (my response variable), internet use rate, income level, level of urbanisation, level of alcohol consumption, CO2 emission level, employment rates, oil consumption levels, polity scores and the life expectancy level.

The SAS Code I used is as under:

LIBNAME mydata "/courses/d1406ae5ba27fe300 " access=readonly; DATA new; set mydata.gapminder;

/* Setting labels for each of the variables */ LABEL suicidegroup = "Level of suicide" internetgroup = "Level of internet use" incomegroup = "Income per person" urbangroup = "Level of urbanisation" alcgroup = "Alcohol consumption" co2group = "CO2 Emission level" employgroup = "Employment level" oilgroup = "Oil Consumption level" politygroup = "Polity score" lifegroup = "Level of life expectancy";  

/* Data management */ IF suicideper100th EQ . THEN suicidegroup = .; ELSE IF suicideper100th LE 8 THEN suicidegroup = '1'; ELSE suicidegroup = '2';

IF internetuserate EQ . THEN internetgroup = .; ELSE IF internetuserate LE 35 THEN internetgroup = '1'; ELSE IF internetuserate LE 60 THEN internetgroup = '2'; ELSE internetgroup = '3';

IF incomeperperson EQ . THEN incomeperperson = .; ELSE IF incomeperperson LE 1000 THEN incomegroup = '1'; ELSE IF incomeperperson LE 15000 THEN incomegroup = '2'; ELSE incomegroup = '3';

IF urbanrate EQ . THEN urbangroup = .; ELSE IF urbanrate LE 50 THEN urbangroup = '1'; ELSE urbangroup = '2';

IF alcconsumption EQ . THEN alcgroup = .; ELSE IF alcconsumption LE 5 THEN alcgroup = 1; ELSE alcgroup = 2;

IF co2emissions EQ . THEN co2group = .; ELSE IF co2emissions LE 200000000 THEN co2group = 1; ELSE co2group = 2;

IF employrate EQ . THEN employgroup = .; ELSE IF employrate LE 50 THEN employgroup = 1; ELSE IF employrate LE 65 THEN employgroup = 2; ELSE employgroup = 3;

IF oilperperson EQ . THEN oilgroup = .; ELSE IF oilperperson LE 1 THEN oilgroup = 1; ELSE oilgroup = 2;

IF polityscore EQ . THEN politygroup = .; ELSE IF polityscore LE 0 THEN politygroup = 1; ELSE politygroup = 2;

IF lifeexpectancy EQ . THEN lifegroup = .; ELSE IF lifeexpectancy LE 65 THEN lifegroup = 1; ELSE IF lifeexpectancy LE 75 THEN lifegroup = 2; ELSE lifegroup = 3;

IF suicidegroup NE .; IF internetgroup NE .; IF incomegroup NE .; IF urbangroup NE .; IF alcgroup NE .; IF co2group NE .; IF employgroup NE .; IF oilgroup NE .; IF politygroup NE .; IF lifegroup NE .;

/* Obtaining the summary of each of the variables */ PROC MEANS data=mydata.gapminder mean min q1 median q3 max; VAR internetuserate suicideper100th incomeperperson urbanrate alcconsumption co2emissions employrate oilperperson polityscore lifeexpectancy;

PROC SORT; BY country;

ODS GRAPHICS ON;

PROC HPSPLIT seed=554321; CLASS suicidegroup internetgroup incomegroup urbangroup alcgroup co2group employgroup oilgroup politygroup lifegroup; MODEL suicidegroup = internetgroup incomegroup urbangroup alcgroup co2group employgroup oilgroup politygroup lifegroup; GROW ENTROPY; PRUNE COSTCOMPLEXITY;

RUN;

The decision tree generated by SAS is as under: 

From the initial node, it is seen that almost 67% of the countries have suicide rates higher than 8 per 100 thousand persons. The first subdivision in the suicide rates has been made on the basis of internet use rate. Higher rates of suicide are more pronounced (84.62%) in group 3, i.e., countries with high internet use rate, than those with low and moderate use wherein almost 53% of the countries have higher suicide rates.

From the confusion matrix above, it is seen that the model error rates are 0.1500 and 0.0750, while the corresponding values of Sensitivity is 0.8500 (85% of type 1 values/low suicide rates are correctly predicted) and Specificity is 0.9250 (92.5% of type 2 values/high suicide rates) are correctly predicted. 

The area under the ROC (Receiver’s Operating Characteristics) curve for suicide level is 0.9438 which indicates high accuracy in correctly predicting higher levels of suicide.

The table shows that the variable employgroup has the largest value. This means that “Employment Rate” is the most important consideration in determining the usability of the sample. 

[reetikasarkar-blog]

Week 1: Introduction to Regression

Data Management and Description

Under this study, I am looking at studying the degree of association (if any) between internet usage and the rate of suicides. Further, the level of income may also have an effect on this probable association and also on each of these variables. 

Here, the explanatory variable is internet penetration which indicates the number of people (per 100 persons) with access to the worldwide web through computers, mobiles and other such devices. The responses from the data set were scaled from very low levels to very high levels. Countries with less than 20 users were taken to have very low usage while those with higher than 80 were grouped as very high internet penetration countries.

The response variable is the suicide rate which indicates the number of suicides per 100,000 persons. In this case, the values were classified into three classes, low, moderate and high. Values below 6 were low and values above 12 were considered high. 

The income level (moderator) indicates the Gross Domestic Product per capita by Purchasing Power Parities. It was classified into five categories ranging from very low to very high. Countries with income below 1000 were very low income countries and those with incomes above 25000 had very high levels of income.

The data management is summarised below:

[reetikasarkar-blog]

Week 1: Introduction to Regression

Data Collection

The purpose of data collection was to promote sustainable global development and achievement of the United Nations Millennium Development Goals by increased use and understanding of statistics and other information about social, economic and environmental development at local, national and global levels. Data is collected through observational studies. There is normally one documentation for each indicator in the GapMinder data sets, but sometimes several indicators are described in one documentation. The data set consists of time series based statistical data collected from several sources such as the World Bank, WHO, etc.

Data on internet penetration (Internet users per 100 people) was collected annually from the year 1990 to 2014 based on the database maintained by International Telecommunication Union, World Telecommunication/ICT Development Report and database, and World Bank estimates. Data was aggregated using weighted averages. 

Data on suicide rates indicates mortality due to self-inflicted injury, per 100 000 standard population (age adjusted). It was obtained from WHO using combination of time series from WHO Violence and Injury Prevention (VIP) and data from WHO Global Burden of Disease 2002 and 2004. Data for the year 2005 was extrapolated for a majority of countries. This was done in order to enable comparisons with other indicators which only have data for 2005. 

Data on income rates was compiled by GapMinder using data from several sources, such as official international statistics, various historical sources and their own estimates. It includes the Gross Domestic Product per capita by Purchasing Power Parities (in international dollars, fixed 2005 prices). The inflation and differences in  the cost of living between countries has been taken into account. The documentation of the indicator may be seen in detail here.

[reetikasarkar-blog]

Week 1: Introduction to Regression

The aim of this week is to choose a suitable dataset for study and to write about the data involving description of the sample, data collection procedure and data management.

Sample

The sample is from the datasets collated and maintained by GapMinder, which is a Swiss non-profit venture promoting sustainable global development and achievement of the United Nations Millennium Development Goals by increased use and understanding of statistics and other information about social, economic and environmental development at local, national and global levels. They make use of time series based statistical data collated from several sources such as the World Bank, WHO, etc., to provide data on several indicators across all the UN member countries of the world. Since the data on each indicator are data combined from more than one source, the analysis was conducted at an aggregate level. The number of observations in the data set are varying across the indicators as for some of the indicators, data is not available for certain countries. 

From the entire set of indicators available in the GapMinder data set, the variables of interest are internet penetration and suicide rates, along with a study of the effect of income on both these variables. The unique identifier for the GapMinder dataset is name of the country.

Internet penetration is shown by a variable namely, ‘internetuserate’ with N=192.

Suicide rate is indicated by a variable called ‘suicideper100th’ with N=191.

Income levels are indicated by ‘incomeperperson’ with N=190.

[reetikasarkar-blog]

Regression Modelling in Practice

This is the third module in the Data Analysis and Interpretation Specialization course of Coursera offered by Wesleyan University. The module is being covered over a period of four weeks involving lessons on Introduction to Regression, Basics of Linear Regression, Multiple Regression and Logistic Regression.

[reetikasarkar-blog]

Week 4: Testing moderation in context of correlation

This week the aim is to test the effect of a moderator on the association between two variables of interest. The dataset chosen is GapMinder and the variables of interest are “internetuserate” and “suicideper100th”. The potential moderator for these variables is “incomeperperson”.

The necessary tests to ascertain the existence of association in the presence of moderators has been performed using both SAS and Python, as shown below.

1) SAS

Program Code

Output

The output generated by the above program may be seen here.

2) Python 2.7

Program Code

# Import necessary libraries import pandas as pd import numpy import scipy.stats import seaborn import matplotlib.pyplot as plt

# Import the entire dataset to memory data = pd.read_csv("gapminder.csv", low_memory=False)

# Set PANDAS to show all columns in DataFrame pd.set_option('display.max_columns', None)

# Set PANDAS to show all rows in DataFrame pd.set_option('display.max_rows', None)

# Upper-case all Dataframe column names data.columns = map(str.upper, data.columns)

# Bug fix for display formats to avoid run-time errors pd.set_option('display.float_format', lambda x: '%f' %x)

# Ensure each of the columns are numeric data['INTERNETUSERATE'] = pd.to_numeric(data['INTERNETUSERATE'], errors='coerce') data['SUICIDEPER100TH'] = pd.to_numeric(data['SUICIDEPER100TH'], errors='coerce') data['INCOMEPERPERSON'] = pd.to_numeric(data['INCOMEPERPERSON'], errors='coerce')

data['INTERNETUSERATE'] = data['INTERNETUSERATE'].replace(' ', numpy.nan) data['SUICIDEPER100TH'] = data['SUICIDEPER100TH'].replace(' ', numpy.nan) data['INCOMEPERPERSON'] = data['INCOMEPERPERSON'].replace(' ', numpy.nan)

data_clean = data.dropna()

print scipy.stats.pearsonr(data_clean['SUICIDEPER100TH'], data_clean['INTERNETUSERATE']) print '       '

def INCOMEGROUP (row):   if row['INCOMEPERPERSON'] <= 1000 : return 1   elif row['INCOMEPERPERSON'] <= 3000 : return 2   elif row['INCOMEPERPERSON'] <= 10000 : return 3   elif row['INCOMEPERPERSON'] <= 25000 : return 4   elif row['INCOMEPERPERSON'] > 25000 : return 5

data_clean['INCOMEGROUP'] = data_clean.apply (lambda row: INCOMEGROUP (row), axis=1)

chk1 = data_clean['INCOMEGROUP'].value_counts(sort=False, dropna=False) print chk1 print '       '

sub1 = data_clean[(data_clean['INCOMEGROUP']== 1)] sub2 = data_clean[(data_clean['INCOMEGROUP']== 2)] sub3 = data_clean[(data_clean['INCOMEGROUP']== 3)] sub4 = data_clean[(data_clean['INCOMEGROUP']== 4)] sub5 = data_clean[(data_clean['INCOMEGROUP']== 5)]

print 'association between suiciderate and internetuserate for VERY LOW income countries' print scipy.stats.pearsonr(sub1['SUICIDEPER100TH'], sub1['INTERNETUSERATE']) print '       ' print 'association between urbanrate and internetuserate for LOW income countries' print scipy.stats.pearsonr(sub2['SUICIDEPER100TH'], sub2['INTERNETUSERATE']) print '       ' print 'association between urbanrate and internetuserate for MODERATE income countries' print scipy.stats.pearsonr(sub3['SUICIDEPER100TH'], sub3['INTERNETUSERATE']) print '       ' print 'association between urbanrate and internetuserate for HIGH income countries' print scipy.stats.pearsonr(sub4['SUICIDEPER100TH'], sub4['INTERNETUSERATE']) print '       ' print 'association between urbanrate and internetuserate for VERY HIGH income countries' print scipy.stats.pearsonr(sub5['SUICIDEPER100TH'], sub5['INTERNETUSERATE']) print '       '

# Scatter Plots scat1 = seaborn.regplot(x="INTERNETUSERATE", y="SUICIDEPER100TH", data=sub1, label="Very Low Income") scat2 = seaborn.regplot(x="INTERNETUSERATE", y="SUICIDEPER100TH", data=sub2, label="Low Income") scat3 = seaborn.regplot(x="INTERNETUSERATE", y="SUICIDEPER100TH", data=sub3, label="Moderate Income") scat4 = seaborn.regplot(x="INTERNETUSERATE", y="SUICIDEPER100TH", data=sub4, label="High Income") scat5 = seaborn.regplot(x="INTERNETUSERATE", y="SUICIDEPER100TH", data=sub5, label="Very High Income") plt.legend(loc='best', frameon=False) plt.xlabel('Internet Penetration') plt.ylabel('Suicide Rate') plt.title('Scatterplot for the Association Between Internet Penetration and Suicide Rate across income groups') plt.show() print scat1 print scat2 print scat3 print scat4 print scat5

Output

(0.068005290321515185, 0.37259750532698521)

1    51 2    46 3    38 4    21 5    18 Name: INCOMEGROUP, dtype: int64

association between suiciderate and internetuserate for VERY LOW income countries (0.065344298510647297, 0.64869906224206764)

association between urbanrate and internetuserate for LOW income countries (0.12813601111511572, 0.39607611429593026)

association between urbanrate and internetuserate for MODERATE income countries (0.33068627104192383, 0.04257658401101342)

association between urbanrate and internetuserate for HIGH income countries (0.36016539267459113, 0.10876213803045545)

association between urbanrate and internetuserate for VERY HIGH income countries (-0.0066689924392987644, 0.97904760099771482)

C:/Users/Reetika/Desktop/Research Project_Coursera/Week3_PythonProgram.py:44: SettingWithCopyWarning: A value is trying to be set on a copy of a slice from a DataFrame. Try using .loc[row_indexer,col_indexer] = value instead

See the caveats in the documentation: http://pandas.pydata.org/pandas-docs/stable/indexing.html#indexing-view-versus-copy  data_clean['INCOMEGROUP'] = data_clean.apply (lambda row: INCOMEGROUP (row), axis=1)

Axes(0.125,0.125;0.775x0.775) Axes(0.125,0.125;0.775x0.775) Axes(0.125,0.125;0.775x0.775) Axes(0.125,0.125;0.775x0.775) Axes(0.125,0.125;0.775x0.775)

Interpretations:

1) Income level 1 (very low): The correlation coefficient between internet penetration and suicide rate is very low indicating that there is no linear relationship between the two. This is further ascertained by the non-significant p-value.

2) Income level 2 (low):  The correlation coefficient between internet penetration and suicide rate is low indicating that there is no linear relationship between the two. This is further ascertained by the non-significant p-value.

3) Income level 3 (moderate):  Although the correlation coefficient between internet penetration and suicide rate is low, the resultant p-value is significant which indicates that for moderate incomes, there may be an association between level of internet penetration and suicide rate.

4) Income level 4 (high): The correlation coefficient between internet penetration and suicide rate is low indicating that there is no linear relationship between the two. This is further ascertained by the non-significant p-value.

5) Income level 5 (very high):  The correlation coefficient between internet penetration and suicide rate is negative and very close to 0 indicating that there is a very slight negative (almost nil) linear association between the two. This is further ascertained by the non-significant p-value.

[reetikasarkar-blog]

Week 3: Generating Pearson Correlation Coefficient

In the third week of the Data Analysis Tools course, the aim is to compute Pearson Correlation coefficient between two quantitative variables. Under the research question chosen, the variables of interest are internetuserate, suicideper100th and incomeperperson from the GapMinder dataset.

Like previous weeks, I will compute the Pearson Correlation Coefficient using both SAS and Python.

1) Using SAS

Program Code

 Output

Model Interpretation for Pearson Correlation Test:

The Pearson correlation coefficient for incomeperperson and internetuserate is 0.75094 indicating that there is a strong positive linear association between income and internet use rate. The p-value for the correlation test is significant (<0.0001). Hence, there is sufficient evidence  to conclude that there is a linear relationship in the population between income level and internet use rate.

The Pearson correlation coefficient for incomeperperson and suicideper100th is 0.00656 which is very close to 0, indicating that suicide rates may not be linearly related to income levels. Also, the p-value for the correlation test is not significant (>0.05). Hence, there is not enough evidence to conclude that there is a linear relationship in the population between income level and suicide rate.

2) Using Python 2.7

Program Code

Output

Model Interpretation for Pearson Correlation Test:

The Pearson correlation coefficient for incomeperperson and internetuserate is 0.80533 indicating that there is a strong positive linear association between income and internet use rate. The p-value for the correlation test is significant (<0.05). Hence, there is sufficient evidence  to conclude that there is a linear relationship in the population between income level and internet use rate.

The Pearson correlation coefficient for incomeperperson and suicideper100th is 0.00327 which is very close to 0, indicating that suicide rates may not be linearly related to income levels. Also, the p-value for the correlation test is not significant (>0.05). Hence, there is not enough evidence to conclude that there is a linear relationship in the population between income level and suicide rate.

[reetikasarkar-blog]

Week 2: Chi-square Test of Independence

This week, the aim is to run a chi-square test of independence for categorical explanatory variable and categorical response variable. The variables chosen are internetuserate and suicideper100th from the ‘GapMinder’ dataset. The hypotheses of interest are:

Null Hypothesis: All suicide rates are similar across all levels of internet penetration.

Alternate Hypothesis: Not all suicide rates are similar across all levels of internet penetration.

Since the variables of interest are all quantitative, they need to be converted into categorical data, after which the chi-square test of independence may be applied to them. I will be carrying out the tests for independence using both SAS and Python.

1) Using SAS

Program Code

Output

The output may be seen here.

2) Using Python 2.7

Program Code

Output

The output may be seen here.

Model Interpretation for Chi-Square Test:

Since the p-value obtained is lower than 0.05 (significant), we reject the null hypothesis and conclude that not all suicide rates (response variable) are similar across all levels of internet penetration (explanatory variable) or equivalently, there is a relationship between levels of internet penetration and suicide rates. Thus, we need to find which pairs of internet penetration levels are associated.

Model Interpretation for post hoc Chi-Square Test results:

In order to minimise the overall Type-I error rate in doing a pair-wise comparison of the 10 pairs of internet penetration levels, we compute the Bonferroni adjustment, i.e., 0.05/10=0.005.

On carrying out the pair-wise post hoc comparisons, it is found that the only pairs with significant p-value is (1,2), thus indicating that the significance is due to very low (category 1) and low (category 2) levels of internet use rate.

NOTE: It may also be noted that a Chi-square test of independence produces invalid results in case any of the cell counts is lesser than 5, and in such cases, a suitable correction needs to be applied (which is beyond the ambit of this week’s course).

[reetikasarkar-blog]

Week 1 : Hypothesis Testing and ANOVA

This is the first week of the second course from the Data Analysis and Interpretation Specialization. In the first course, my research question was to ascertain the existence of any association between Internet use rate and Suicide rate. I used the GapMinder dataset for this purpose and data visualisation for the chosen variables have also been recorded on this blog. Moving on, the aim in this course is to an in-depth analysis of the association of the two variables considering the effect of income rate and urban rate on both these variables. 

This week, I will carry out ANOVA tests and post hoc tests (if required) to study whether there is association between:

a) income level (explanatory) and internet use rate (response)

Null hypothesis: There is no association between income level and internet use rate

Alternate hypothesis:  There is an association between income level and internet use rate

b) level of urbanisation (explanatory) and internet use rate (response)

Null hypothesis: There is no association between urbanisation and internet use rate

Alternate hypothesis:  There is an association between urbanisation and internet use rate

c) income level (explanatory) and suicide rate (response)

Null hypothesis: There is no association between income level and suicide rate  

Alternate hypothesis: There is an association between income level and suicide rate 

d)  level of urbanisation (explanatory) and suicide rate (response)

Null hypothesis: There is no association between urbanisation and suicide rate

Alternate hypothesis:  There is an association between urbanisation and suicide rate

For this purpose, I would first group the explanatory variables, viz., income rates and urban rates from the dataset into five categories and then perform ANOVA on each of them for the variables ‘internetuserate’ and ‘suicideper100th’.

I would be doing the exercise using both SAS and Python.

1) Using SAS

Program Code

Output

The output may be seen here.

2) Using Python 2.7

Program Code

Output

The output may be seen here.

Interpretations

a) income and internet use rate

Model Interpretation for ANOVA: The value of F statistic is highly significant, thus indicating strong association between income and internet use rate. Since the mean internet use increases with increase in income, we further conclude that there is strong positive association.

Model Interpretation for post hoc ANOVA results: The Tukey HSD test reveals that there is significant association between income and internet use rate at all levels of income.

b) urbanisation and internet use rate

Model Interpretation for ANOVA: The value of F statistic is significant, thus indicating strong association between urbanisation and internet use rate. Also since the mean internet use increases with increase in urbanisation, we may further conclude that there is positive association.

Model Interpretation for post hoc ANOVA results: The Tukey HSD test reveals that the pairs of groups with significant association are (1,4), (1,5), (2,4), (2,5), (3,4), (3,5) and (4,5).

c) income and suicide rate

Model Interpretation for ANOVA: Since the F statistic is not significant, we do not have sufficient evidence to reject the null hypothesis. Thus, we may conclude that any indication of an effect of income levels on suicide rates cannot be attributed to an association between them. They may have crept in by chance. This is further substantiated by the mean suicide rates for different levels of income which are not distributed uniformly. Since we fail to reject null hypothesis, the need for post hoc tests does not arise.

d) urbanisation and suicide rate

Model Interpretation for ANOVA:  Since the F statistic is not significant, we do not have sufficient evidence to reject the null hypothesis. Thus, we may conclude that any indication of an effect of urbanisation on suicide rates cannot be attributed to an association between them. They may have crept in by chance. Since we fail to reject null hypothesis, the need for post hoc tests does not arise.

[reetikasarkar-blog]

Week 4 : Visualising Data

This week, the goal is to graphically represent the managed data in order to understand it better. Moving ahead, I have also included another variable from the GapMinder dataset, namely, ‘urbanrate’ as I would also like to examine whether internet usage and suicide rates vary with the degree of urbanisation, in addition to their association with the income levels and whether or not there is any association between internet usage rates and suicide rates. The codebook for this analysis is given below:

Like previous weeks, this week too, I am going to do the graphical representation of the data using both SAS and Python.

1) Using SAS

Program Code

Output

The output may be seen here.

Observations

From the scatter plot between internet use rate and suicide rate, it cannot be conclusively stated whether the two variables are associated or not as the points are widely scattered and no definite pattern can be observed. This is substantiated by the bar chart which does not depict any particular pattern between the two variables. However, it is seen that suicide rates are comparatively higher for high and very high income groups.

From the univariate bar charts, it is observed that most of the countries have very low number of internet users. As regards the suicide rates, majority of the countries show very low to low rates. Most of the countries from the dataset belong to very low income groups while majority of the countries recorded comparatively higher rates of urbanisation.

From the bivariate bar charts, it is seen that higher incomes and higher rate of urbanisation indicate higher internet use rate. Additionally, suicide rates decrease with increase in income but they are highest for very high income groups. In relation to urban rates, suicide rate are highest for very low level of urbanisation and lowest for very high urban population, although there does not appear to exist a specific downward trend.

2) Using Python 2.7

Program Code

Output

The output may be seen here.

Observations

From the scatter plot between internet use rate and suicide rate, it cannot be conclusively stated whether the two variables are associated or not as the points are widely scattered and no definite pattern can be observed. 

From the univariate bar charts, it is observed that most of the countries have very low number of internet users (below the 20th percentile). As regards suicide rates, majority of the countries are below the 40th percentile and there is a downward trend in the data. Most of the countries from the dataset belong to very low income groups (below 20th percentile) while majority of the countries were between 40th and 80th percentile.

From the bivariate bar charts, it is seen that higher incomes and higher rate of urbanisation indicate higher internet use rate. Additionally, suicide rates increase with increase in income and then decrease. In relation to urban rates, suicide rate are highest for lower level of urbanisation and there does not appear to exist a specific trend in the plot.

[reetikasarkar-blog]

Week 3 : Making Data Management Decisions

Since the dataset I have chosen is based on data from countries across the world and does not have ‘no response’ categories, data management in terms of excluding such values is not pertinent. Also, the value of the variables being quantitative and not categorical, reversing code values might not be very helpful for the purpose of this study. 

I generated the frequency distributions for the three variables, viz., internetuserate, suiciderate and incomeperperson using both SAS as well as Python 2.7 in order to optimize my learning. Below are the screenshots of the program code and the output for each of the two methods. 

1) Using SAS

The categories defined for the purpose of studying the association among the chosen are as follows:

SAS Program Code

Output

Observations

I have collapsed the three variables into three new variables namely, internetpenetration, suiciderate and incomegroup to form a class based frequency distribution. 

From the frequency table for internet users, it may be seen that 39.58% of the countries have very low internet penetration and 8.85% have very high internet penetration. Majority of the countries (75.52%) have very low to moderate internet penetration.

As regards the suicide rates, majority of the countries (72.25%) have very low to low rates while only 4.19% of the countries recorded very high suicide rates.

From the frequency table for income groups, it can be seen that most of the countries lie in the very low to low income groups while almost equal number of countries have high and very high incomes.

2) Using Python 2.7

I have divided the variables into groups using the pd.cut function from the pandas library and generated 5 classes based on 20th, 40th, 60th and 80th percentiles.

Python Program Code

Output

Observations

The describe function returns the summary statistics related to the variables under study. However, on running the program the quantiles would show up as nan which is something I have bee unable to figure out so far.

The three variables under study have been recoded as INTERNETPENETRATION, SUICIDERATE and INCOMEGROUP to form a frequency distribution.

From the frequency table for internet users, it may be seen that majority (73 or 38%) of the countries are below the first quartile.

As regards the suicide rates, majority of the countries (157 or approx. 82.2%) are below the 40th percentile.

From the frequency table for income groups, it can be seen that most of the countries lie below the 20th percentile (162 or 85.26%).

Conclusion

While running the SAS program, it was found that the summary statistics for the variables under study are displayed without any glitches using the proc means statement. Subsequently, the frequency distribution was generated using the IF and ELSE statements in order to group the data into five classes ranging from very low to very high.

Analogously, on attempting to display the same using Python, it was found that the quartiles showed up as nan values, although the class intervals based on quartiles displayed consistent results.