
** [example_school.Rmd](example_school.Rmd) **

# Example school

Real data from

[source](https://www.youtube.com/watch?v=LzAwEKrn2Mc)

[data source](https://cran.r-project.org/web/packages/merTools/merTools.pdf)

```
hsb A subset of data from the 1982 High School and Beyond survey used
as examples for HLM software

Format
A data frame with 7,185 observations on the following 8 variables.
schid a numeric vector, 160 unique values
mathach a numeric vector for the performance on a standardized math assessment
female a numeric vector coded 0 for male and 1 for female
ses a numeric measure of student socio-economic status
minority a numeric vector coded 0 for white and 1 for non-white students
schtype a numeric vector coded 0 for public and 1 for private schools
meanses a numeric, the average SES for each school in the data set
size a numeric for the number of students in the school
```

schid: school id
minority:
female: level 1 predictor
ses: level 1 predictor, ?socioeconomic status 
mathach: mathematics achievement, level 1 outcome variable
schtype: level 2 predictor
meanses: level 2 predictor

mathach is our reaction variable

```{r}
# we load library merTools that includes dataset hsb
library(merTools)
# we define dd, a data object containing the data of hsb, just for convenience. 
# We could have used hsb also, a data object that exists in the merTools namespace
dd <- hsb
# we take a look at our reaction variable using describe() of library(psych)
psych::describe(dd$mathach - min(dd$mathach))

# we inspect the nominal scaled variables and their cominations
unique(dd[,c("minority", "female", "schtype")])
# hm, this score seems to qualify the better the achievement the higher the score
# ... but what means a negative score?
# anyway
# to calculate a variation index to be explained, we use the distance to 0
dd$diff.a <- dd$mathach - min(dd$mathach)
var.a <- sum(dd$diff.a ^ 2) / length(dd$diff.a)
sd.a <- sqrt(var.a)
var.mean <- var(dd$mathach)
sd.mean <- sqrt(var.mean)

# we see a reduction
cat(sd.a, " > ", sd.mean)
# we could conclude, that arithmetic mean is a better estimator than 0 to predict our reaction variable in new observations
# ... but this comparison is sort of arbitrarily, as there is no accepted way of a value to calculate the individual diff 

# in statisctics we therefore use the usual variance/sd and the arithmetic mean to compare the introduction of more parameters in the model
# we call this often our NULL-model
# in model formula this is defined by a explanatory variable that is constant, usually 1
m.0 <- lm(dd$mathach ~ 1)
sd(m.1$residuals)
# ... standard deviation of the residuals (unexplained variation of reaction variable) is the same as above
cat(sd.a, " > ", sd(m.0$residuals))

# we plot that pseudo reduction
library(tidyverse)
psych::describe(m.1$residuals)
#tibble(sd.a = sqrt(dd$diff.a ^ 2), sd.0 = sqrt(m.0$residuals ^ 2)) %>% 
se <- function(x) sqrt(var(x)/length(x))
tibble(d.sd.a = dd$diff.a, d.sd.0 = m.0$residuals) -> dd.g0
#tibble(sd.a = sqrt(dd$diff.a ^ 2), sd.0 = m.0$residuals) %>% 

dd.g <- tibble(
  group = as.factor(c('a', 'm.0')),
	mm = c(0, 0),
  se = c(sum(sqrt(dd$diff.a ^ 2)) / length(dd$diff.a), sd(m.0$residuals))
)
  
ggplot(dd.g, aes(x=group)) +
  geom_point(data=dd.g, stat='identity', aes(group, y=mm, size=10), show.legend=F) +
	geom_errorbar(data=dd.g, stat='identity', aes(x=group, ymax = mm + se, ymin = mm - se), width = 0.2)


```



Version: `r format(Sys.time(), '%d %B, %Y %H:%M')`
