• Last week: ANOVA and its associated nightmares.
• Today all correlation and regression.

A correlation is mostly the first step. You want to find out the strength of an association

Sample correlation coefficient. These are the formulae.

\(r = \frac{\Sigma(x_i - \bar{x})(y_i - \bar{y})}{\sqrt{\Sigma(x_i - \bar{x})^2\Sigma(y_i - \bar{y})^2}}\)

is equivalent to

\(r = \frac{cov (x,y)}{\sqrt{var(x)*var(y)}}\)

Covariance divided by the geometric mean of variances.

Read back your notes or check here

Correlation is not causation (but some work in computing science, can infer this in basic associations!)

Check this

# custom confidence interval 90%
cor.test(iris$Sepal.Length, iris$Sepal.Width, conf.level = .90)

## 
##  Pearson's product-moment correlation
## 
## data:  iris$Sepal.Length and iris$Sepal.Width
## t = -1.4403, df = 148, p-value = 0.1519
## alternative hypothesis: true correlation is not equal to 0
## 90 percent confidence interval:
##  -0.24846981  0.01754741
## sample estimates:
##        cor 
## -0.1175698

The Pearson correlation between sepal length and sepal width was not significant (r(148)= -.12, p=.15).

For further info. see the tutorial referenced at the end by Willems.

??corrplot

This is based on Fox (2015) and Fox (2002)

library(carData)
prestige<-carData::Prestige

education: The average number of years of education for occupational incumbents.

income: The average income of occupational incumbents, in dollars.

women: The percentage of women in the occupation.

prestige:The average prestige rating for the occupation.

census: The code of the occupation used in the survey.

type: Professional and managerial(prof), white collar(wc), blue collar(bc), or missing(NA)

Load the Prestige dataset

Make a correlation matrix with all the continuous variables. (use filter from dplyr to remove strings)

And: either export a .docx or make a correlation plot.

When would you use linear regression?

model_education<-lm(prestige~education, data=prestige)
summary(model_education)

## 
## Call:
## lm(formula = prestige ~ education, data = prestige)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -26.0397  -6.5228   0.6611   6.7430  18.1636 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -10.732      3.677  -2.919  0.00434 ** 
## education      5.361      0.332  16.148  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 9.103 on 100 degrees of freedom
## Multiple R-squared:  0.7228, Adjusted R-squared:   0.72 
## F-statistic: 260.8 on 1 and 100 DF,  p-value: < 2.2e-16

plot_edu

plot(model_education, which = 2)

plot(model_education, which = 3)

plot(model_education, which = 4)

Perhaps for basic models but it will quickly become complex.

Some of these plots give you information which you could not get from just scatter plots.

require(lmtest)
bptest(model_education)

## 
##  studentized Breusch-Pagan test
## 
## data:  model_education
## BP = 0.70577, df = 1, p-value = 0.4009

dwtest(model_education) #Durbin-Watson test.

## 
##  Durbin-Watson test
## 
## data:  model_education
## DW = 1.4392, p-value = 0.001598
## alternative hypothesis: true autocorrelation is greater than 0

Also note that there is one other alternative typically used in econometrics (Breusch-Godfroy). Check here.

Box.test(residuals(model_education), type="Ljung-Box") #Box test

## 
##  Box-Ljung test
## 
## data:  residuals(model_education)
## X-squared = 7.9574, df = 1, p-value = 0.004789

Visual inspection suggested that a linear fit was appropriate. Visual inspection also suggested that the errors from the OLS regression model were approximately normally distributed. A Breusch-Pagan test showed that the assumption of homoscedasticity could be upheld (\(\chi^2\)(1)= 0.71, p=.40). However, a Durbin-Watson test indicated that the assumption of independence of errors was violated (Durbin-Watson test= 1.44, p=.002).

Table X describes a hierarchical regression on prestige. The first model showed a statistically significant, positive, effect of years of education on prestige. In the second model, the effect of education was upheld, while controlling for the % of women. The effect of % women was also statistically significant.

Note be wary with interaction terms! (center them beforehand - values might not be correct from QuantPsyc). Alternative function: lm.beta.

require(QuantPsyc)
lm.beta(model_education)

## education 
## 0.8501769

# view results
results_coeff_B

## 
## ORDINARY NONPARAMETRIC BOOTSTRAP
## 
## 
## Call:
## boot(data = prestige, statistic = coeff_B, R = 10000)
## 
## 
## Bootstrap Statistics :
##     original      bias    std. error
## t1* 5.427796 0.003518518    0.288098

log1<-glm(diabetes~mass, data=pima_data, family=binomial)
log2<-glm(diabetes~mass+age, data=pima_data, family=binomial)
summary(log1)

## 
## Call:
## glm(formula = diabetes ~ mass, family = binomial, data = pima_data)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.0094  -0.9184  -0.6598   1.2254   1.9107  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -3.99682    0.42885   -9.32  < 2e-16 ***
## mass         0.10250    0.01261    8.13 4.31e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 981.53  on 756  degrees of freedom
## Residual deviance: 904.89  on 755  degrees of freedom
##   (11 observations deleted due to missingness)
## AIC: 908.89
## 
## Number of Fisher Scoring iterations: 4

summary(log2)

## 
## Call:
## glm(formula = diabetes ~ mass + age, family = binomial, data = pima_data)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.0306  -0.8907  -0.5692   1.1011   2.1555  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -5.740908   0.537331 -10.684  < 2e-16 ***
## mass         0.107810   0.013071   8.248  < 2e-16 ***
## age          0.045746   0.007017   6.519 7.07e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 981.53  on 756  degrees of freedom
## Residual deviance: 860.61  on 754  degrees of freedom
##   (11 observations deleted due to missingness)
## AIC: 866.61
## 
## Number of Fisher Scoring iterations: 4

• Load last week's Salaries dataset. If necessary, calculate monthly_salary again.

• calculate a Spearman and Pearson correlation between years since Ph.D. and monthly salary. Conduct a significance test for one of them. Pick one and discuss the outcome

• Build a hierarchical regression model, in step 1 use monthly salary as dependent and year since Ph.D. as independent, in step 2 add gender (keep all variables from step 1), in step 3 add rank (keep all variables from step 2). Export the hierarchical regression table.

• Test the assumptions for the final model, how would you combat any issues that arise in terms of violations of these assumptions?

• Calculate the standardized coefficients for model 1 (\(\beta\)). Interpret this coefficient. Bootstrap this coefficient and report as you would do in a paper.

• What does the final model predict in terms of monthly salary shift, when there is a shift of being 3 years post-Ph.D.

• Install and load the 'caret' package. load the GermanCredit dataset. ??GermanCredit to find out more also here. It contains data on credit scores for a 1,000 cases and has 62 variables. As done before, you might want to explore the descriptive statistics.

• Make a corrplot() for up to 10 continuous variables in the GermanCredit dataset. Export a correlation matrix (up to 10 variables).

• Conduct a hierarchical logistic regression on 'Class' as dependent, in each step keeping the previous variables. In step 1 add Amount. In step 2 add CreditHistory.PaidDuly. In step 3 add Age. In step 4 add Telephone. Export a table with your key results.

• Choose some categorical variables and make a (beautiful) mosaic plot with credit class as dependent.

• Calculate all odds ratios for model 4. Calculate the odds ratio for a shift in 1,000 DM (Deutsche Mark - pre-euro era) in credit status for model 4.

• BONUS: Make a mosaic plot via ggmosaic() with credit class as dependent.

• BONUS: Have a look at ROC curves and plot one for the credit data.

Also check the reading list! (many more than listed here)

• Crawley, M. J. (2013). The R book: Second edition. New York, NY: John Wiley & Sons.
• Fox, J. (2015). Applied regression analysis and generalized linear models. London, UK: Sage.
• Torres-Reyna, O. (2010). Regression 101.
• Wickham, H., & Grolemund, G. (2017). R for data science. Sebastopol, CA: O’Reilly.
• Willems, K. (2017). R correlation tutorial
• Zuur, A., Ieno, E.N., Walker, N., Saveliev, A.A., & Smith, G.M. (2009). Mixed effects models and extensions in ecology with R. New York, NY: Springer.