---
title: "Homework 10 Solutions"
date: "11/17/2021"
output: 
  pdf_document:
    keep_tex: true
fig_caption: true
header-includes:
  - \usepackage{amsmath}
  - \usepackage{amsthm}
  - \DeclareMathOperator{\Cov}{Cov}
  - \DeclareMathOperator{\Corr}{Corr}
  - \DeclareMathOperator{\Var}{Var}
  - \DeclareMathOperator{\upmodels}{\perp \!\!\! \perp}
---


# Problem 1

```{r message = F, warning = F} 
library(tidyverse)
library(arm)
library(lme4)
library(latex2exp) # for nice plot labeling (using latex notation)

cdi <- read.table(file = "../cdi.dat")

cdi$pci <- rescale(2*cdi$per.cap.income) 
cdi$phg <- rescale(2*cdi$pct.hs.grad)
```

## (a)

See Figure \ref{f:fig1a}.

```{r message = F, warning = F}
# pooled regression
lm_all <- lm(pci ~ phg, data = cdi)

lm_all_rep <- data.frame(intercept = lm_all[["coefficients"]][1],
                         slope = lm_all[["coefficients"]][2],
                         state = unique(cdi$state))

# unpooled regression
state_lm_func <- function(df){
    mod <- lm(pci ~ phg, data = df)
    return(data.frame(intercept = mod[["coefficients"]][1], 
                      slope = mod[["coefficients"]][2]))
}

state_lm_df <- cdi %>%
  group_by(state) %>% do(., state_lm_func(.))

## unpooled regression version II
lm_unpooled <- lm(pci ~ state + phg:state - 1, data = cdi) # single model

coef_unpooled <- coef(lm_unpooled)
state_loc <- substring(names(coef_unpooled), 6,7)
slope_loc <- substring(names(coef_unpooled), 9,11)

intercepts <- coef_unpooled[slope_loc == ""] 
intercepts_df <- data.frame(state = substring(names(intercepts), 6,7),
                            intercept = intercepts)
slopes <- coef_unpooled[slope_loc == "phg"]
slopes_df <- data.frame(state = substring(names(slopes), 6,7),
                        slope = slopes)

state_lm_df2 <- intercepts_df %>% left_join(slopes_df, by = 'state')

testthat::expect_equivalent(state_lm_df, state_lm_df2) # no error - so true

## unpooled regression version III
lm_unpooled2 <- lm(pci ~ state*phg - 1 - phg, data = cdi) # single model

coef_unpooled2 <- coef(lm_unpooled2)
state_loc2 <- substring(names(coef_unpooled2), 6,7)
slope_loc2 <- substring(names(coef_unpooled2), 9,11)

intercepts2 <- coef_unpooled2[slope_loc2 == ""] 
intercepts_df2 <- data.frame(state = substring(names(intercepts2), 6,7),
                             intercept = intercepts2)
slopes2 <- coef_unpooled2[slope_loc2 == "phg"]
slopes_df2 <- data.frame(state = substring(names(slopes2), 6,7),
                         slope = slopes2)

state_lm_df3 <- intercepts_df2 %>% left_join(slopes_df2, by = 'state')

testthat::expect_equivalent(state_lm_df, state_lm_df3) # no error - so true


# multilevel model
multilevel_lmer <- lmer(formula = pci ~ 1 + phg + (1 | state) + (0 + phg | state), 
                        data = cdi)

multilevel_lmer_df <- coef(multilevel_lmer)$state %>%
  tibble::rownames_to_column() %>%
  rename(state = "rowname",
         intercept = "(Intercept)",
         slope = "phg")

# multilevel model version II
betas <- fixef(multilevel_lmer)
etas <- ranef(multilevel_lmer)$state

combo_multilevel_lmer_df <- etas + 
  matrix(rep(betas, each = nrow(etas)), nrow = nrow(etas))

testthat::expect_equivalent(combo_multilevel_lmer_df,
                            coef(multilevel_lmer)$state) # doesn't error if true

# combined
combined <- rbind(mutate(lm_all_rep[, c("state", "intercept", "slope")], 
                         id = "pooled regression"),
                  mutate(state_lm_df[, c("state", "intercept", "slope")], 
                         id = "unpooled regression"),
                  mutate(multilevel_lmer_df[, c("state", "intercept", "slope")],
                         id = "multilevel model")) %>%
  mutate(id = factor(id, levels = c("pooled regression",
                                    "unpooled regression",
                                    "multilevel model")))
        
```


```{r fig.height=7.5, message= F, warning = F, fig.cap = "Problem 1a: Different linear models to compare per capital income to percent high school graduates \\label{f:fig1a}."}
ggplot(cdi) +
  geom_point(aes(x = phg, y= pci)) +
  geom_abline(data = combined,
              aes(slope = slope, intercept= intercept, group = state,
                  color = id))+
  facet_wrap(~ state) +
  theme_bw() +
  labs(color = "model type")
```

## (b)

See figure \ref{f:fig1b}.

```{r}
# multilevel model
multilevel_res_fit <- data.frame(fitted = fitted(multilevel_lmer), 
                                 residual = resid(multilevel_lmer), 
                                 state = cdi$state)

# multilevel model, part II
source("../residual-functions.r")
multilevel_res_fit2 <- data.frame(fitted = yhat.cond(multilevel_lmer),
                                  residual = r.cond(multilevel_lmer),
                                  state = cdi$state)

testthat::expect_equivalent(multilevel_res_fit, multilevel_res_fit2) # no error if true
```

```{r, fig.height=7.5, message = F, warning = F, fig.cap = "Problem 1b: Examinig the conditional residual vs conditional $\\hat{y}$ for multilevel model \\label{f:fig1b}."}
multilevel_res_fit %>%
ggplot() +
  geom_point(aes(x = fitted, y = residual)) +
  facet_wrap(~ state) +
  geom_hline(yintercept = 0) +
  theme_bw() +
  labs(y = "conditional residual",
       x = TeX("conditional $\\hat{y}$"))
```

\clearpage

# Problem 2

NOTES: 

(1) Instead of giving "full" solutions for #2, I'm just going to provide some suggestive code...  Mostly I am not providing interpretations-- you will need to do that in your IDMRAD paper for the college, after you are satisfied with your technical appendix!

(2) There are many pages here because I couldn't suppress some of the output from fitLMER.fnc() from library(LMERconveniencefunctions).

```{r}
library(lme4)
library(arm)
library(ggplot2)

##
## Read the data in wide and tall formats...
ratings <- read.csv("ratings.csv",header=T)
tall <-  read.csv("tall.csv",header=T)

##
## Some preliminary exploration of the data with View(), summary(), table(), etc.,
## omitted here -- BUT it is important that you do it, so you can catch things 
## like missing data (however it is coded), not all raters using all categories, etc.

##
## Take care that all ratings run from 1 to 4, 
## whether or not rater used all categories...
tall$Rating <- factor(tall$Rating,levels=1:4)
for (i in unique(tall$Rubric)) {
  ratings[,i] <- factor(ratings[,i],levels=1:4)
}

##
## Note that in the "ratings" data frame, the missing "Sex"
## value is "--" while in the "tall" data frame it is ""
## (a string of length 0).
##
## Make the "tall" be consistent with the "ratings" coding.
tall$Sex[nchar(tall$Sex)==0] <- "--"

##
## Extract the reduced data set with the 13 artifacts that all 3 raters saw...
ratings.13 <- ratings[grep("O",ratings$Artifact),]
tall.13 <-  tall[grep("O",tall$Artifact),]
```


## 2(a) Is the distribution of ratings for each rubrics pretty much indistinguishable from the other rubrics, or are there rubrics that tend to get especially high or low ratings? Is the distribution of ratings given by each rater pretty much indistinguishable from the other raters, or are there raters that tend to give especially high or low ratings?


Here are some ideas to compare distributions across Rubrics.  I leave other possibilities, pretty formatting (wtih kable, etc.), and smart interpretations up to you...
```{r}
##
## Bar plots for the reduced data set
g <- ggplot(tall.13,aes(x = Rating)) +
  facet_wrap( ~ Rubric) +
  geom_bar()

g

##
## Table of counts might make a nice supplement
tmp <- data.frame(lapply(split(tall.13$Rating,tall.13$Rubric),summary))
row.names(tmp) <- paste("Rating",1:4)

tmp

##
## Barplots for full data set
g <- ggplot(tall,aes(x = Rating)) +
  facet_wrap( ~ Rubric) +
  geom_bar()

g

##
## Table of counts again.  A bit pesky since there are NA's...
tmp0 <- lapply(split(tall$Rating,tall$Rubric),summary)
tmp <- data.frame(matrix(0,nrow=5,ncol=7))  ## seven rubrics...
names(tmp) <- names(tmp0)
row.names(tmp) <- c(paste("Rating",1:4),"<NA>")
for (i in names(tmp0)) {
  tmp[,i] <- tmp[,i] + c(tmp0[[i]],0)[1:5]
}

tmp

```
And here are some idea to compare distributions across Raters.  Again, I leave other possibilities, attractive formatting (with kable, changing height of graphs, etc.), and smart intepretations up to you...

```{r}

##
## Needed to make the title of each facet more human-readable...
rater.name <- function(x) { paste("Rater",x) }

## 
## Barplots for reduced data...
g <- ggplot(tall.13,aes(x = Rating)) +
  facet_wrap( ~ Rater, labeller=labeller(Rater=rater.name)) +
  geom_bar()

g

##
## Corresponding table of counts...
tmp <- data.frame(lapply(split(tall.13$Rating,tall.13$Rater),summary))
row.names(tmp) <- paste("Rating",1:4)
names(tmp) <- paste("Rater",1:3)

tmp

##
## Barplots for full data...
g <- ggplot(tall,aes(x = Rating)) +
  facet_wrap( ~ Rater, labeller=labeller(Rater=rater.name)) +
  geom_bar()

g

##
## Corresponding table of counts...
tmp0 <- lapply(split(tall$Rating,tall$Rater),summary)
tmp <- data.frame(matrix(0,nrow=5,ncol=3))  ## three raters...
names(tmp) <- names(tmp0)
row.names(tmp) <- c(paste("Rating",1:4),"<NA>")
for (i in names(tmp0)) {
  tmp[,i] <- tmp[,i] + c(tmp0[[i]],0)[1:5]
}
names(tmp) <- paste("Rater",1:3)
tmp

```

Now, those NA's have me curious...

```{r}
tall[apply(tall,1,function(x){any(is.na(x))}),]

ratings[ratings$Sex=="--",]
```
First, note that none of the missing values occur in the smaller 13-rubric data set (how can we tell this from the output above?).  So we don't have to worry about missing data at all in analyses that just involve this smaller data set.

Second, in any modeling that we do, the "Rating" is the outcome variable, so R will just drop the two observations with missing Rating values.  This will mean that the "full" data sets may be different for models that involve different rubrics: For models involving five of the rubrics we will get all the data from all the raters, but for models involving CritDes we would be missing a rating from Rater 2, and for models involving VisOrg we would be missing a rating from Rater 1.  We need to be vigilant about when these differences actually occur, since they could undermine some model comparisons (different data sets).

Third, we will also have to be careful of the missing "Sex" value (currently coded as "--".  If we coded it as NA, then R would drop it from models that have Sex as a predictor, which would make comparing models with and without Sex as a predictor more difficult (different data sets!). We could just drop this student from all analyses, but it seems like a waste to lose that data.  We could code it as "F" or "M" if we had a convincing justification for doing so, but since I don't have convincing justification (do you??), I'm just going to leave it as a third "Sex" category for now...    

## 2(b) For each rubric, do the raters generally agree on their scores?  If not, is there one rater who disagrees with the others?  Or do they all disagree?

Here are some coding ideas.  NOTE: Once again, I leave other ideas, attractive formatting, and interpretations (actually answering the question!) to you...

```{r}

##
## useful preliminaries
Rubric.names <- sort(unique(tall$Rubric))

##
## First we examine the 13 "common" artifacts that all 3 raters saw...

## 
## Note: extracting sig^2 and tau^2 from the fitted lmer() object took a little
## spelunking...  First I fitted a "test model"
##
## tmp <- lmer(as.numeric(Rating) ~ 1 + (1|Artifact),data=tall.13[tall.13$Rubric=="RsrchQ",])
##
## then I looked at names(summary(tmp)) and tried to see what was under each name
## by looking at summary(tmp)$methTitle, summary(tmp)$objClass, etc. for all the 
## names I found.  I quickly found summary(tmp)$sigma, which can be squared to get 
## sig^2.  It took more exploring with summary(tmp)$varcor to get tau^2...
##

ICC.vec <- NULL
for (i in Rubric.names) {
  
  tmp <- lmer(as.numeric(Rating) ~ 1 + (1|Artifact), data=tall.13[tall.13$Rubric==i,])
  sig2 <- summary(tmp)$sigma^2
  tau2 <- attr(summary(tmp)$varcor[[1]],"stddev")^2 
  ICC <- tau2 / (tau2 + sig2)
  ICC.vec <- c(ICC.vec,ICC) 
}
names(ICC.vec) <- Rubric.names

agreement.results <- cbind(ICC.common=ICC.vec,"       a12"=0,a23=0,a13=0)

agreement.tables <- as.list(rep(NA,7))
names(agreement.tables) <- Rubric.names

for (i in Rubric.names) {
  r12 <- data.frame(r1=factor(ratings.13[ratings.13$Rater==1,i],levels=1:4),
                    r2=factor(ratings.13[ratings.13$Rater==2,i],levels=1:4),
                    a1=ratings.13[ratings.13$Rater==1,"Artifact"],
                    a2=ratings.13[ratings.13$Rater==2,"Artifact"])
  if(any(r12[,3]!=r12[,4])) { stop(paste("Rater 1-2 Artifact mismatch on rubric",i)) }
  a12 <- mean(r12[,1]==r12[,2])
  r12 <- table(r12[,1:2])  ## print this to see how much agreement there is among raters 1-2

  r23 <- data.frame(r2=factor(ratings.13[ratings.13$Rater==2,i],levels=1:4),
                    r3=factor(ratings.13[ratings.13$Rater==3,i],levels=1:4),
                    a2=ratings.13[ratings.13$Rater==2,"Artifact"],
                    a3=ratings.13[ratings.13$Rater==3,"Artifact"])
  if(any(r23[,3]!=r23[,4])) { stop(paste("Rater 2-3 Artifact mismatch on rubric",i)) }
  a23 <- mean(r23[,1]==r23[,2])
  r23 <- table(r23[,1:2])  ## print this to see how much agreement there is among raters 2-3

  r13 <- data.frame(r1=factor(ratings.13[ratings.13$Rater==1,i],levels=1:4),
                    r3=factor(ratings.13[ratings.13$Rater==3,i],levels=1:4),
                    a1=ratings.13[ratings.13$Rater==1,"Artifact"],
                    a3=ratings.13[ratings.13$Rater==3,"Artifact"])
  if(any(r13[,3]!=r13[,4])) { stop(paste("Rater 1-3 Artifact mismatch on rubric",i)) }
  a13 <- mean(r13[,1]==r13[,2])
  r13 <- table(r13[,1:2])  ## print this to see how much agreement there is among raters 1-3
  
  agreement.results[i,2:4] <- c(a12,a23,a13)
  
  agreement.tables[[i]] <- list(r12,r23,r13) 
  
}
round(agreement.results,2)

##

if (F) { print(agreement.tables) } 

## change to "if (T)" to get a crude printout of all tables...

## 
## Now add in ICC's calculated from all the data...

ICC.vec <- NULL
for (i in Rubric.names) {
  
  tmp <- lmer(as.numeric(Rating) ~ 1 + (1|Artifact), data=tall[tall$Rubric==i,])
  sig2 <- summary(tmp)$sigma^2
  tau2 <- attr(summary(tmp)$varcor[[1]],"stddev")^2 
  ICC <- tau2 / (tau2 + sig2)
  ICC.vec <- c(ICC.vec,ICC) 
}
names(ICC.vec) <- Rubric.names

agreement.results <- cbind(ICC.alldata=ICC.vec,agreement.results)

round(agreement.results,2)

```



## 2(c) More generally, how are the various factors in this experiement (Rater, Semester, Sex, Repeated, Rubric) related to the ratings?  Do the factors interact in any interesting ways?

There is a lot here, so I will break it up into "subsections":

* 2(c)(i): Adding fixed effects to the seven rubric-specific models using just the data from the 13 common artifacts that al three raters saw

* 2(c)(ii): Adding fixed effects to the seven rubric-specific models using all the data

* 2(c)(iii): Trying interactions and new random effects for the seven rubric specific models using all the data

* 2(c)(iv): Trying to add fixed effects, interactions, and new random effects to the "combined" model Rating ~ 1 + (0 + Rubric|Artifact), using all the data.

I want to note that some of the output below is pretty verbose (that's why there are so many pages!), and some of it runs off the edge of the page.  I dislike both of these (verbose output, and output you can't read because it runs off the page) but I haven't been able to fix it.  Some R functions just do not respect page layout attributes, and some functions don't have a 'turn off verbose output' option.

Also, some of the models take lmer considerable time to fit.  You may think your computer has hung.  It probably has not hung; some models take a minute or two to fit.  If anything is taking more than 5 minutes, that's the time to consider whether you really need that model. 
 
### 2(c)(i): Adding fixed effects to the seven rubric-specific models using just the data from the 13 common artifacts that al three raters saw

First, we try to add fixed effects to our seven rubric-specific models...  In principle it will matter whether we use only the data reduced to the 13 common artifacts, or the full data set. 

I will start with the reduced data (so of course I can't check "repeated" on this reduced data---why not?).

```{r message=FALSE, warning=FALSE}
library(LMERConvenienceFunctions)
library(RLRsim)

## Experiments before trying "production" code:
##
## I started by fitting a single model and trying fitLMER.fnc() on it. 
## 
## As stated on Piazza, fitLMER.fnc() doesn't seem to like intecept-only
## as the final model, and so for models including rater, I removed 
## the intercept, so that (effectively) rater would always be in the 
## model, and fitLMER.fnc() wouldn't complain.
##
## So my starting model for experimenting was

tmp <- lmer(as.numeric(Rating) ~ -1 + as.factor(Rater) + 
              Semester + Sex + (1|Artifact),
            data=tall.13[tall.13$Rubric=="RsrchQ",],REML=FALSE)

##
## Since backwards-elimination always involves nested models, I could
## use t-tests, F-tests or likelihood ratio tests to eliminate fixed
## effects.  The default for fitLMER.fnc() is to use t-tests with a 
## threshold of 2 (cutoff for the t-statistic, rather than a cutoff
## like 0.05 for the p-value). This is good enough for me.  I will also 
## force fitLMER.fnc() to fit using ML rather than REML, so the
## t-tests are as close to correct as I can get.
##
## So a typical function call would be 

tmp.back_elim <- fitLMER.fnc(tmp,set.REML.FALSE = TRUE,log.file.name = FALSE)

## (It would be nice to eliminate the verbose narrative here, but I don't see 
## a way to do that.  Just stuck with it I guess...)

## Anyway, backwards elimination with fitLMER.fnc() yields a model 
## with raters only:

formula(tmp.back_elim)

## (Note that raters will always remain in the model, since raters have to rate
## in a category >= 1, and the t-tests compares each rater's average rating to 0. 
## Unless the sample size is really small, this should always yield a significant 
## coefficient for each rater dummy variable in the model.)

## The estimates for raters don't look that different from each other, 
## so we can test to see if they are different by comparing with the 
## intercept-only model

tmp.int_only <- update(tmp.back_elim, . ~ . + 1 - as.factor(Rater))

anova(tmp.int_only,tmp.back_elim)

## Again the models are nested so I really only need to look at the p-value
## from the likelihood ratio chi-squared test.  A little fooling around with
## names(anova(tmp.int_only,tmp.back_elim)), etc., shows that I can get this 
## as 

anova(tmp.int_only,tmp.back_elim)$"Pr(>Chisq)"[2]

## it looks like the intercept-only model is adequate here (the p-value 
## is much greater than 0.05 or any other common significance level).

## Note: since no main effects were retained, there's really no reason to 
## check for interactions.

## Now, I need to code this into a loop so I don't have to do everything by hand...


Rubric.names <- sort(unique(tall$Rubric))

model.formula.13 <- as.list(rep(NA,7))
names(model.formula.13) <- Rubric.names

## There will be a lot of output from fitLMER.fnc() here... Sorry!

for (i in Rubric.names) {
  
  ## fit each base model
  rubric.data <- tall.13[tall.13$Rubric==i,]
  tmp <- lmer(as.numeric(Rating) ~ -1 + as.factor(Rater) + 
              Semester + Sex + (1|Artifact),
            data=rubric.data,REML=FALSE)
  
  ## do backwards elimination
  tmp.back_elim <- fitLMER.fnc(tmp,set.REML.FALSE = TRUE,log.file.name = FALSE)
  
  ## check to see if the raters are significantly different from one another
  tmp.single_intercept <- update(tmp.back_elim, . ~ . + 1 - as.factor(Rater))
  pval <- anova(tmp.single_intercept,tmp.back_elim)$"Pr(>Chisq)"[2]
  
  ## choose the best model
  if (pval<=0.05) {
    tmp_final <- tmp.back_elim
  } else {
    tmp_final <- tmp.single_intercept
  }

  ## and add to list...
  model.formula.13[[i]] <- formula(tmp_final)
  
}

## see what "final models" we got...
model.formula.13
```
So, it looks like we don't need to add any fixed effects or interactions to the models for each rubric, using only the data reduced to the 13 common rubrics.

### 2(c)(ii): Adding fixed effects to the seven rubric-specific models using all the data

Now let's try with the full data...

```{r}
Rubric.names <- sort(unique(tall$Rubric))

## Note: Now the missing ratings become important.  We want to use the same data
## set for every model fit and model comparison.  I am going to eliminate by 
## hand the two observations with missing data, and only do fitting and comparison
## on this "slightly" reduced data set.  

tall[c(161,684),] ## just to check that these are the rows with missing ratings...
tall.nonmissing <- tall[-c(161,684),]  ## now delete them...

## I can't think of a good justification for imputing the "Sex" of the student who
## didn't report this to either M or F, and leaving it as "--" makes the models
## harder to interpret.  So I will eliminate that person from the data set also...

tall.nonmissing[tall.nonmissing$Sex=="--",]  ## check which rows will be eliminated
tall.nonmissing <- tall.nonmissing[tall.nonmissing$Sex!="--",] ## eliminate them

model.formula.alldata <- as.list(rep(NA,7))
names(model.formula.alldata) <- Rubric.names

## There will be a lot of output from fitLMER.fnc() here... Sorry!

for (i in Rubric.names) {
  
  ## fit each base model
  rubric.data <- tall.nonmissing[tall.nonmissing$Rubric==i,]
  tmp <- lmer(as.numeric(Rating) ~ -1 + as.factor(Rater) + 
              Semester + Sex + (1|Artifact),
            data=rubric.data,REML=FALSE)
  
  ## do backwards elimination
  tmp.back_elim <- fitLMER.fnc(tmp,set.REML.FALSE = TRUE,log.file.name = FALSE)
  
  ## check to see if the raters are significantly different from one another
  tmp.single_intercept <- update(tmp.back_elim, . ~ . + 1 - as.factor(Rater))
  pval <- anova(tmp.single_intercept,tmp.back_elim)$"Pr(>Chisq)"[2]
  
  ## choose the best model
  if (pval<=0.05) {
    tmp_final <- tmp.back_elim
  } else {
    tmp_final <- tmp.single_intercept
  }

  ## and add to list...
  model.formula.alldata[[i]] <- formula(tmp_final)
  
}

## see what "final models" we got...
model.formula.alldata
```

### 2(c)(iii): Trying interactions and new random effects for the seven rubric specific models using all the data

Now we see there are some differences among the models: For `InitEDA`, `RsrchQ` and `SelMeth`, the models are just the simple random-intercept models.  For the other four, the models are a little more complex.  We should examine each of these 4 models to see (a) if the fixed effects make sense to us; and (2) if there are any interactions or additional random effects to consider. 

I will illustrate a bit of this with the model for `SelMeth`, and leave the rest to you.

```{r error=TRUE}


## refit the model and check on the t-statistics -- do all the variables matter?
fla <- formula(model.formula.alldata[["SelMeth"]])
tmp <- lmer(fla,data=tall.nonmissing[tall.nonmissing$Rubric=="SelMeth",])
round(summary(tmp)$coef,2)  ## fixed effects and their t-values

## apparently they do.

## now check to make sure we really need "Rater" as a factor...
tmp.single_intercept <- update(tmp, . ~ . + 1 - as.factor(Rater))
anova(tmp.single_intercept,tmp)

## looks like we do, so we keep "tmp" as our best model so far...

## now let's check for fixed-effect interactions... Since only Rater and Semester
## are involved, we only need to examine Rater*Semester

tmp.fixed_interactions <- update(tmp, . ~ . + as.factor(Rater)*Semester - Semester)
## I've specified the model so that I can see (a) a different intercept for each
## rater, and (b) a different semester effect for each rater.

anova(tmp,tmp.fixed_interactions)

## Looks like the fixed-effect interactions are not needed; again we keep 
## "tmp" as our best model so far... 

## Finally we check for random effects.  We should only add random effects that 
## are also present as fixed effects.  This means, for this model, we should try
## (Rater|Artifact) and (Semester|Artifact).

## I will show how to test these with exactRLRT()...

## Testing (Semester|Artifact)...

m0 <- tmp                                       ## Null hypothesis
mA <- update(m0, . ~ . + (Semester|Artifact))   ## Alternative hypotheses
m  <- update(mA, . ~ . - (1|Artifact))          ## Model with only the new R.E.

exactRLRT(m0=m0,mA=mA,m=m)

## Many error messages!  But note what the first one, for model mA is: there are
## more random effects than there are observations in the data set!  As explained
## on Piazza, this means lmer() cannot fit a model.  Thus, the model
##
## as.numeric(Rating) ~ -1 + as.factor(Rater) + Semester + 
##                           (1 | Artifact) + (Semseter | Artifact)
##
## isn't even possible, so no testing is needed.  



## Testng (as.factor(Rater)|Artifact)

m0 <- tmp                                       ## Null hypothesis
mA <- update(m0, . ~ . + (as.factor(Rater)|Artifact))   ## Alternative hypotheses
m  <- update(mA, . ~ . - (1|Artifact))          ## Model with only the new R.E.

exactRLRT(m0=m0,mA=mA,m=m)

## Same thing happened!  Again, the model
##
## as.numeric(Rating) ~ -1 + as.factor(Rater) + Semester + 
##                           (1 | Artifact) + (as.factor(Rater) | Artifact)
##
## isn't even possible, so no testing is needed.


## Thus, we weren't able to add or take away anything from the model "tmp", 
## so this is our final model for SelMeth:

summary(tmp)

## You would need to do something similar with the other models that are not 
## just intercept-only models (i.e the models for CritDes and InterpRes)

## ALSO, please don't forget:  I am not giving interpretations to the model fits
## or coefficient estimates here. That is something I'm leaving for you, as you
## complete your analyses and write an IDMRAD paper for the college.

```

### 2(c)(iv): Trying to add fixed effects, interactions, and new random effects to the "combined" model Rating ~ 1 + (0 + Rubric|Artifact), using all the data.


Now we try something similar with the "combined" model suggested on p. 4 of the project assignment sheet.

```{r}
## Start with the "combined" intercept-only model...

comb.0 <- lmer(as.numeric(Rating) ~ 1 + (0 + Rubric | Artifact),
               data=tall.nonmissing)

summary(comb.0)

## R complains that we have a "boundary (singular) fit", i.e. the
## variance-covariance matrix for the random effects is singular 
## (not of full rank), or nearly singular.  
##
## What this typically means is that some of the random effects are highly
## correlated with one another.  We can see this in the "Random effects" 
## block of summary(comb.0):  
##
## * The random effects for VisOrg and TxtOrg seem highly correlated with 
##   each other and with everything except for the rand. effect for SelMeth
##
## * The random effects for InterpRes and InitEDA are highly correlated
## 
## * The random effects for RsrchQ and InterpRes are highly correlated
## 
## etc.
##
## In some ways we should not be surprised: these rubrics all represent 
## features of a good research report, and we would expect that if someone
## is good at one or two of these features, they are probably good at the 
## others.

## Although the random effects are highly correlated, we can still proceed with 
## our variable selection...

## Try adding fixed effects with no interactions...

comb.full <- update(comb.0, . ~ . + as.factor(Rater) + Semester + 
                      Sex + Repeated + Rubric)

summary(comb.full)
##
## It's interesting to note that comb.full is no longer a boundary (singular) 
## fit.  Adding the fixed effects changed the residuals enough that the 
## variance-covariance matrix for the random effects is no longer (nearly) 
## singular.

comb.back_elim <- fitLMER.fnc(comb.full, log.file.name = FALSE)

summary(comb.back_elim)

## The final model fit is a boundary fit again, but we will proceed to try
## interactions

comb.inter <- update(comb.back_elim, . ~ . + as.factor(Rater)*Semester*Rubric)

## This didn't quite converge, so we will try switching optimizers and increasing
## the number of iterations allowed...

ss <- getME(comb.inter,c("theta","fixef"))
comb.inter.u<- update(comb.inter,start=ss,
             control=lmerControl(optimizer="bobyqa",
                                  optCtrl=list(maxfun=2e5)))

## it takes a few seconds to fit, but at least we got a converged fit.
## again, boundary fit (near-singular random effects variance-covariance mtx)

summary(comb.inter.u)
## If you compare with summary(comb.inter) you will see that
## there wasn't much difference in the fitted values; we could 
## probably have just proceeded wth the model comb.inter.  But 
## since we have the converged model we will use it for fixed  
## effects selection

comb.inter_elim <- fitLMER.fnc(comb.inter.u, log.file.name = FALSE)

summary(comb.inter_elim)

## it's a little hard to compare summaries for such big models, so let's look 
## at the highlights:

formula(comb.inter.u)
formula(comb.inter_elim)
formula(comb.back_elim)

summary(comb.inter.u)$varcor
summary(comb.inter_elim)$varcor
summary(comb.back_elim)$varcor

anova(comb.back_elim,comb.inter_elim,comb.inter.u)
## the models are nested so we can use AIC, BIC or likelihod ratio (deviance) 
## tests... AIC and the LRT agree on comb.inter_elim; BIC likes the simpler
## comb.back_elim.

## Interestingly, comb.inter_elim adds a rater x rubric interaction to 
## the main-effects model comb.back_elim.  This suggests that the raters 
## do not all use the rubrics in the same way.
```


```{r fig.height=8}
## In addition to looking at the fixed effect coefficients in
## summary(comb.inter_elim)$coef, we could also see if there's
## a pattern in an appropriate facets plot

g <- ggplot(tall.nonmissing, aes(x=Rating)) + 
  geom_bar() + 
  facet_wrap( ~ Rubric + Rater, nrow=7)

g

## and it does look as if the 3 raters have different ways of scoring the 7 rubrics,
## so the interaction we found in comb.inter_elim makes sense.  (Clearly it is
## not the case that one rater is simply more harsh than another, or something 
## like that.  Can you describe the different patterns of scoring among the 3 
## raters?  Can you relate it to coefficient estimates in summary()$coef?)
```


```{r error=TRUE}
## Finally, we consider adding random effects to what seems like the 
## best model so far, comb.inter_elim...

## The fixed-effects terms we have to work with are:
## 
## as.factor(Rater) 
## Semester
## as.factor(Rater):Rubric
##
## We want to add each of these *without* a random intercept, to preserve the 
## structure of the model (separate random interepts for each rubric)
##
## In all cases, there is more than one random effect to test (3 for raters, 
## 2 for semesters, 7 for rubrics, and 21 for the interaction).  Since exactRLRT()
## can only test single random effects, we can't use it.  Instead we inspect AIC 
## andBIC from anova() tables for these...

## Fitting some of these models produces various errors and warnings; I am not
## going to worry about them too much, in order to get an idea of what random
## effects I may want...

m0 <- comb.inter_elim 
mA <- lmer(as.numeric(Rating) ~ (0 + Rubric | Artifact) + 
             (0 + as.factor(Rater) | Artifact) + as.factor(Rater) + 
             Semester + Rubric + as.factor(Rater):Rubric, data=tall.nonmissing)

anova(m0,mA)
## 
## AIC and BIC both like including (0 + as.factor(Rater) | Artifact) in the model

m0 <- comb.inter_elim 
mA <- lmer(as.numeric(Rating) ~ (0 + Rubric | Artifact) + 
             (0 + Semester | Artifact) + as.factor(Rater) + 
             Semester + Rubric + as.factor(Rater):Rubric, data=tall.nonmissing)

anova(m0,mA)
##
## AIC and BIC do not like (0 + Semester | Artifact) in the model...


m0 <- comb.inter_elim 
mA <- lmer(as.numeric(Rating) ~ (0 + Rubric | Artifact) + 
             (0 + as.factor(Rater) | Artifact) + 
             (0 + as.factor(Rater):Rubric | Artifact) + as.factor(Rater) + 
             Semester + Rubric + as.factor(Rater):Rubric, data=tall.nonmissing)

## anova(m0,mA)     -- Not needed!
## 
## There are not enough observations to fit mA here, so we need not do any 
## formal model comparison...

## So, to summarize, the "final" model appears to be

comb.final <- lmer(as.numeric(Rating) ~ (0 + Rubric | Artifact) + 
             (0 + as.factor(Rater) | Artifact) + as.factor(Rater) + 
             Semester + Rubric + as.factor(Rater):Rubric, data=tall.nonmissing)

formula(comb.final)
summary(comb.final)$varcor
summary(comb.final)$coef

## if we accept comb.final as our final model, we can interpret the pieces as 
## follows:
##
## (0 + as.factor(Rater) | Artifact) + as.factor(Rater)
##   * There is a kind of Rater x Artifact interaction: each Rater's
##     rating on each Artifact differs from what we would expect (from the 
##     fixed effects alone) by a small random effect that depends on the Artifact  
##
## Rubric + as.factor(Rater) + as.factor(Rater):Rubric
##   * There is a Rater x Rubric interaction: each Rater uses each
##     Rubric in a way that is not like, or even parallel to, other rater's 
##     Rubric usage.  (we saw that in the facets plot above also).
##
## (0 + Rubric | Artifact) + Rubric
##   * There is a kind of Rubric x Artifact interaction: There are 
##     different average scores on each rubric, but the rubric averages also
##     vary a bit from one Artifact to the next, by a small random effect that
##     depends on Artifact

## In all of this, the fact that Rubric scores depend on Artifact (that is, 
## there is a kind of Rubric x Artifact interaction) is what we might expect: 
## the artifacts aren't all of equal quality on each rubric, and so we should 
## expect the average scores on each Rubric to vary from one Artifact to the next.
##
## More troubling are the Rater x Rubric interaction and the "kind of" 
## Rater x Artifact interaction.  The Rater x Rubric interaction suggests 
## that the Raters are not all interpreting the Rubrics in the same way.  The 
## "kind of" Rater x Artifact interaction suggests that the Raters are not 
## interpreting the evidence in the artifacts in the same way.  These 
## interactions suggest that perhaps the raters should be trained more, to
## make the raters' ratings more similar to each other.


```


## 2(d) Is there anything else interesting to say about this data?

I leave this entirely up to you.  By now you are pretty familiar with the data (probably sick of it actually!), so you may already have some ideas here.

Here are some things you could try (you don't need to try all of them, and there may be other, better things to try instead.  Just try to come up with something that is interesting and useful to the college):

* Do some additional EDA.  We really haven't tried to do a very complete EDA on this data yet.

* Are there better ways to illustrate or interpret the models we have fitted?

* Are there ways (through EDA, plots of model fits, etc.) to understand the differences in the models fitted to the data from the 13 common items, vs fitting to all the data?

* Is there any way to justify, based on the data and on where the data come from, imputing the missing Sex value as "M" or "F"?

* For some of the models, residual diagnostic plots might be interesting.

* What other models could you fit?  More multi-level models?  Some ordinary regression models?  glm's?  etc...

* etc.

I look forward to anything interesting that you can come up with!


