---
title: Regression
author: Blanca Calvo
date: "Generation date: `r format(Sys.time(), '%b %d, %Y - %H:%M:%S')`"
output:
  html_document:
    toc: true
    code_folding: show
    toc_float:
        collapsed: false
        smooth_scroll: true
    number_sections: true
---

```{r setup, message=F, echo=F}
# some custom options to generate a nice html file
options(digits = 3)
options(width=100)
packages <- c("knitr","pander")
if (length(setdiff(packages, rownames(installed.packages()))) > 0) {
  install.packages(setdiff(packages, rownames(installed.packages())), repos='http://cran.us.r-project.org')
}


library(knitr)
opts_chunk$set(cache=F, comment='#', tidy = T, tidy.opts=list(blank=T, width.cutoff=100), fig.align='center', message=F, warning=T, fig.width=6, fig.height=6, dev='svg',eval=T)

library(pander)
panderOptions('table.alignment.default',function(df) ifelse(sapply(df, is.numeric), 'right', 'left'))

library(dplyr)
library(ggplot2)
library(cowplot)

```

# Regression

## Introduction

The data we will use in the first part is a hypothetical study on child language acquisition (created by a [colleague of mine](http://coltekin.net/cagri/R/r-exercisesse8.html) and slightly adapted by me). We want to investigate the effects of amount of time spend in front of TV to two-year-old children's language development. The response variable in this data set, `cdi`, is a standard measure of children's language abilities based on parental reports. The predictor of interest is `tv.hours`, which is the weekly hours of TV time for each child.

You will have to fill in most commands yourself, but this should be feasible given the slides of the lecture, which can be viewed here: http://www.let.rug.nl/wieling/Statistics/Regression-large. While you can just enter the commands in RStudio, it is also possible to modify the source of this so-called R-markdown file directly in RStudio and press the "Knit HTML" button to generate an html file which contains both the commands you've used and their output. You can download the file to your current working directory in R by pasting the following command: `download.file('http://www.let.rug.nl/wieling/Statistics/Regression-large/lab/lab.Rmd', 'lab.Rmd')`. You can then open this file in RStudio. In this file, all R commands which are located within chunks (beginning and ending with three backticks) will be evaluated. Creating an R markdown file is very useful as your analysis becomes reproducible and easy to check for others. Note that chunks have options, with which you can customize the output. See for more information: http://www.rstudio.com/wp-content/uploads/2016/03/rmarkdown-cheatsheet-2.0.pdf

Note that the significance threshold $\alpha$ is set to 0.05 as usual.

## Importing the data
We will first download the data (which was saved as a dataframe in R: an `.rda` file), and load it into R.

```{r import}
download.file('http://www.let.rug.nl/wieling/Statistics/Regression-large/lab/tv.rda', 'tv.rda')
load('tv.rda') # an rda file can be loaded with the command load
```

## Structure of the data

Investigate the data with descriptive statistics and plots. Also calculate the correlations between the numerical variables.

### Summary of the data
```{r}
summary(tv)
```
From the summary we see it's a table of 80 observations with 40 female and 40 male. The cdi mean is 99.9 and the median 100 and the value go from 97 to 103. The mean of the hours spent in front of tv is 11.7 and the median is 11.0, being 2 hours the minmun value and 23 hours the maximum.

### Point plots for relation of the different columns and cdi score

Here we want to explore the possible effects of the collected data on the cdi score based on gender, by plotting in a points graph.
```{r}
plot1 <- tv%>%
  ggplot(aes(x = tv.hours, y = cdi, col = gender))+
  geom_point()+
  #geom_line(aes (y = mean(tv$cdi)), col = "grey") +
  labs(y = "cdi score", x = "hours in front of tv") +
  ggtitle('hours in tv')+
  theme(text = element_text(size=10))

plot2 <- tv%>%
  ggplot(aes(x = daycare.hours, y = cdi, col = gender))+
  geom_point()+
  #geom_line(aes (y = mean(tv$cdi)), col = "grey") +
  labs(y = "cdi score", x = "daycare hours") +
  ggtitle('daycare hours')+
  theme(text = element_text(size=10))

plot3 <- tv%>%
  ggplot(aes(x = book.reading, y = cdi, col = gender))+
  geom_point()+
  #geom_line(aes (y = mean(tv$cdi)), col = "grey") +
  labs(y = "cdi score", x = "book reading hours") +
  ggtitle('reading hours')+
  theme(text = element_text(size=10))

plot4 <- tv%>%
  ggplot(aes(x = mot.education, y = cdi, col = gender))+
  geom_point()+
  #geom_line(aes (y = mean(tv$cdi)), col = "grey") +
  labs(y = "cdi score", x = "mot.education") +
  ggtitle('mot.education')+
  theme(text = element_text(size=10))

plot_grid(plot1, plot2, plot3, plot4)
```
The first thing that seems clear by these plots is that male cdi scores are lower than female cdi scores. The increase of the hours in front of the tv seem to have a negative effect on the cdi score. The daycare hours don't show a clear effect on cdi score in this plot. Both the reading hours and the mot.education seem to have a positive effect on the cdi score.


### Differences between gender

Seeing the previous results, we want to confirm the existence of a difference between cdi scores based on gender. We calculate the mean and the median grouped by gender and boxplot the grouped data.
```{r}
# mean of cdi per gender

tv%>%
  group_by(gender)%>%
  summarise(mean= mean(cdi),
            median = median(cdi))

#quartiles

tapply(tv$cdi, tv$gender, quantile)

#boxplot

boxplot(tv$cdi ~ tv$gender, xlab="Gender", ylab="cdi")

```

We see that both the mean and the median are higher for female kids than for male kids. We can also see that quantile 1 of cdi score for females is the same value of quantile 3 for male kids.

### Correlation of the different values

```{r}

cor(tv[, c("cdi","tv.hours", "mot.education", "daycare.hours", "book.reading")])

```

With the new data we can confirm that the correllation between cdi and hours in tv is negative. We also see that mot.education and book.reading have a positive correlation with cdi and daycare.hours has a weak correlation with cdi. For the rest of the correlations we can higlight the positive relation between mot.education and book.reading (the meaning of the variable mot.education is not known).

## Regression modeling

In this section, we will fit the best model for the data, predicting the `cdi` language development score on the basis of various predictors.

### The best model without interactions

Fit the best model without taking into account interactions.
```{r lm}
m1 <- lm (cdi~tv.hours, tv)
m2 <- lm (cdi~tv.hours + gender, tv)
AIC(m1) - AIC(m2)
m3 <- lm (cdi~tv.hours + mot.education + gender, tv)
AIC(m2) - AIC(m3)
m4 <- lm (cdi~tv.hours + mot.education + book.reading + gender, tv)
AIC(m3) - AIC(m4)
m5 <- lm (cdi~tv.hours + mot.education + daycare.hours + book.reading + gender, tv)
AIC(m4) - AIC(m5)

summary(m4)
```
The AIC comparison is positive and higher than 2 for all the more complex models except from the last, in which daycare.hours is included. These results are in line with our expectations from the exploration part. Fort this reason, we exclude daycare.hours from our model and keep all the other variables.


### Visualize the results

```{r vis}
library(visreg)
visreg(m4)

```

As expected in the exploration of the data, tv.hours has a negative and significant effect in cdi score, while mot.education and book.reading have a positive effect in cdi score. Being a female has a positive effect in the cdi score.


### Assess the regression assumptions discussed in the lecture

#### Linearity
```{r linearity}
# linearity -> que la relació tingui una forma molt diferent a la lineal. ex: parabòlica
par(mfrow = c(2, 2))
library(car)
crPlot(m4, var = "tv.hours")
crPlot(m4, var = "mot.education")
crPlot(m4, var = "book.reading")
crPlot(m4, var = "gender")

```

The variables tv.hours and book.reading seem to have a linear effect on cdi score. Nontheless, the variable mot.education suggests a stagnation after 17. How to assess linearity in categorical variables??

#### Autocorrelation
```{r autocorrelation}
# autocorrelation in residuals
durbinWatsonTest(m4)

acf(resid(m4))
```

Autocorrelation in mainly important for time series data.

#### Multicollinearity
```{r multicollinearity}
# multicollinearity -> que unes variables independents expliquin les altres, llavors no són independents
vif(m4) #interpretar i arreglar si cal

summary(lm(tv.hours ~ mot.education, data = tv))$r.squared
summary(lm(tv.hours ~ book.reading, data = tv))$r.squared
summary(lm(book.reading ~ mot.education, data = tv))$r.squared
summary(lm(tv.hours ~ gender, data = tv))$r.squared
summary(lm(mot.education ~ gender, data = tv))$r.squared
summary(lm(book.reading ~ gender, data = tv))$r.squared

```
None of the variables used in the model explains more than 50% of another variable. For this reason, there is no collinearity.


#### Homoscedacity
```{r homoscedacity}

# homoscedasticity

plot(fitted(m4), resid(m4)) # interpretar i arreglar si cal
ncvTest(m4)
```
There seems to be no heteroscedacity looking at the graph. The p-value is too big to asses a significant heteroscedacity.

#### Distribution of residuals
```{r distribution}
# distribution of residuals

shapiro.test(resid(m4))$p.value
qqnorm(resid(m4))
qqline(resid(m4))

```
Residuals are normally distributed.

### Are interactions necessary?

Fit the best model while taking into account interactions. For simplicity and speed, we'll only investigate potential interactions with gender. Also visualize the results if an interaction is necessary.

```{r inter}

m6 <- lm (cdi~tv.hours + mot.education + book.reading + gender + tv.hours:gender, tv)
AIC(m4)-AIC(m6)
m7 <- lm (cdi~tv.hours + mot.education + book.reading + gender + mot.education:gender, tv)
AIC(m4)-AIC(m7)
m8 <- lm (cdi~tv.hours + mot.education + book.reading + gender + book.reading:gender, tv)
AIC(m4)-AIC(m8)

```
None of the AIC is lower than the one obtained by our best model. This indicates that interactions are not needed.

### Effect size

Obtain the effect size, both of the full model and the individual predictors

```{r eff}
library(lmSupport)
modelEffectSizes(m4)
summary(m4)
```
The model overall effect size of the model is of 0.636 (adjusted), which is a large effect. Individually, gender is the only predictor with a large effect (0.229).


### Model criticism

Apply model criticism by excluding observations with residuals outside 2.5 SD.

```{r mc}
tv2 <- tv[abs(scale(resid(m4))) < 2.5, ]
sum(abs(scale(resid(m4))) >= 2.5)

m42 <- lm(cdi~tv.hours + mot.education + book.reading + gender, tv2)
summary(m42)

```
Excluding observations with absolute residuals larger than 2.5 standard deviations above the mean deletes 2 observations from out data. The p-value doesn't change after excluding the outliers, and the patterns look the same.


### Bootstrap sampling

Conduct bootstrap sampling with 1000 repetitions.

```{r bs}
summary(m4.boot <- Boot(m42, R = 1000))


```
The results of the bootstrap sampling are very similar to the original correlation results of the model, which proofs that the model is robust.

# Logistic regression

## Introduction
In this assignment we will use logistic regression to assess the differences in pronunciation at the end of syllables in New York during the 1950s and 1960s depending on the socioeconomic class (phoneme [r] vs [&#601;] (schwa)).

We use the following data:


<table style="width:28%;"><colgroup><col width="13%"></col><col width="13%"></col></colgroup><thead><tr class="header"><th align="left">outcome</th><th align="left">status</th></tr></thead><tbody><tr class="odd"><td align="left">r</td><td align="left">upper</td></tr><tr class="even"><td align="left">r</td><td align="left">upper</td></tr><tr class="odd"><td align="left">r</td><td align="left">upper</td></tr><tr class="even"><td align="left">r</td><td align="left">upper</td></tr><tr class="odd"><td align="left">r</td><td align="left">upper</td></tr><tr class="even"><td align="left">r</td><td align="left">upper</td></tr><tr class="odd"><td align="left">r</td><td align="left">upper</td></tr><tr class="even"><td align="left">r</td><td align="left">upper</td></tr><tr class="odd"><td align="left">r</td><td align="left">upper</td></tr><tr class="even"><td align="left">r</td><td align="left">upper</td></tr><tr class="odd"><td align="left">r</td><td align="left">upper</td></tr><tr class="even"><td align="left">r</td><td align="left">upper</td></tr><tr class="odd"><td align="left">r</td><td align="left">upper</td></tr><tr class="even"><td align="left">r</td><td align="left">upper</td></tr><tr class="odd"><td align="left">r</td><td align="left">upper</td></tr><tr class="even"><td align="left">r</td><td align="left">upper</td></tr><tr class="odd"><td align="left">r</td><td align="left">upper</td></tr><tr class="even"><td align="left">r</td><td align="left">upper</td></tr><tr class="odd"><td align="left">r</td><td align="left">upper</td></tr><tr class="even"><td align="left">r</td><td align="left">upper</td></tr><tr class="odd"><td align="left">r</td><td align="left">upper</td></tr><tr class="even"><td align="left">r</td><td align="left">upper</td></tr><tr class="odd"><td align="left">r</td><td align="left">upper</td></tr><tr class="even"><td align="left">r</td><td align="left">upper</td></tr><tr class="odd"><td align="left">r</td><td align="left">upper</td></tr><tr class="even"><td align="left">r</td><td align="left">upper</td></tr><tr class="odd"><td align="left">r</td><td align="left">upper</td></tr><tr class="even"><td align="left">r</td><td align="left">upper</td></tr><tr class="odd"><td align="left">r</td><td align="left">upper</td></tr><tr class="even"><td align="left">r</td><td align="left">upper</td></tr><tr class="odd"><td align="left">schwa</td><td align="left">upper</td></tr><tr class="even"><td align="left">schwa</td><td align="left">upper</td></tr><tr class="odd"><td align="left">schwa</td><td align="left">upper</td></tr><tr class="even"><td align="left">schwa</td><td align="left">upper</td></tr><tr class="odd"><td align="left">schwa</td><td align="left">upper</td></tr><tr class="even"><td align="left">schwa</td><td align="left">upper</td></tr><tr class="odd"><td align="left">r</td><td align="left">middle</td></tr><tr class="even"><td align="left">r</td><td align="left">middle</td></tr><tr class="odd"><td align="left">r</td><td align="left">middle</td></tr><tr class="even"><td align="left">r</td><td align="left">middle</td></tr><tr class="odd"><td align="left">r</td><td align="left">middle</td></tr><tr class="even"><td align="left">r</td><td align="left">middle</td></tr><tr class="odd"><td align="left">r</td><td align="left">middle</td></tr><tr class="even"><td align="left">r</td><td align="left">middle</td></tr><tr class="odd"><td align="left">r</td><td align="left">middle</td></tr><tr class="even"><td align="left">r</td><td align="left">middle</td></tr><tr class="odd"><td align="left">r</td><td align="left">middle</td></tr><tr class="even"><td align="left">r</td><td align="left">middle</td></tr><tr class="odd"><td align="left">r</td><td align="left">middle</td></tr><tr class="even"><td align="left">r</td><td align="left">middle</td></tr><tr class="odd"><td align="left">r</td><td align="left">middle</td></tr><tr class="even"><td align="left">r</td><td align="left">middle</td></tr><tr class="odd"><td align="left">r</td><td align="left">middle</td></tr><tr class="even"><td align="left">r</td><td align="left">middle</td></tr><tr class="odd"><td align="left">r</td><td align="left">middle</td></tr><tr class="even"><td align="left">r</td><td align="left">middle</td></tr><tr class="odd"><td align="left">schwa</td><td align="left">middle</td></tr><tr class="even"><td align="left">schwa</td><td align="left">middle</td></tr><tr class="odd"><td align="left">schwa</td><td align="left">middle</td></tr><tr class="even"><td align="left">schwa</td><td align="left">middle</td></tr><tr class="odd"><td align="left">schwa</td><td align="left">middle</td></tr><tr class="even"><td align="left">schwa</td><td align="left">middle</td></tr><tr class="odd"><td align="left">schwa</td><td align="left">middle</td></tr><tr class="even"><td align="left">schwa</td><td align="left">middle</td></tr><tr class="odd"><td align="left">schwa</td><td align="left">middle</td></tr><tr class="even"><td align="left">schwa</td><td align="left">middle</td></tr><tr class="odd"><td align="left">schwa</td><td align="left">middle</td></tr><tr class="even"><td align="left">schwa</td><td align="left">middle</td></tr><tr class="odd"><td align="left">schwa</td><td align="left">middle</td></tr><tr class="even"><td align="left">schwa</td><td align="left">middle</td></tr><tr class="odd"><td align="left">schwa</td><td align="left">middle</td></tr><tr class="even"><td align="left">schwa</td><td align="left">middle</td></tr><tr class="odd"><td align="left">schwa</td><td align="left">middle</td></tr><tr class="even"><td align="left">schwa</td><td align="left">middle</td></tr><tr class="odd"><td align="left">schwa</td><td align="left">middle</td></tr><tr class="even"><td align="left">schwa</td><td align="left">middle</td></tr><tr class="odd"><td align="left">schwa</td><td align="left">middle</td></tr><tr class="even"><td align="left">schwa</td><td align="left">middle</td></tr><tr class="odd"><td align="left">schwa</td><td align="left">middle</td></tr><tr class="even"><td align="left">schwa</td><td align="left">middle</td></tr><tr class="odd"><td align="left">schwa</td><td align="left">middle</td></tr><tr class="even"><td align="left">schwa</td><td align="left">middle</td></tr><tr class="odd"><td align="left">schwa</td><td align="left">middle</td></tr><tr class="even"><td align="left">schwa</td><td align="left">middle</td></tr><tr class="odd"><td align="left">schwa</td><td align="left">middle</td></tr><tr class="even"><td align="left">schwa</td><td align="left">middle</td></tr><tr class="odd"><td align="left">schwa</td><td align="left">middle</td></tr><tr class="even"><td align="left">schwa</td><td align="left">middle</td></tr><tr class="odd"><td align="left">schwa</td><td align="left">middle</td></tr><tr class="even"><td align="left">schwa</td><td align="left">middle</td></tr><tr class="odd"><td align="left">schwa</td><td align="left">middle</td></tr><tr class="even"><td align="left">schwa</td><td align="left">middle</td></tr><tr class="odd"><td align="left">schwa</td><td align="left">middle</td></tr><tr class="even"><td align="left">schwa</td><td align="left">middle</td></tr><tr class="odd"><td align="left">schwa</td><td align="left">middle</td></tr><tr class="even"><td align="left">schwa</td><td align="left">middle</td></tr><tr class="odd"><td align="left">schwa</td><td align="left">middle</td></tr><tr class="even"><td align="left">schwa</td><td align="left">middle</td></tr><tr class="odd"><td align="left">schwa</td><td align="left">middle</td></tr><tr class="even"><td align="left">schwa</td><td align="left">middle</td></tr><tr class="odd"><td align="left">schwa</td><td align="left">middle</td></tr><tr class="even"><td align="left">schwa</td><td align="left">middle</td></tr><tr class="odd"><td align="left">schwa</td><td align="left">middle</td></tr><tr class="even"><td align="left">schwa</td><td align="left">middle</td></tr><tr class="odd"><td align="left">schwa</td><td align="left">middle</td></tr><tr class="even"><td align="left">schwa</td><td align="left">middle</td></tr><tr class="odd"><td align="left">schwa</td><td align="left">middle</td></tr><tr class="even"><td align="left">schwa</td><td align="left">middle</td></tr><tr class="odd"><td align="left">schwa</td><td align="left">middle</td></tr><tr class="even"><td align="left">schwa</td><td align="left">middle</td></tr><tr class="odd"><td align="left">schwa</td><td align="left">middle</td></tr><tr class="even"><td align="left">schwa</td><td align="left">middle</td></tr><tr class="odd"><td align="left">schwa</td><td align="left">middle</td></tr><tr class="even"><td align="left">schwa</td><td align="left">middle</td></tr><tr class="odd"><td align="left">schwa</td><td align="left">middle</td></tr><tr class="even"><td align="left">schwa</td><td align="left">middle</td></tr><tr class="odd"><td align="left">schwa</td><td align="left">middle</td></tr><tr class="even"><td align="left">schwa</td><td align="left">middle</td></tr><tr class="odd"><td align="left">schwa</td><td align="left">middle</td></tr><tr class="even"><td align="left">schwa</td><td align="left">middle</td></tr><tr class="odd"><td align="left">schwa</td><td align="left">middle</td></tr><tr class="even"><td align="left">schwa</td><td align="left">middle</td></tr><tr class="odd"><td align="left">schwa</td><td align="left">middle</td></tr><tr class="even"><td align="left">schwa</td><td align="left">middle</td></tr><tr class="odd"><td align="left">schwa</td><td align="left">middle</td></tr><tr class="even"><td align="left">schwa</td><td align="left">middle</td></tr><tr class="odd"><td align="left">schwa</td><td align="left">middle</td></tr><tr class="even"><td align="left">schwa</td><td align="left">middle</td></tr><tr class="odd"><td align="left">schwa</td><td align="left">middle</td></tr><tr class="even"><td align="left">schwa</td><td align="left">middle</td></tr><tr class="odd"><td align="left">r</td><td align="left">lower</td></tr><tr class="even"><td align="left">r</td><td align="left">lower</td></tr><tr class="odd"><td align="left">r</td><td align="left">lower</td></tr><tr class="even"><td align="left">r</td><td align="left">lower</td></tr><tr class="odd"><td align="left">schwa</td><td align="left">lower</td></tr><tr class="even"><td align="left">schwa</td><td align="left">lower</td></tr><tr class="odd"><td align="left">schwa</td><td align="left">lower</td></tr><tr class="even"><td align="left">schwa</td><td align="left">lower</td></tr><tr class="odd"><td align="left">schwa</td><td align="left">lower</td></tr><tr class="even"><td align="left">schwa</td><td align="left">lower</td></tr><tr class="odd"><td align="left">schwa</td><td align="left">lower</td></tr><tr class="even"><td align="left">schwa</td><td align="left">lower</td></tr><tr class="odd"><td align="left">schwa</td><td align="left">lower</td></tr><tr class="even"><td align="left">schwa</td><td align="left">lower</td></tr><tr class="odd"><td align="left">schwa</td><td align="left">lower</td></tr><tr class="even"><td align="left">schwa</td><td align="left">lower</td></tr><tr class="odd"><td align="left">schwa</td><td align="left">lower</td></tr><tr class="even"><td align="left">schwa</td><td align="left">lower</td></tr><tr class="odd"><td align="left">schwa</td><td align="left">lower</td></tr><tr class="even"><td align="left">schwa</td><td align="left">lower</td></tr><tr class="odd"><td align="left">schwa</td><td align="left">lower</td></tr><tr class="even"><td align="left">schwa</td><td align="left">lower</td></tr><tr class="odd"><td align="left">schwa</td><td align="left">lower</td></tr><tr class="even"><td align="left">schwa</td><td align="left">lower</td></tr><tr class="odd"><td align="left">schwa</td><td align="left">lower</td></tr><tr class="even"><td align="left">schwa</td><td align="left">lower</td></tr><tr class="odd"><td align="left">schwa</td><td align="left">lower</td></tr><tr class="even"><td align="left">schwa</td><td align="left">lower</td></tr><tr class="odd"><td align="left">schwa</td><td align="left">lower</td></tr><tr class="even"><td align="left">schwa</td><td align="left">lower</td></tr><tr class="odd"><td align="left">schwa</td><td align="left">lower</td></tr><tr class="even"><td align="left">schwa</td><td align="left">lower</td></tr><tr class="odd"><td align="left">schwa</td><td align="left">lower</td></tr><tr class="even"><td align="left">schwa</td><td align="left">lower</td></tr><tr class="odd"><td align="left">schwa</td><td align="left">lower</td></tr><tr class="even"><td align="left">schwa</td><td align="left">lower</td></tr><tr class="odd"><td align="left">schwa</td><td align="left">lower</td></tr><tr class="even"><td align="left">schwa</td><td align="left">lower</td></tr><tr class="odd"><td align="left">schwa</td><td align="left">lower</td></tr><tr class="even"><td align="left">schwa</td><td align="left">lower</td></tr><tr class="odd"><td align="left">schwa</td><td align="left">lower</td></tr><tr class="even"><td align="left">schwa</td><td align="left">lower</td></tr><tr class="odd"><td align="left">schwa</td><td align="left">lower</td></tr><tr class="even"><td align="left">schwa</td><td align="left">lower</td></tr><tr class="odd"><td align="left">schwa</td><td align="left">lower</td></tr><tr class="even"><td align="left">schwa</td><td align="left">lower</td></tr><tr class="odd"><td align="left">schwa</td><td align="left">lower</td></tr><tr class="even"><td align="left">schwa</td><td align="left">lower</td></tr><tr class="odd"><td align="left">schwa</td><td align="left">lower</td></tr><tr class="even"><td align="left">schwa</td><td align="left">lower</td></tr><tr class="odd"><td align="left">schwa</td><td align="left">lower</td></tr><tr class="even"><td align="left">schwa</td><td align="left">lower</td></tr><tr class="odd"><td align="left">schwa</td><td align="left">lower</td></tr><tr class="even"><td align="left">schwa</td><td align="left">lower</td></tr></tbody></table>

```{r, results='asis', echo=F, eval=F}
status <- c(rep("upper", 36), rep("middle",94), rep("lower",54))
pronunciation <- c(rep("r", 30), rep("schwa", 6), rep("r", 20), rep("schwa", 74), rep("r", 4), rep("schwa", 50))

# create the data frame
dat <- data.frame(pronunciation, status)
pandoc.table(dat)
```

## Create data frame

Use R code to generate this data without having to type all the observations (there are 184 observations...). Hint: you can use the function `c()` to concatenate information and function `rep()` to repeat one element multiple times, e.g. `rep("upper",36)` will create a vector in which `"upper"` is repeated 36 times. In the end you should have a data frame with 1 dependent variable (`pronunciation`, with 2 levels: `schwa` and `r`) and 1 predictor (`status`, with 3 levels: `lower`, `middle` and `upper`).

```{r createtable}
status <- c(rep("upper", 36), rep("middle",94), rep("lower",54))
pronunciation <- c(rep("r", 30), rep("schwa", 6), rep("r", 20), rep("schwa", 74), rep("r", 4), rep("schwa", 50))

# create the data frame
dat <- data.frame(pronunciation, status)
pandoc.table(dat)

```

Furthermore, make sure that both variables are coded as factors and show that this indeed is the case. Set `lower` as your reference level for the predictor `status`. Now add an additional predictor `IsR` which is equal to 1 if `pronunciation` equal [r] and 0 otherwise. This will be our dependent variable.

```{r factors}
levels(dat$status)
levels(dat$pronunciation)

library(stringr)

for (i in dat){
  if (str_detect(pronunciation[i], "r")){
    dat$Isr <- 1
  } else{
    dat$Isr <- 0
  }
}

summary(dat)

```

## Logistic regression model

Now fit a binomial logistic regression model to the data.

```{r logmod}
#m <- glm(SpeciesVirginica ~ Petal.Area, data = iris, family = "binomial")

```

## Visualize the results

Now visualize the results.

```{r vislog}
# your code goes here
```


## Assumptions

Check the logistic regression assumptions.

```{r checkassumptions}
# your code goes here
```

## Goodness of fit

Evaluate how well the model performs, by calculating the index of concordance.

```{r goodnessoffit}
# your code goes here
```

## Interpretation

Given your model, what are the odds of `r` to `schwa` when the predictor `status` is at its reference level (`lower`)?
Hint: the estimates in a logistic regression model are in logits, which is the natural logarithm of the odds. To back transform to odds, you'll need to invert the logartithm (how?). Add the R code which calculates these odds. Also calculate the *probability* of the outcome being `r` when status is at its reference level (`lower`).

```{r interpretation1}
# your code goes here
```

Which `status` is associated with a greater probability of pronouncing the [r], instead of [&#601;]? What are the odds of hearing [r] with this status as opposed to the `lower` status?

```{r interpretation2}
# your code goes here
```

# Answers
The answers to the questions in this file can be viewed here: http://www.let.rug.nl/wieling/Statistics/Regression-large/lab/answers. The associated R markdown file can be downloaded here: http://www.let.rug.nl/wieling/Statistics/Regression-large/lab/answers/answers.Rmd

# Replication

From within RStudio, you can simply download this file using the commands:
```{r knitrs, eval=F}
# download original file if not already exists (to prevent overwriting)
#if (!file.exists('lab.Rmd')) {
#  download.file('http://www.let.rug.nl/wieling/Statistics/Regression-large/lab/lab.Rmd', 'lab.Rmd')
#}
```
Subsequently, open it in the editor and use the Knit HMTL button to generate the html file.

If you use plain R, you first have to install [Pandoc](http://johnmacfarlane.net/pandoc/). Then copy the following lines to the most recent version of R.

```{r knitr, eval=F}
# install rmarkdown package if not installed
#if(!"rmarkdown" %in% rownames(installed.packages())) {
#   install.packages("rmarkdown")
#}
#library(rmarkdown) # load rmarkdown package

# download original file if not already exists (to prevent overwriting)
#if (!file.exists('lab.Rmd')) {
#  download.file('http://www.let.rug.nl/wieling/Statistics/Regression-large/lab/lab.Rmd', 'lab.Rmd')
#}

# generate output
#render('lab.Rmd') # generates html file with results

# view output in browser
#browseURL(paste('file://', file.path(getwd(),'lab.html'), sep='')) # shows result
```