---
title: "Homework 09 Solutions"
date: "2022-11-15"
output: 
  pdf_document:
    keep_tex: true
    fig_caption: true
header-includes:
  - \usepackage{longtable}
  - \usepackage{booktabs}
  - \usepackage{dcolumn}
  - \usepackage{amsmath}
  - \usepackage{amsthm}
  - \usepackage{hyperref}
  - \usepackage{xcolor}
  - \DeclareMathOperator{\Cov}{Cov}
  - \DeclareMathOperator{\Corr}{Corr}
  - \DeclareMathOperator{\Var}{Var}
  - \DeclareMathOperator{\upmodels}{\perp \!\!\! \perp}
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, include = TRUE, message = FALSE)
# library(knitr)
```

```{r}
library(arm)
library(ggplot2)
library(HLMdiag)    ## easy extraction of MLM residuals
library(gridExtra)  ## For the grid.arrange() function

library(texreg)  ## for pretty regression tables
```

\def\red#1{\textcolor{red}{#1}}
\def\blue#1{\textcolor{blue}{#1}}
\let\bar\overline
\let\ul\underline

\def\theenumi{\roman{enumi}}

# Problem 1.

Explore the data, looking for usual or unusual relationships, usual or unusual data values, etc.  Perform any data cleaning, transformations, etc. that may be needed; justify your actions.  Submit a brief but interesting summary, with justifications.  

\it

This summary is a bit long, but hopefully it is still clear...
\label{p:1}

I begin by reading in the data\ldots

```{r}
data <- read.csv("ratings.csv",header=T)

str(data)
```


The following tables show that the researchers really did produce a balanced design:
```{r}
table(data$Harmony)
table(data$Instrument)
table(data$Voice)

with(data, table(paste(Harmony,Voice,Instrument)))
```

While exploring the distribution of the response variables we find that, although they are supposed to be integer ratings from 1 to 10, in fact we find some unusual values for them:
```{r }
##
## Counts of Classical responses:
rbind(count=table(data$Classical))

##
## Counts of Popular responses:
rbind(count=table(data$Popular))
```
We see
\begin{itemize}
\item A Classical rating of 0 was given 8 times, and a Popular rating of 0 was given 25 times
\item There are single ratings of 3.5, 4.2,4.6, 6.8, 9.5 and 19 in each response variable
\end{itemize}
I decided not to recode the ratings of 0 or the non-integer ratings between 1 and 10.  It is reasonable to guess that they reflect the intended responses of listeners who may not have entirely understood the "rate 1 to 10" instructions.

I will recode the 19's to 10's, for the following reasons:
\begin{itemize}
\item The Classical rating is paired with a modest Popular rating, and vice versa, so it's likely that the respondents intended to respond with 10's for these ratings:
\end{itemize}
```{r}

## The code is a little complicated because of the way R treats NA's in expressions. 
##
## Basically, if any part of an expression evaluates to NA, then the whole expression
## evaluates to NA.  In the code below, I first filter out the NA's with
## !is.na(data$Classical) for example, and then test the condition that I'm 
## interested in (data$Classical=="19") on the remaining rows.  Otherwise I would 
## get a whole set of rows with NA's instead of selecting the row with the "19" in it.
##
## Similarly for the "Popular" code.

##
## The Classical rating of 19:
c19 <-data$Subject[!is.na(data$Classical) & data$Classical=="19"]
x <- data[data$Subject==c19,c("Subject", "Classical","Popular")]
x[!is.na(x$Classical) & x$Classical=="19",]

##
## The Popular rating of 19:
p19 <- data$Subject[!is.na(data$Popular) & data$Popular=="19"]
x <- data[data$Subject==p19,c("Subject", "Classical","Popular")]
x[!is.na(x$Popular) & x$Popular=="19",]
```
\begin{itemize}
\item The "9" key is next to the "0" key on a computer keyboard, so it seems likely that the 19's are mis-typed 10's.
\end{itemize}

Next we see that several variables have a significant portion of missing values:
```{r}
##
## Count and Percent of Missing Values:
x <- data.frame(count=sapply(data,function(x) sum(is.na(x))),
      percent=sapply(data,function(x) round(100*mean(is.na(x)),2)))
x[x$percent!=0,]
```

\begin{itemize}

\item The missing values in the response variables Classical and Popular cannot (and should not) be imputed, so I will delete those rows.  Further exploration (not shown here) shows that the same rows with missing Classical values are missing the Popular values, and vice versa.
\item The variables X1stInstr and X2ndInstr have 60\% and 80\% missing values, respectively, which makes them hard to use and hard to impute so I will delete those variables.
\end{itemize}

Taking into account all of the above, our full data set can be constructed as follows
```{r}
fulldata <- data[,!(names(data) %in% c("X1stInstr","X2ndInstr"))]
fulldata <- fulldata[!is.na(fulldata$Classical),]
fulldata <- fulldata[!is.na(fulldata$Popular),]
fulldata$Classical[fulldata$Classical=="19"]<-"10"
fulldata$Popular[fulldata$Popular=="19"]<-"10"
```

From the tables and counts below we see that we have lost 27 rows (deleting rows with missing responses) and 2 columns (deleting the two variables X1stInstr and X2ndInstr), and we no longer have a perfectly balanced experiment.   
```{r}
dim(data)
dim(fulldata)
table(fulldata$Harmony)
table(fulldata$Instrument)
table(fulldata$Voice)

with(fulldata, table(paste(Harmony,Voice,Instrument)))

table(fulldata$Classical)

table(fulldata$Popular)

length(unique(fulldata$Subject))
```

The "fulldata" set still has many missing values. For model selection, when I need to compare models on the same data set, it may be useful to have a data set with no missing values.  A simple approach is to construct a reduced data set by removing rows with missing values:
```{r}
nomissdata <- fulldata[apply(fulldata,1,function(x) !any(is.na(x))),]
dim(fulldata)
dim(nomissdata)

table(nomissdata$Harmony)
table(nomissdata$Instrument)
table(nomissdata$Voice)

length(unique(nomissdata$Subject))

with(nomissdata, table(paste(Harmony,Voice,Instrument)))

table(nomissdata$Classical)

table(nomissdata$Popular)


```
We have lost nearly 1000 observations and we are down to 43 subjects from the original 70, but the data is still fairly balanced among the experimental conditions, and most of the 43 subjects saw all 36 combinations of Harmony, Instrument and Voice.
\begin{itemize}
\item I will use ``fulldata'' for Problem \#2 below, since it only involves response variables and design variables, and there is no missingness in these variables in the ``fulldata'' data set.
\item I will begin with ``nomissdata'' for Problems \#3 and \#4 below but if it becomes obvious that some variables with missing values are not useful in variable selection, I may be able to delete those variables and have a somewhat larger data set for model comparisons.
\end{itemize}

Finally it is worthwhile to look at histograms of the variables (a) to see if the distributions are similar in "fulldata" and "nomissdata" and (b) to check for unusual distributions of any kind.

In the plots on the next two pages, we see that

\begin{itemize}

\item The distributions of the numeric variables are broadly similar across the two data sets

\item Most variables do not seem to need transformations, either because they are really discrete/categorical variables (e.g. all the variables on a scale from 0 to 5) or they are not sufficiently skewed to warrant transformation.

\item Possible exceptions to the previous point are the first three variables plotted, Selfdeclare, OMSI and X16.minus.17.  Selfdeclare looks like another categorical variable so I would probably leave it alone.  
\begin{itemize}

\item  OMSI looks continuous and right-skewed so I could consider a log or similar transformation.  

\item X16.minus.17 could be shifted and then transformed.
\end{itemize}
but it's not clear how meaningful the transformed variables would be, and so I'm inclined to leave these two variables alone as well.

\end{itemize}

\clearpage

```{r fig.width=8,fig.height=10}
par(mfrow=c(5,4))
for (i in names(fulldata)[-c(1:5,24)]) { ## select only the numeric variables
  hist(as.numeric(fulldata[,i]),xlab=i,main="",cex=0.5)
}
```

```{r fig.width=8,fig.height=10}
par(mfrow=c(5,4))
for (i in names(nomissdata)[-c(1:5,24)]) {  ## select only the numeric variables
  hist(as.numeric(nomissdata[,i]),xlab=i,main="",cex=0.5)
}
```

\rm


In your work below, \texttt{Classical} should not appear in any models for \texttt{Popular}, nor vice-versa.  When evaluating models, please remember to look at appropriate plots as well as test statistics.

# Problem 2.

\ul{\textit{The three main experimental factors.}}

## Problem 2(a).

Examine the influence of the three main experimental 
factors (Instrument, Harmony \& Voice) on Classical ratings, using 
conventional linear models and/or analysis of variance models.  
Comment briefly on your findings, providing suitable brief evidence for each
result. \textit{Hint: Since there are only three factors here, it is worth considering interactions of all orders.}
\label{p:2a}

\it

```{r}


## Build some basic models to check whether we need 3-way or 2-way interactions 

lm.all.interactions <- lm(Classical ~ Instrument * Harmony * Voice, data=fulldata)

lm.no.threeway <- update(lm.all.interactions, . ~ . - Instrument : Harmony : Voice)

lm.no.interactions <- lm(Classical ~ Instrument + Harmony + Voice, data=fulldata)

## Here's the "full model" will all interactions.   Not too much is significant

display(lm.all.interactions)

## Let's compare the models with F tests...

anova(lm(Classical ~ 1, data=fulldata),lm.no.interactions,lm.no.threeway,lm.all.interactions)

## From this it seems like just the additive model with no interactions is best, 
## but let's at least check the pairwise two-way interactions... 

anova(lm.no.interactions, update(lm.no.interactions, . ~ . + Instrument:Harmony))

anova(lm.no.interactions, update(lm.no.interactions, . ~ . + Instrument:Voice))

anova(lm.no.interactions, update(lm.no.interactions, . ~ . + Voice:Harmony))

## It seems like having Voice:Harmony is a good idea; we can summarize as follows:


lm.vh <-  update(lm.no.interactions, . ~ . + Voice:Harmony)

anova(lm(Classical ~ 1, data=fulldata),lm.no.interactions,lm.vh,lm.no.threeway,lm.all.interactions)

## We could also look at AIC and BIC, where we don't need to worry 
## so much about nested models....

lm.0 <- lm(Classical ~ 1, data=fulldata)
lm.iv <-  update(lm.no.interactions, . ~ . + Instrument:Voice)
lm.ih <-  update(lm.no.interactions, . ~ . + Instrument:Harmony)

IC <- function(M) {
    x <-c(AIC=AIC(M),BIC=BIC(M))
    return(x)
}

print(x <- rbind(lm.0=IC(lm.0),
      lm.no.interactions=IC(lm.no.interactions),
      lm.iv=IC(lm.iv),
      lm.ih=IC(lm.ih),
      lm.vh=IC(lm.vh),
      lm.no.threway=IC(lm.no.threeway),
      lm.all.interactions=IC(lm.all.interactions))
)

x <- as.data.frame(x)

x[x$AIC==min(x$AIC),]

x[x$BIC==min(x$BIC),]

## So the result is much the same as for the F-test comparisons. 
## The best BIC model is the no-interactions model, and the best 
## AIC model adds the Voice:Harmony interaction to that model.

```
There is some evidence in favor of keeping the Voice:Harmony interaction (AIC, and ANOVA comparing this model with the no-interactions model), but also some evidence against (BIC, and ANOVA comparing the model with all two-way interactions with the no-interactions model).  

One reason for keeping the interaction is that it give you more to talk about with Dr. Jimenez.  A reason for not keeping it is that the evidence is that, in the face of ambiguous evidence, simpler is better.

In what follows I will drop the Voice:Harmony interaction term, but you could argue it either way. 

\rm

## Problem 2(b).

Since we have approximately 36 ratings from each 
participant, we can expand the model in part~2(a) by including a random intercept for each participant.  Such a model is sometimes called a "repeated measures" model. 

### 2(b)(i) 

Write this expanded model in mathematical terms as a multi-level
model, using notation we have used in class.  Be sure to carefully identify what are the groups, and what are the individual observations.
\label{p:2b.i}


\it

In this case the groups $j$ are the listeners in the experiment and the individual observations $i$ are the individual ratings made by each listener.  The two-level model is then
\def\c#1{\mbox{\scriptsize\it #1}}
\begin{eqnarray*}
\mbox{Classical}_i & = & \alpha_{0j[i]} + 
\beta_{\c{guitar}} \cdot 1_{\{\c{Instrument}_i=\c{guitar}\}} + 
\beta_{\c{piano}} \cdot 1_{\{\c{Instrument}_i=\c{piano}\}} + 
\beta_{\c{string}} \cdot 1_{\{\c{Instrument}_i=\c{string}\}} + \\ & &  
\beta_{\c{contrary}} \cdot 1_{\{\c{Voice}_i=\c{contracy}\}} + 
\beta_{\c{par3rd}} \cdot 1_{\{\c{Voice}_i=\c{par3rd}\}} + 
\beta_{\c{par5th}} \cdot 1_{\{\c{Voice}_i=\c{par5th}\}} + \\ & & 
\beta_{\c{I-IV-V}} \cdot 1_{\{\c{Harmony}_i=\c{I-IV-V}\}} + 
\beta_{\c{I-V-IV}} \cdot 1_{\{\c{Harmony}_i=\c{I-V-IV}\}} + \\ & & 
\beta_{\c{I-V-VI}} \cdot 1_{\{\c{Harmony}_i=\c{I-V-VI}\}} + 
\beta_{\c{IV-I-V}} \cdot 1_{\{\c{Harmony}_i=\c{IV-I-V}\}} + \epsilon_i ~,~~ \epsilon_i \stackrel{iid}{\sim} N(0,\sigma^2)\\
\alpha_{0j} & = & \beta_0 + \eta_{0j} ~,~~ \eta_j \stackrel{iid}{\sim} N(0,\tau_0^2)
\end{eqnarray*}
where $1_{\{c{condition}\}}$ are dummy variables that equal 1 if $\{\mbox{\it condition}\}$ is true, and equal 0 if $\{\mbox{\it condition}\}$ is false.  Note that because of collinearity with the intercept, the first level of every factor will be dropped by lmer().

\bf Note:\it{} All the fixed-effect terms for the experimental factors have to be in the level-1 model since they are all level-1 covariates (that is, they vary with $i$ rather than $j$, and the level-2 model can only deal with terms that vary with $j$).

\rm


### 2(b)(ii) 

Re-examine the influence of the three main experimental factors 
(Instrument, Harmony \& Voice) on Classical ratings (that is, repeat part 2(a)), using models that always include the random intercept for participants (Please make sure the models are fitted with REML=FALSE (i.e. maximum likelihood) to ensure that the model comparisons make sense.). 
\label{p:2b.ii}

\it

```{r}

## Again, we build some basic models to check whether we need 3-way or 2-way interactions.

## Note that in the data frame, "Classical" and "Popular" are, for some reason, string
## variables.  lm() was able to deal with this by converting the strings to their 
## obvious numeric equivalents, but lmer() can't deal. So we have to convert them, before
## moving forward.  We'll go ahead and convert them in both "fulldata" and"nomissdata"

fulldata$Classical <- as.numeric(fulldata$Classical)
fulldata$Popular <- as.numeric(fulldata$Popular)
nomissdata$Classical <- as.numeric(nomissdata$Classical)
nomissdata$Popular <- as.numeric(nomissdata$Popular)

lmer.all.interactions <- lmer(Classical ~ Instrument * Harmony * Voice + (1|Subject),
                              data=fulldata)

lmer.no.threeway <- update(lmer.all.interactions, . ~ . - Instrument : Harmony : Voice)

lmer.no.interactions <- lmer(Classical ~ Instrument + Harmony + Voice + (1|Subject),
                             data=fulldata)

display(lmer.all.interactions)

## So here, \hat\tau = 1.30 and \hat\sigma = 1.87, whereas for lm.all.interactions,
## \hat\sigma = 2.28.  It's interesting (and not surprising!! why??) to note that
## sqrt(1.30^2 + 1.87^2) = 2.277477: Taking away the square roots, the two variance 
## components from the lmer() model add up to the residual variance in the lm() model.

lmer.0 <- lmer(Classical ~ 1 + (1|Subject),data=fulldata)

lmer.iv <- update(lmer.no.interactions, . ~ . + Instrument:Voice)
lmer.ih <- update(lmer.no.interactions, . ~ . + Instrument:Harmony)
lmer.vh <- update(lmer.no.interactions, . ~ . + Voice:Harmony)

## For AIC and BIC we just need to make sure that the lmer() model is fitted 
## with REML=F

IC <- function(M) {
    M.ml <- update(M,REML=F)  ## ensure the model is fitted with REML=F
    x <-c(AIC=AIC(M.ml),BIC=BIC(M.ml))
    return(x)
}


print(x <- rbind(lmer.0=IC(lmer.0),
      lmer.no.interactions=IC(lmer.no.interactions),
      lmer.iv=IC(lmer.iv),
      lmer.ih=IC(lmer.ih),
      lmer.vh=IC(lmer.vh),
      lmer.no.threway=IC(lmer.no.threeway),
      lmer.all.interactions=IC(lmer.all.interactions))
)

x <- as.data.frame(x)

x[x$AIC==min(x$AIC),]

x[x$BIC==min(x$BIC),]

```
There are two interesting things to note here;
\begin{itemize}
\item Each lmer() model fits better than the corresponding lm() model, in the sense that the AIC and BIC are always lower for the lmer() model than the lm() model.  This already suggests that the random effect is useful (though we should check this with BIC also).
\item Once again, AIC picks the model with Voice:Harmony interaction, and BIC picks the no-interaction model.
\end{itemize}


\rm

### 2(b)(iii) 

Is the random effect needed?  Justify 
your answer with evidence.  You may (or may not) find valid evidence in fitted models, in statistical tests, and/or in fit indices.
\label{p:2b.iii}

\it


```{r}

## We begin by defining AIC, BIC, and DIC so that they apply to both 
## lmer() models with mle's and to lm() models.

## From now on I'll assume that AIC, BIC and DIC are evaluated on lmer() 
## models with REML=F

## AIC -- no change needed

## BIC -- no change needed

## DIc -- need to check to see if the model is "lm" or "lmer"...


DIC <- function(M) {
         if(class(M)=="lm") { x <- AIC(M) }
         else { x <- unname(extractDIC(M)) }   
         return(x)
       }


IC <- function(M) {
    if (class(M)=="lm") { M.ml <- M }  ## don't have to refit lm() models.
    else { M.ml <- update(M,REML=F) }  ## ensure that lmer() is fitted with REML=F
    x <-c(AIC=AIC(M.ml),BIC=BIC(M.ml),DIC=DIC(M.ml))
    return(x)
}

## I'm going to go ahead and calculate all the IC's for all the models, 
## in case it is interesting to look at them later

x1 <- rbind(lm.0=IC(lm.0),
      lm.no.interactions=IC(lm.no.interactions),
      lm.iv=IC(lm.iv),
      lm.ih=IC(lm.ih),
      lm.vh=IC(lm.vh),
      lm.no.threway=IC(lm.no.threeway),
      lm.all.interactions=IC(lm.all.interactions))

x2 <- rbind(lmer.0=IC(lmer.0),
      lmer.no.interactions=IC(lmer.no.interactions),
      lmer.iv=IC(lmer.iv),
      lmer.ih=IC(lmer.ih),
      lmer.vh=IC(lmer.vh),
      lmer.no.threway=IC(lmer.no.threeway),
      lmer.all.interactions=IC(lmer.all.interactions))

## All I am really interested in is the difference in DIC values between 
## the lm() models (which do not have a random effect) and the lmer() models
## (which have a random intercept for each lister).

cbind("DIC(lm)-DIC(lmer)" = x1[,3] - x2[,3])


```

The differences in DIC values between models without a random intercept (the lm() models)
 and models with a random intercept (the lmer() models) is never less than 400, and always favors the lmer() model with the random intercept. For the two models that seem to be working best, the no-interactions model and the model that adds Voice:Harmony to that model, the difference is around 800.  
 
Since we know a difference of 3 is interesting and a difference of 10 is really something to pay attention to, clearly the random intercept for each person is important to have in the model, regardless of which fixed effects we use.
\rm


## Problem 2(c).

The random intercept in a repeated measures model 
can account for "personal biases" in ratings: for example, perhaps person A is more
inclined to rate everything as classical, and person B is more inclined 
to rate everything as popular.  This can be accounted for by the random 
intercept.  Alternatively, perhaps personal biases vary with the type 
of instrument, type of harmony, and/or type of voice leading.  For 
example, perhaps people vary in the degree to which they are inclined 
to call music played by a string quartet "classical".  This suggests, 
e.g., a random effect of the form \texttt{(Instrument|Subject)} or \texttt{(0 + Instrument|Subject)} (and similarly for \verb|Harmony| and \verb|Voice|).  One could argue for a similar random effect for each person/harmony combination, and for each person/voice 
leading combination.   

### 2(c)(i)

Determine whether the model in 2(b)(ii) is improved by adding one or more of these three new random effect terms.  Find the best combination of these terms (A good rule of thumb is to only include a term like (Voice$|$Subject) in the model if Voice is already a fixed effect.) and provide suitable evidence for your answer.
\label{p:2c.i}

\it

So far the two models 
\begin{itemize}
\item \verb|lmer.no.interactions|, with formula \verb/Classical ~ Instrument + Harmony + Voice + (1 | Subject)/
\item \verb|lm.vh|, with formula \verb/Classical ~ Instrument + Harmony + Voice + (1 | Subject) + Harmony:Voice/
\end{itemize}
seem to be the best from 2(b)(ii).  I will show what we get for each of these two models.

There are two things to note before I begin:
\begin{itemize}
\item When I am considering more than one random effect, I can model without correlations, for example like this:
\begin{verbatim}
Classical ~ Inst + Voic + Harm + (1|Subject) + (0+Inst|Subject) + (0+Voic|Subject) 
\end{verbatim}
or with correlations, for example like this:
\begin{verbatim}
Classical ~ Inst + Voic + Harm + (1 + Inst + Voic | Subject)
\end{verbatim}
I am going to first try to model \bf with\it{} correlations.  But it turns out that the models are a little more interpretable if we force some correlations to be zero, so I will show final models like that also.
\item In addition to the original models, there are seven modified models to consider (3 models with one additional random effect, 3 models with two additional random effects, and 1 model with all three random effects).
\end{itemize}

#### Let's first look at `lmer.no.interactions`:

```{r}

formula(lmer.no.interactions)

lmer.ni <- lmer.no.interactions ## just to have a shorter name to type!!

lmer.ni.i <- update(lmer.ni, . ~ . - (1|Subject) + (1 + Instrument|Subject))

lmer.ni.v <- update(lmer.ni, . ~ . - (1|Subject) + (1 + Voice|Subject))

lmer.ni.h <- update(lmer.ni, . ~ . - (1|Subject) + (1 + Harmony|Subject))

lmer.ni.iv <- update(lmer.ni, . ~ . - (1|Subject) + (1 + Instrument + Voice|Subject))

lmer.ni.ih <- update(lmer.ni, . ~ . - (1|Subject) + (1 + Instrument + Harmony|Subject))

lmer.ni.vh <- update(lmer.ni, . ~ . - (1|Subject) + (1 + Voice + Harmony|Subject))

lmer.ni.ivh <- update(lmer.ni, . ~ . - (1|Subject) + (1 + Instrument + Voice + Harmony|Subject))

DIC.ml <- function(M) DIC(update(M,REML=F))

print(x <- rbind(lmer.ni=c(DIC=DIC.ml(lmer.ni)),
      lmer.ni.i=c(DIC=DIC.ml(lmer.ni.i)),
      lmer.ni.v=c(DIC=DIC.ml(lmer.ni.v)),
      lmer.ni.h=c(DIC=DIC.ml(lmer.ni.h)),
      lmer.ni.iv=c(DIC=DIC.ml(lmer.ni.iv)),
      lmer.ni.ih=c(DIC=DIC.ml(lmer.ni.ih)),
      lmer.ni.vh=c(DIC=DIC.ml(lmer.ni.vh)),
      lmer.ni.ivh=c(DIC=DIC.ml(lmer.ni.ivh)))
)

```
Three observations are important to make here. 

First, when we specify a term like (Harmony|Subject), we get a random coefficient for \bf each level\it{} of the discrete Harmony variable, as can be seen in this summary of the model \verb|lmer.ni.h|:
```{r}
display(lmer.ni.h)
```
Second, six of the models complain about "boundary (singular) fit". As we discussed in class, this can mean either (a) a $\tau^2$ is essentially 0; or (b) a $\rho$ is essentially $-1$or +1 (more generally this complaint happens when the variance covariance matrix of the random vector of $\eta$ is not of full rank).  For example, in the summary of \verb|lmer.ni.h| above, we can see correlations of +1 and $-$1.  

Correlations of +1 and $-1$mean that there are really fewer random effects than we have tried to specify.  In the case of \verb|lmer.ni.h| we see that the random effect for \verb|Harmony = I-V-IV| is identical to the random intercept (since their correlation is $-1$), and the random effect for \verb|Harmony = I-V-VI| is exactly $-1$ times the random effect for \verb|Harmony = IV-I-V| (since their correlation is $-1$).

These correlations aren't necessarily bad (especially because there were no convergence warnings); they are just telling us we could have a simpler model if some of the coefficients shared the same random component $\eta$.  It is not so easy to implement a constraint like ``sharing $\eta$'s'' in lmer(), so we won't try.  We'll just observe that the model seems to be a bit over-specified.

(It is conceptually easier to implement a ``sharing $\eta$'s'' constraint when building a Bayesian version of the model, and maybe we will take a look at this later in the semester).

Finally, there is a clear winner in terms of DIC.  Remarkably enough, it is the most complicated model.  Here is a summary of the model:

```{r}
display(lmer.ni.ivh)

## There is a warning here about "maxfun" being less than the recommended
## number.  "maxfun" is the number of iterations needed for the non-linear
## maximization (that we talked about briefly in class).  If we were also
## getting a convergence warning, I would want to increase maxfun (you can
## figure out how with a google search). But since there doesn't seem to be 
## a convergence problem, I'm not going to worry about maxfun.
```
This model also has a singular variance-covariance matrix for the $\eta$'s, but it is not because any correlation is +1 or $-$1.  Instead, the particular combination of correlations makes the variance-covariance matrix less than full rank. That makes the model a little harder to interpret, and I will show an alternative below that is a little easier to interpret.

#### Now let's look at `lmer.ni.vh`:


```{r}

formula(lmer.vh)

lmer.vh.i <- update(lmer.vh, . ~ . - (1|Subject) + (1 + Instrument|Subject))

lmer.vh.v <- update(lmer.vh, . ~ . - (1|Subject) + (1 + Voice|Subject))

lmer.vh.h <- update(lmer.vh, . ~ . - (1|Subject) + (1 + Harmony|Subject))

lmer.vh.iv <- update(lmer.vh, . ~ . - (1|Subject) + (1 + Instrument + Voice|Subject))

lmer.vh.ih <- update(lmer.vh, . ~ . - (1|Subject) + (1 + Instrument + Harmony|Subject))

lmer.vh.vh <- update(lmer.vh, . ~ . - (1|Subject) + (1 + Voice + Harmony|Subject))

lmer.vh.ivh <- update(lmer.vh, . ~ . - (1|Subject) + (1 + Instrument + Voice + Harmony|Subject))

print(x <- rbind(lmer.vh=c(DIC=DIC.ml(lmer.vh)),
      lmer.vh.i=c(DIC=DIC.ml(lmer.vh.i)),
      lmer.vh.v=c(DIC=DIC.ml(lmer.vh.v)),
      lmer.vh.h=c(DIC=DIC.ml(lmer.vh.h)),
      lmer.vh.iv=c(DIC=DIC.ml(lmer.vh.iv)),
      lmer.vh.ih=c(DIC=DIC.ml(lmer.vh.ih)),
      lmer.vh.vh=c(DIC=DIC.ml(lmer.vh.vh)),
      lmer.vh.ivh=c(DIC=DIC.ml(lmer.vh.ivh)))
)

```
All of the comments that applied to the models we built above that elaborated on the `lmer.ni` model also  apply to these models as well.

And once again, the most complex model, \verb|lmer.vh.ivh|, is preferred by DIC.  Let's take a quick look at it:
```{r}
display(lmer.vh.ivh)
```
As before, the singularity in this model is not due to a specific correlation of +1 or $-1$, but rather due to the fact that the estimated variance-covariance matrix for the vector of $\eta$'s is less than full rank. 

#### Alternative models with correlations constrained to be zero.

The last thing I want to show is what these models look like if we force correlations between (groups of) random effects to be zero. 

```{r}
lmer.ni.ivh.c0 <- lmer(Classical ~ Instrument + Harmony + Voice + (1|Subject) + 
    (0+Instrument|Subject) + (0+Voice|Subject) + (0+Harmony|Subject),
    data=fulldata)

## or equivalently lmer.ni.ivh.c0 <- lmer(Classical ~ Instrument + Harmony + Voice + 
##                                (Instrument+Harmony+Voice||Subject),data=fulldata)
display(lmer.ni.ivh.c0)

lmer.vh.ivh.c0 <- lmer(Classical ~ Instrument + Harmony*Voice + 
                         (Instrument+Harmony+Voice||Subject),data=fulldata)

display(lmer.vh.ivh.c0)
```

These models do not fit as well as the models where the correlations are free (and the differences are pretty dramatic):
```{r}
x <- t(sapply(list(lmer.ni.ivh,lmer.ni.ivh.c0,lmer.vh.ivh,lmer.vh.ivh.c0),IC))
dimnames(x)[[1]] <-c("lmer.ni.ivh","lmer.ni.ivh.c0","lmer.vh.ivh","lmer.vh.ivh.c0")
x
```
However, the random effects seem to make a bit more sense.  For example, comparing the $\hat\tau$'s and $\hat\rho$'s for the \verb|lmer.ni.ivh| model vs the \verb|lmer.ni.ivh.c0| model,
```{r}
VarCorr(lmer.ni.ivh)
VarCorr(lmer.ni.ivh.c0)
```
we see that $\hat\sigma$ is almost the same for the two models, but instead of an overall intercept $\eta_{0j}$ (whose $\hat\tau_0$ is so small we could consider omitting it), we are getting a distinct $\eta$ for every level of every design variable, and the correlations between the $\eta$'s make some sense.  For example, referring to `lmer.ni.ivh.c0`, 

\begin{itemize}

\item Among the $\eta$'s for Instrument, $\hat\rho(\eta_{Instrumentguitar},\eta_{Instrumentstring})=-0.968$ .  This suggests that the listeners' personal tendency to call a piece classical when the instrument is strings is equal to and opposite their personal tendency to call it classical when the instrument is guitar.  On the other hand  $\hat\rho(\eta_{Instrumentguitar},\eta_{Instrumentpiano})=0.291$ and $\hat\rho(\eta_{Instrumentpiano},\eta_{Instrumentstring})=-0.282$, suggesting that any personal tendency to call something classical if it is piano is much more ambiguous than for the other two instrments.

\item Looking at the $\hat\rho$'s among the $\eta$'s for Voice, it seems like the tendency to call something classical doesn't change that much with the type of voice leading; all three $\eta$'s here are pretty highly correlated (estimated correlations of 0.919, 0.967 and 0.989) 

\item Looking at the $\hat\rho$'s among the $\eta$'s for Harmony, the personal tendency to call something classical is pretty much the same for the I-IV-V, I-V-IV and IV-I-V harmonies, but how listeners hear the I-V-VI harmony is somewhat different.

\end{itemize}
Very similar interpretations can be made for the $\hat\rho$'s and $\hat\tau$'s for the `lmer.vh.ivh.c0` model.  Further analysis using the ranef() function to extract estimated $\eta$ values would be useful to confirm these ideas.   

As modelers we would now have to decide (in consultation with Dr. Jimenez, of course) whether the greater interpretability is worth the substantially worse fits of these models.

In the remainder of this document I will focus on the model `lmer.ni.ivh.c0` (no fixed-effect interactions, and zero correlations between groups of $\eta$'s) because
\begin{itemize}
\item The results for fixed effects are similar to what we would get if we allowed all the $\rho$'s to be estimated, and the random effect correlations are more interpretable. 
\item lmer() runs a little faster on them, and we have a lot of models to fit below.
\end{itemize}
There is an argument for keeping the Voice:Harmony in the model, but it just makes the calculations, and especially writing out multilevel and hierarchical models. In addition, as we will see below, the model with the Voice:Harmony interaction is consistently selected by by AIC, and the model without this interaction is consistently selected by ABF I may return to this interaction later.

\rm



### 2(c)(ii)

Re-examine the influence of the three main experimental factors 
(Instrument, Harmony \& Voice) on Classical ratings, using the 
random effects you found in part (i). Comment briefly on 
your findings, providing suitable brief evidence for each result. In 
addition, comment on the sizes of the variances of whichever random effects you added to the model, 
with respect to each other and with respect to the estimated residual 
variance.
\label{p:2c.ii}

\it


This is basically going back to the kind of analysis we did in 2(a), but now with the random effects we have found in 2(c)(i). 

First we re-consider the fixed effects with fully-correlated $\eta$'s, `(Instrument + Voice + Harmony | Subject)`:

```{r}
## We do not consider
##  lmer.0 <- lmer(1 + (Instrument + Voice + Harmony | Subject), data=fulldata) 
## because it does not have fixed effects corresponding to the \eta's for 
## Instrument, Voice and Harmony.
##

lmer.ni.ivh <- lmer(Classical ~ Instrument + Voice + Harmony + 
                               (Instrument + Voice + Harmony | Subject), data=fulldata)

lmer.iv.ivh <- lmer(Classical ~ Instrument + Voice + Harmony + Instrument:Voice +  
                               (Instrument + Voice + Harmony | Subject), data=fulldata)

lmer.ih.ivh <- lmer(Classical ~ Instrument + Voice + Harmony + Instrument:Harmony + 
                               (Instrument + Voice + Harmony | Subject), data=fulldata)

lmer.vh.ivh <- lmer(Classical ~ Instrument + Voice + Harmony + Voice:Harmony + 
                               (Instrument + Voice + Harmony | Subject), data=fulldata)

## There is a convergence warning here, but the convergence measure 0.0088 is reasonably 
## close to the tolerance 0.002, so I won't worry too much about it. I could try other
## optimizers, or just increase maxfun, to see if I can get to convergence but I will 
## let this be for now.

lmer.no.threeway.ivh <- lmer(Classical ~ Instrument * Voice * Harmony - 
                           Instrument:Voice:Harmony +  
                               (Instrument + Voice + Harmony | Subject), data=fulldata)

lmer.all.interactions.ivh <- lmer(Classical ~ Instrument * Voice * Harmony + 
                               (Instrument + Voice + Harmony | Subject), data=fulldata)


print(x <- rbind(
      lmer.ni.ivh=IC(lmer.ni.ivh),
      lmer.iv.ivh=IC(lmer.iv.ivh),
      lmer.ih.ivh=IC(lmer.ih.ivh),
      lmer.vh.ivh=IC(lmer.vh.ivh),
      lmer.no.threeway.ivh=IC(lmer.no.threeway.ivh),
      lmer.all.interactions.ivh=IC(lmer.all.interactions.ivh))
)

x <- as.data.frame(x)
x[x$AIC==min(x$AIC),]
x[x$BIC==min(x$BIC),]

```
AIC likes the model `lmer.vh.ivh` with only the interaction Voice:Harmony, while BIC likes the simpler no-interactions model `lmer.ni.ivh`.

Now, let's try the same models forcing correlations beween groups of $\eta$'s to be zero, `(Instrument + Voice + Harmony || Subject)`:

```{r}

lmer.ni.ivh.c0 <- lmer(Classical ~ Instrument + Voice + Harmony + 
                               (Instrument + Voice + Harmony || Subject), data=fulldata)

lmer.iv.ivh.c0 <- lmer(Classical ~ Instrument + Voice + Harmony + Instrument:Voice +  
                               (Instrument + Voice + Harmony || Subject), data=fulldata)

lmer.ih.ivh.c0 <- lmer(Classical ~ Instrument + Voice + Harmony + Instrument:Harmony + 
                               (Instrument + Voice + Harmony || Subject), data=fulldata)

lmer.vh.ivh.c0 <- lmer(Classical ~ Instrument + Voice + Harmony + Voice:Harmony + 
                               (Instrument + Voice + Harmony || Subject), data=fulldata)

lmer.no.threeway.ivh.c0 <- lmer(Classical ~ Instrument * Voice * Harmony - 
                              Instrument:Voice:Harmony +  
                               (Instrument + Voice + Harmony || Subject), data=fulldata)

lmer.all.interactions.ivh.c0 <- lmer(Classical ~ Instrument * Voice * Harmony + 
                               (Instrument + Voice + Harmony || Subject), data=fulldata)

print(x <- rbind(
      lmer.ni.ivh.c0=IC(lmer.ni.ivh.c0),
      lmer.iv.ivh.c0=IC(lmer.iv.ivh.c0),
      lmer.ih.ivh.c0=IC(lmer.ih.ivh.c0),
      lmer.vh.ivh.c0=IC(lmer.vh.ivh.c0),
      lmer.no.threeway.ivh.c0=IC(lmer.no.threeway.ivh.c0),
      lmer.all.interactions.ivh.c0=IC(lmer.all.interactions.ivh.c0))
)

x <- as.data.frame(x)
x[x$AIC==min(x$AIC),]
x[x$BIC==min(x$BIC),]

```
For these models, AIC prefers `lmer.vh.ivh.c0`, while BIC prefers `lmer.ni.ivh.c0`.  These are the same fixed effects as found with the `(Instrument + Voice + Harmony | Subject)` random effects.    

For simplicity I will prefer the `lmer.ni.ivh.c0` model, but the Voice:Harmony interaction keeps coming up; it may be worth returning to the model(s) with the Voice:Harmony interaction later.

\rm



### 2(c)(iii)

Carefully write this model in mathematical terms as a hierarchical
linear model, using notation from class. Indicate values of estimated parameters (and their SE's, if available).
\label{p:2c.iii}

\it

I am taking 
```{r eval=F}
lmer.ni.ivh.c0 <- lmer(Classical ~ Instrument + Voice + Harmony + 
                               (Instrument + Voice + Harmony || Subject), data=fulldata)
```
to be the my preferred model, for now.

First, I will write this as a multilevel model (it will look very similar to the model for 2(b)(i)), and then I will re-write it as a hierarchical model.

I will omit the correlations between $\eta$'s below when I write the multi-level model, but I will include them for clariry when I write the hierarchical model. We have seen in 2(b)(ii) what the correlation structure will look like, depending on whether we are using `(Instrument + Voice + Harmony | Subject)` or `(Instrument + Voice + Harmony || Subject)`.  

#### Multilevel model:


\def\c#1{\mbox{\scriptsize\it #1}}
\begin{eqnarray*}
\mbox{Classical}_i & = & \alpha_{0j[i]} + 
\alpha_{\c{guitar},j[i]} \cdot 1_{\{\c{Instrument}_i=\c{guitar}\}} + 
\alpha_{\c{piano},j[i]} \cdot 1_{\{\c{Instrument}_i=\c{piano}\}} + 
\alpha_{\c{string},j[i]} \cdot 1_{\{\c{Instrument}_i=\c{string}\}} + \\ & &  
\alpha_{\c{contrary},j[i]} \cdot 1_{\{\c{Voice}_i=\c{contracy}\}} + 
\alpha_{\c{par3rd},j[i]} \cdot 1_{\{\c{Voice}_i=\c{par3rd}\}} + 
\alpha_{\c{par5th},j[i]} \cdot 1_{\{\c{Voice}_i=\c{par5th}\}} + \\ & & 
\alpha_{\c{I-IV-V},j[i]} \cdot 1_{\{\c{Harmony}_i=\c{I-IV-V}\}} + 
\alpha_{\c{I-V-IV},j[i]} \cdot 1_{\{\c{Harmony}_i=\c{I-V-IV}\}} + \\ & & 
\alpha_{\c{I-V-VI},j[i]} \cdot 1_{\{\c{Harmony}_i=\c{I-V-VI}\}} + 
\alpha_{\c{IV-I-V},j[i]} \cdot 1_{\{\c{Harmony}_i=\c{IV-I-V}\}} + \epsilon_i ~,~~ \epsilon_i \stackrel{iid}{\sim} N(0,\sigma^2)\\
\alpha_{0j} & = & \beta_0 + \eta_{0j} ~,~~ \eta_{0j} \stackrel{iid}{\sim} N(0,\tau_0^2) \\
\alpha_{\c{guitar},j} & = & \beta_{\c{guitar}} + \eta_{\c{guitar},j}~,~~ 
   \eta_{\c{guitar},j} \stackrel{iid}{\sim} N(0,\tau_{\c{guitar}}^2) \\
\alpha_{\c{piano},j} & = &  \beta_{\c{piano}} + \eta_{\c{piano},j}~,~~ 
   \eta_{\c{piano},j} \stackrel{iid}{\sim} N(0,\tau_{\c{piano}}^2) \\
\alpha_{\c{string},j} & = & \beta_{\c{string}} + \eta_{\c{string},j}~,~~ 
   \eta_{\c{string},j} \stackrel{iid}{\sim} N(0,\tau_{\c{string}}^2) \\
\alpha_{\c{contrary},j} & = & \beta_{\c{contrary}} + \eta_{\c{contrary},j}~,~~ 
   \eta_{\c{contrary},j} \stackrel{iid}{\sim} N(0,\tau_{\c{contrary}}^2) \\
\alpha_{\c{par3rd},j} & = & \beta_{\c{par3rd}} + \eta_{\c{par3rd},j}~,~~ 
   \eta_{\c{par3rd},j}  \stackrel{iid}{\sim} N(0,\tau_{\c{par3rd}}^2) \\
\alpha_{\c{par5th},j} & = & \beta_{\c{par5th}} + \eta_{\c{par5th},j}~,~~ 
   \eta_{\c{par5th},j}  \stackrel{iid}{\sim} N(0,\tau_{\c{par5th}}^2) \\
\alpha_{\c{I-IV-V},j} & = & \beta_{\c{I-IV-V}} + \eta_{\c{I-IV-V},j}~,~~ 
   \eta_{\c{I-IV-V},j} \stackrel{iid}{\sim} N(0,\tau_{\c{I-IV-V}}^2) \\
\alpha_{\c{I-V-IV},j} & = & \beta_{\c{I-V-IV}} + \eta_{\c{I-V-IV},j}~,~~ 
   \eta_{\c{I-V-IV},j} \stackrel{iid}{\sim} N(0,\tau_{\c{I-V-IV}}^2) \\
\alpha_{\c{I-V-VI},j} & = & \beta_{\c{I-V-VI}} + \eta_{\c{I-V-VI},j}~,~~ 
   \eta_{\c{I-V-VI},j} \stackrel{iid}{\sim} N(0,\tau_{\c{I-V-VI}}^2) \\
\alpha_{\c{IV-I-V},j} & = & \beta_{\c{IV-I-V}} + \eta_{\c{IV-I-V},j}~,~~ 
   \eta_{\c{IV-I-V},j} \stackrel{iid}{\sim} N(0,\tau_{\c{IV-I-V}}^2) 
\end{eqnarray*}


#### Hierarchical model:

It's actually somewhat difficult (or at least tedious) to write the Hierarchical model efficiently.  One way would be to define
\begin{eqnarray*}
m_i & = & \alpha_{0j[i]} + 
\alpha_{\c{guitar},j[i]} \cdot 1_{\{\c{Instrument}_i=\c{guitar}\}} + 
\alpha_{\c{piano},j[i]} \cdot 1_{\{\c{Instrument}_i=\c{piano}\}} + 
\alpha_{\c{string},j[i]} \cdot 1_{\{\c{Instrument}_i=\c{string}\}} + \\ & &  
\alpha_{\c{contrary},j[i]} \cdot 1_{\{\c{Voice}_i=\c{contracy}\}} + 
\alpha_{\c{par3rd},j[i]} \cdot 1_{\{\c{Voice}_i=\c{par3rd}\}} + 
\alpha_{\c{par5th},j[i]} \cdot 1_{\{\c{Voice}_i=\c{par5th}\}} + \\ & & 
\alpha_{\c{I-IV-V},j[i]} \cdot 1_{\{\c{Harmony}_i=\c{I-IV-V}\}} + 
\alpha_{\c{I-V-IV},j[i]} \cdot 1_{\{\c{Harmony}_i=\c{I-V-IV}\}} + \\ & & 
\alpha_{\c{I-V-VI},j[i]} \cdot 1_{\{\c{Harmony}_i=\c{I-V-VI}\}} + 
\alpha_{\c{IV-I-V},j[i]} \cdot 1_{\{\c{Harmony}_i=\c{IV-I-V}\}}
\end{eqnarray*}
and then write, ignoring the correlations again,
\begin{eqnarray*}
Classical_i & \sim & N(m_i,\sigma^2)\\
\alpha_{0j} & \sim & \beta_0 + \eta_{0j} ~,~~ \eta_{0j} \stackrel{iid}{\sim} N(0,\tau_0^2) \\
\alpha_{\c{guitar},j} & \sim & N(\beta_{\c{guitar}},\tau_{\c{guitar}}^2) \\
\alpha_{\c{piano},j} & \sim &  N(\beta_{\c{piano}},\tau_{\c{piano}}^2) \\
\alpha_{\c{string},j} & \sim & N(\beta_{\c{string}},\tau_{\c{string}}^2) \\
\alpha_{\c{contrary},j} & \sim & N(\beta_{\c{contrary}},\tau_{\c{contrary}}^2) \\
\alpha_{\c{par3rd},j} & \sim & N(\beta_{\c{par3rd}},\tau_{\c{par3rd}}^2) \\
\alpha_{\c{par5th},j} & \sim & N(\beta_{\c{par5th}},\tau_{\c{par5th}}^2) \\
\alpha_{\c{I-IV-V},j} & \sim & N(\beta_{\c{I-IV-V}},\tau_{\c{I-IV-V}}^2) \\
\alpha_{\c{I-V-IV},j} & \sim & N(\beta_{\c{I-V-IV}},\tau_{\c{I-V-IV}}^2) \\
\alpha_{\c{I-V-VI},j} & \sim & N(\beta_{\c{I-V-VI}},\tau_{\c{I-V-VI}}^2) \\
\alpha_{\c{IV-I-V},j} & \sim & N(\beta_{\c{IV-I-V}},\tau_{\c{IV-I-V}}^2) 
\end{eqnarray*}
A fully correct specification would express the random vector of $\alpha$'s as a multivariate normal, like this for the model with random effects `(Instrument + Voice + Harmony || Subject)`, for example:

\small
\[
\hspace*{-0.75in}\left[\begin{array}{c}
\alpha_{0j} \\
\alpha_{\c{guitar},j} \\
\alpha_{\c{piano},j}  \\
\alpha_{\c{string},j} \\
\alpha_{\c{contrary},j} \\
\alpha_{\c{par3rd},j} \\
\alpha_{\c{par5th},j} \\
\alpha_{\c{I-IV-V},j} \\
\alpha_{\c{I-V-IV},j} \\
\alpha_{\c{I-V-VI},j} \\
\alpha_{\c{IV-I-V},j} 
\end{array}\right] 
\stackrel{iid}{\sim}
N\left(
\left[\begin{array}{c}
\beta_{0j} \\
\beta_{\c{guitar}} \\
\beta_{\c{piano}} \\
\beta_{\c{string}} \\
\beta_{\c{contrary}} \\
\beta_{\c{par3rd}} \\
\beta_{\c{par5th}} \\
\beta_{\c{I-IV-V}} \\
\beta_{\c{I-V-IV}} \\
\beta_{\c{I-V-VI}} \\
\beta_{\c{IV-I-V}} 
\end{array}\right] ,
\left[\begin{array}{ccccccccccc}
\tau_{0j}  \\
0 & \tau_{\c{guitar}}^2  \\
0 & \rho_{gp} & \tau_{\c{piano}}^2  \\
0 & \rho_{gs} &\rho_{ps} & \tau_{\c{string}}^2 & \multicolumn{7}{c}{\mbox{\em(symmetric)}} \\
0 & 0 & 0 & 0 & \tau_{\c{contrary}}^2 \\
0 & 0 & 0 & 0 & \rho_{c3} & \tau_{\c{par3rd}}^2 \\
0 & 0 & 0 & 0 & \rho_{c5} & \rho_{35} & \tau_{\c{par5th}}^2 \\
0 & 0 & 0 & 0 & 0 & 0 & 0 & \tau_{\c{I-IV-V}}^2 \\
0 & 0 & 0 & 0 & 0 & 0 & 0 & \rho_{145,154} & \tau_{\c{I-V-IV}}^2 \\
0 & 0 & 0 & 0 & 0 & 0 & 0 & \rho_{145,156} & \rho_{154,156} & \tau_{\c{I-V-VI}}^2 \\
0 & 0 & 0 & 0 & 0 & 0 & 0 & \rho_{145,415} & \rho_{154,415} & \rho_{156,415} & \tau_{\c{IV-I-V}}^2 
\end{array}\right] 
\right)
\]
\normalsize

\noindent
and the model with random effects `(Instrument + Voice + Harmony | Subject)` would look the same, except the 0's in the variance-covariance matrix above would be replaced with additional correlations. 


The estimated $\beta$'s, $\tau$'s and $\rho$'s can be read off the display() output below.
\begin{itemize}
\item The $\hat\beta$'s and their SE's can be read from the ``coef.est'' and ``coef.se'' columns in the first table in the output
\item The $\hat\tau$'s and $\hat\rho$'s can be read from the second table, labelled ``Error terms:'', with the $\hat\tau$'s in the ``Std.Dev'' column, and the nonzero $\hat\rho$'s in the columns labelled "corr".
\end{itemize}

```{r}
display(lmer.ni.ivh.c0)
```


\rm



Because they are design variables in the experiment, the main effects for the three experimental
factors, Instrument, Harmony, and Voice, should be included in all models 
for the remainder of this homework, regardless of what you found about their
influence or lack of influence on ratings.  (Interactions among the design variables need not be included if you do not find them important, but the main effects for the design variables should \textbf{always} be in your models.)




# Problem 3.

\ul{\textit{Person covariates.}}  For this problem, begin 
with your best model from problem 2.

## Problem 3(a).

Determine which person covariates should be added to the 
model as fixed effects.    
Show a suitable summary of your work, 
and list the final set of variables that you would 
include in the model.  \textit{Hint: Some covariates that are actually
factor variables are coded as numeric.  Be careful to treat them as factors!}
\label{p:3a}

\it

So, again, I'm going to start with the model 
```{r eval=F}
lmer.ni.ivh.c0 <- lmer(Classical ~ Instrument + Voice + Harmony + 
                               (Instrument + Voice + Harmony || Subject), data=fulldata)
```
and I note that the "person covariates" are level-2 covariates (persons are the "groups" in this analysis) so they will only enter as fixed effects, not random effects in the `(...|Subject)` notation.

Here is my strategy:
\begin{itemize}
\item I will work with the smaller data set `nomissdata` since I need a data set that will be the same no matter which person covariates I put in the model, and in the larger data set `fulldata` many of the person covariates have missing values (so the actual data set would change, depending on what covariates are in the model and where the missing values for those covariates were, as I am comparing models).


\item I will start with `fitLMER.fnc` from `library(LMERConvenienceFunctions)`.  This does backwards selection of fixed effects, followed by forward selection of random effects, and then one more pass of backwards selection of fixed effects.  I will have it skip the forward selection of random effects (by not suggesting any new random effects for it to try).  I will also need to specify that `fitLMER.fnc` should use AIC or BIC for variable selection.

\item If `fitLMER.fnc` doesn't produce useful results, I will try some sort of forward selection of fixed effects by hand.  (But it turns out that `fitLMER.fnc` does just fine with my data set `nomissdata`, so I didn't try forward selection by hand...)

\end{itemize}

```{r}
library(LMERConvenienceFunctions)

nomissdata$Classical <- as.numeric(nomissdata$Classical)
nomissdata$Popular <- as.numeric(nomissdata$Popular)
```

Here's what happens if I treat all the predictors (except for the three design factors Instrument, Voice, and Harmony) as continuous.  This is not the right way to handle this problem (the categorical variables should be converted to factors!), but it is good "practice" before trying to convert predictors to factors..
```{r}

big.predictors <- names(nomissdata)[-c(1,2,24,25,26)]

## I am omitting the following variables from the set of predictors 
## that I want to start with, for backward selection:

names(nomissdata)[c(1,2,24,25,26)]

## Now build the model formula...

big.fla <- paste("Classical ~", paste(big.predictors,collapse=" + "),
                 "+ (Instrument + Voice + Harmony || Subject)")

lmer.big.cts <- lmer(big.fla, data=nomissdata, REML=F)

## We can already get a pretty good idea of what's going on by looking at the 
## display() function and the fixed effects coefficient table from summary()

display(lmer.big.cts)

summary(lmer.big.cts)$coef

## now let's see what fitLMER.fnc does with it...

lmer.AIC.cts <- fitLMER.fnc(lmer.big.cts, backfit.on="t", method="AIC", prune.ranefs=F,
                        log.file.name="lmerAICcts.log")

display(lmer.AIC.cts)

summary(lmer.AIC.cts)$coef

lmer.BIC.cts <- fitLMER.fnc(lmer.big.cts, backfit.on="t", method="BIC", prune.ranefs=F,
                        log.file.name="lmerBICcts.log")

display(lmer.BIC.cts)

summary(lmer.BIC.cts)$coef


```
Both AIC and BIC kept all three design factors (I would have added them back in if they hadn't!), and in fact they also kept the same person covariates in this case: "Selfdeclare", "OMSI", "X16.minus.17", "ConsNotes", "ClsListen" and "PianoPlay".  Also the correlations among the random effects look very similar to the correlations from our work from problem 2.

Now let's try this after converting the predictors that should be categorical to categorical variables:

```{r}
nomissdata.cat <- nomissdata

names(nomissdata.cat)

## Let's see which variables "look like" they should be categorical:

apply(nomissdata.cat[,8:24],2,table)

## It's a little hard to judge exactly which variables should be categorical for sure.
## I think the ones that definitely should be are: "PachListen", "ClsListen",
## "KnowRob", "KnowAxis", "CollegeMusic", "APTheory", "Composing",
## "PianoPlay" and "GuitarPlay". The variable "NoClass" (number of music classes 
## taken) should definitely be quantitative.   The others are kind of ambiguous to 
## me (either because Jimemez applied arithmetic to them, or they take on values
## other than non-negative integers), and I am going to leave them as numeric/quantitative.

cat.vars <- c("PachListen", "ClsListen", "KnowRob", "KnowAxis", "CollegeMusic",
              "APTheory", "Composing", "PianoPlay", "GuitarPlay")

nomissdata.cat[,cat.vars] <- apply(nomissdata.cat[,cat.vars],2,as.factor)

str(nomissdata.cat)
```
Now I am ready to try `fitLMER.fnc` on the new data set with categorical variables...
```{r}
lmer.big.cat <- lmer(big.fla, data=nomissdata.cat, REML=F)

## It's harder to see what's going on by just listing the fitted parameters, 
## because of all the levels of the categorical variables...

display(lmer.big.cat)

summary(lmer.big.cat)$coef

## so let's see what fitLMER.fnc does with it...

lmer.AIC.cat <- fitLMER.fnc(lmer.big.cat, backfit.on="t", method="AIC", prune.ranefs=F,
                        log.file.name="lmerAICcat.log")

display(lmer.AIC.cat)

summary(lmer.AIC.cat)$coef

lmer.BIC.cat <- fitLMER.fnc(lmer.big.cat, backfit.on="t", method="BIC", prune.ranefs=F,
                        log.file.name="lmerBICcat.log")

display(lmer.BIC.cat)

summary(lmer.BIC.cat)$coef

formula(lmer.AIC.cat)
formula(lmer.BIC.cat)

```
Once again, AIC and BIC picked the same models, but they are somewhat different from the models chosen when we treated all the predictors as continuous.  The person-level covariates chosen were:
\begin{itemize}
\item X16.minus.17
\item ConsNotes
\item PachListen
\item ClsListen
\item KnowRob
\item KnowAxis
\item X1990s2000s
\item X1990s2000s.minus.1960s1970s
\item NoClass
\item APTheory
\item Composing
\item GuitarPlay 
\end{itemize}
Based on the `fitLMER.fnc` work, these are the person level covariates that I would add to the model.


#### Aside #1: 

When we look at display(lmer.AIC.cat) or display(lmer.BIC.cat), we that there is really only one $\eta$ shared among the levels of Voice (because of the correlations of +1 and $-1$), and so it would be nice to implement this restriction in the model.  Unfortunately it is not straightforward to use lmer() to restrict $\eta$'s in this way... 

#### Aside #2:

We haven't looked at any residual plots.  Let's at least look at residual plots for `lmer.AIC.cat`, that we are reasonably happy with

Note: For some reason, the hlm_resid function in library(HLMdiag) wasn't producing the 2nd level residuals ($\eta$'s), and so I just extracted raw residuals directly:
\begin{itemize}
\item Conditional residuals are in `residuals(fitted.lmer.model)`
\item $\eta$'s are in `ranef(fitted.lmer.model)`
\end{itemize}

```{r fig.height=10,fig.width=8}
res.cond <- residuals(lmer.AIC.cat)
fit.cond <- fitted(lmer.AIC.cat)

eta <- ranef(lmer.AIC.cat)$Subject
fit.grp <- sapply(split(fit.cond,nomissdata.cat$Subject),mean)

par(mfcol=c(4,3))

## 1

plot(res.cond ~ fit.cond,xlab="Conditional Fitted Values Per Listener", ylab="Conditional Residuals",
     main="Level 1 Residuals")
abline(h=0)

qqnorm(res.cond,ylab="Conditional Residual QUantiles", main="Level 1 Residuals")
qqline(res.cond)

## 2

plot(eta$"(Intercept)" ~ fit.grp,xlab="Mean Fitted Values Per Listener",ylab=expression(eta),
     main=expression(paste("Intercept ",eta)))
abline(h=0)

qqnorm(eta$"(Intercept)",ylab=expression(paste(eta," QUantiles")), 
       main=expression(paste("Intercept ",eta)))
qqline(eta$"(Intercept)")

## 3

plot(eta$Instrumentguitar ~ fit.grp,xlab="Mean Fitted Values Per Listener",ylab=expression(eta),
     main=expression(paste("Instrument=guitar ",eta)))
abline(h=0)

qqnorm(eta$Instrumentguitar,ylab=expression(paste(eta," QUantiles")), 
       main=expression(paste("Instrument=guitar ",eta)))
qqline(eta$Instrumentguitar)

## 4 

plot(eta$Instrumentpiano ~ fit.grp,xlab="Mean Fitted Values Per Listener",ylab=expression(eta),
     main=expression(paste("Instrument=piano ",eta)))
abline(h=0)

qqnorm(eta$Instrumentpiano,ylab=expression(paste(eta," QUantiles")), 
       main=expression(paste("Instrument=piano ",eta)))
qqline(eta$Instrumentpiano)

## 5 

plot(eta$Instrumentstring ~ fit.grp,xlab="Mean Fitted Values Per Listener",ylab=expression(eta),
     main=expression(paste("Instrument=string ",eta)))
abline(h=0)

qqnorm(eta$Instrumentstring,ylab=expression(paste(eta," QUantiles")), 
       main=expression(paste("Instrument=string ",eta)))
qqline(eta$Instrumentstring)

## 6

plot(eta$Voicecontrary ~ fit.grp,xlab="Mean Fitted Values Per Listener",ylab=expression(eta),
     main=expression(paste("Voice=contrary ",eta)))
abline(h=0)

qqnorm(eta$Voicecontrary,ylab=expression(paste(eta," QUantiles")), 
       main=expression(paste("Voice=contrary ",eta)))
qqline(eta$Voicecontrary)

```


```{r fig.height=10,fig.width=8}

par(mfcol=c(4,3))


## 7

plot(eta$Voicepar3rd ~ fit.grp,xlab="Mean Fitted Values Per Listener",ylab=expression(eta),
     main=expression(paste("Voice=par3rd ",eta)))
abline(h=0)

qqnorm(eta$Voicepar3rd,ylab=expression(paste(eta," QUantiles")), 
       main=expression(paste("Voice=par3rd ",eta)))
qqline(eta$Voicepar3rd)

## 8

plot(eta$Voicepar5th~ fit.grp,xlab="Mean Fitted Values Per Listener",ylab=expression(eta),
     main=expression(paste("Voice=par5th ",eta)))
abline(h=0)

qqnorm(eta$Voicepar5th,ylab=expression(paste(eta," QUantiles")), 
       main=expression(paste("Voice=par5th ",eta)))
qqline(eta$Voicepar5th)

## 9

plot(eta$"HarmonyI-IV-V"  ~ fit.grp,xlab="Mean Fitted Values Per Listener",ylab=expression(eta),
     main=expression(paste("Harmony=I-IV-V ",eta)))
abline(h=0)

qqnorm(eta$"HarmonyI-IV-V",ylab=expression(paste(eta," QUantiles")), 
       main=expression(paste("Harmony=I-IV-V ",eta)))
qqline(eta$"HarmonyI-IV-V")

## 10

plot(eta$"HarmonyI-V-IV" ~ fit.grp,xlab="Mean Fitted Values Per Listener",ylab=expression(eta),
     main=expression(paste("Harmony=I-V-IV ",eta)))
abline(h=0)

qqnorm(eta$"HarmonyI-V-IV",ylab=expression(paste(eta," QUantiles")), 
       main=expression(paste("Harmony=I-V-IV ",eta)))
qqline(eta$"HarmonyI-V-IV")

## 11

plot(eta$"HarmonyI-V-VI" ~ fit.grp,xlab="Mean Fitted Values Per Listener",ylab=expression(eta),
     main=expression(paste("Harmony=I-V-VI ",eta)))
abline(h=0)

qqnorm(eta$"HarmonyI-V-VI",ylab=expression(paste(eta," QUantiles")), 
       main=expression(paste("Harmony=I-V-VI ",eta)))
qqline(eta$"HarmonyI-V-VI")

## 12

plot(eta$"HarmonyIV-I-V" ~ fit.grp,xlab="Mean Fitted Values Per Listener",ylab=expression(eta),
     main=expression(paste("Harmony=IV-I-V ",eta)))
abline(h=0)

qqnorm(eta$"HarmonyIV-I-V",ylab=expression(paste(eta," QUantiles")), 
       main=expression(paste("Harmony=IV-I-V ",eta)))
qqline(eta$"HarmonyIV-I-V")

```
All these residual plots adn qq plots look pretty good.  The only one that looks a bit strange is the first plot, i.e., the plot of the level-1 conditional residuals.  It is seriated (points fall along lines within the plot), but that is essentially because all of the fixed effect predictors are discrete. 

\rm


## Problem 3(b).

Once the fixed effects are settled, go back and check to see whether
there should be any change in the random effects; focus on random slopes for design variables and their interactions.  Provide suitable 
evidence to justify your results.
\label{p:3b}

\it

Let's go ahead and just do the same model selection exercise we did before, for the random effects...

```{r}
lmer.AIC.i <- update(lmer.AIC.cat, . ~ . - (0+Voice|Subject) - (0+Harmony|Subject)) 
lmer.AIC.v <- update(lmer.AIC.cat, . ~ . - (0+Instrument|Subject) - (0+Harmony|Subject))
lmer.AIC.h <- update(lmer.AIC.cat, . ~ . - (0+Instrument|Subject) - (0+Voice|Subject))
lmer.AIC.iv <- update(lmer.AIC.cat, . ~ . - (0+Harmony|Subject))
lmer.AIC.ih <- update(lmer.AIC.cat, . ~ . - (0+Voice|Subject))
lmer.AIC.vh <- update(lmer.AIC.cat, . ~ . - (0+Instrument|Subject))
lmer.AIC.ivh <- lmer.AIC.cat

print(x <- rbind(
  lmer.AIC.i=c(DIC=DIC.ml(lmer.AIC.i)),
  lmer.AIC.v=c(DIC=DIC.ml(lmer.AIC.v)),
  lmer.AIC.h=c(DIC=DIC.ml(lmer.AIC.h)),
  lmer.AIC.iv=c(DIC=DIC.ml(lmer.AIC.iv)),
  lmer.AIC.ih=c(DIC=DIC.ml(lmer.AIC.ih)),
  lmer.AIC.vh=c(DIC=DIC.ml(lmer.AIC.vh)),
  lmer.AIC.ivh=c(DIC=DIC.ml(lmer.AIC.ivh)))
  )

```
We just end up with `lmer.AIC.ivh`, which is the same model that we got from backwards elimination in part (a).  


\rm



## Problem 3(c).

Briefly interpret the effect of each variable 
kept in the final model, on Classical ratings.
\label{p:3c}

\it


From the summary of the fixed effect estimates (the $\hat\beta$'s),
```{r}
round(summary(lmer.AIC.ivh)$coef,2)

library(car)

vif(lmer.AIC.ivh)
```
There don't appear to be serious collinearities, so we can infer from the fixed effects coeffient table above that
\begin{itemize}
\item The harmony I-V-IV adds nearly a point on average to the Classical rating, vs. the "baseline" harmony I-IV-V.  The other harmonies do not really change the Classical rating much. \bf Listeners generally thing I-V-IV is a more classical sounding harmony structure.\it
\item The instrument "piano" adds over 1.5 points to the Classical rating, vs. the baseline instrument "guitar".  The instrument "string" adds over 3.5 points to the Classical rating, vs. "guitar".  \bf Listeners think piano, and especially strings, sound more classical than guitar.\it
\item The voice-leading patterns "parallel 3rds" and "parallel 5ths" reduce the Classical rating by roughly 0.4 and 0.3 points, respectively, relative to the baseline "contrary motion". \bf Listeners think contrary motion sounds mildly more classical than the other two voice-leading patterns.\it
\item We don't know what `X16.minus17` is, but it has a small negative effect on Classical ratings.
\item People who concentrate on the notes when listening have classical ratings which are on average about 0.3 lower ($\hat\beta_{ConsNotes}=-0.3$), than people who do not concentrate on the notes.
\item The pattern with PachListen (familiarity with Pachelbel's Canon) and ClsListen (how much do you listen to classical music) is not straightforward.  For example, listeners who scored 4 out of 5 on these two measures added nearly 2, and nearly 7 points, respectively, to their Classical ratings, on average, relative to the baseline category of 0 ( = `` not at all'') on each variable.  But listeners who scored 5 out of 5 on these measures actually reduced their Classical ratings somewhat, on average, relative to category 0.  However, these may be due in part to a small-sample bias issue, since only one listener scored 4 on either PachListen or ClsListen. \bf(See calculations below.)\it
\item The patterns with KnowRob (have you heard Rob Paravonian’s Pachelbel Rant) and KnowAxis (have you heard Axis of Evil’s Comedy bit on the 4 Pachelbel chords in popular music?) is more straightforward.  Relative to baseline category 0 (= ``not at all''),  listeners who express modest familiarity with the Rob Paravonian bit have lower Classical ratings by about 0.36 points, whereas those with modest familiarity with the Axis of Awesome bit, reduce their classical ratings by more than six points. The giant jump ($\hat\beta_{KNowAxis=1}-6.30$) may be mainly a small-sample bias issue, since only one listener scored in category 1 on this variable.   \bf(See calculations below.)\it
\item Each additional point of familiarity with music from the 1990's conveys a small reduction in Classical ratings on average, but each additional point in the contrast between familiarity with music from the 1990's vs familiarity withmusic from the 1960's conveys a similar-sized increase in Classical ratings; it seems likely that these two effects basically cancel each other out.
\item Each additional Music class that a listener has taken increases their Classical scores by about 0.74; on the other hand, having taken AP Theory conveys a 2.64 increase in classical ratings. 
\item The pattern of estimated coefficients for Composing is again strange, with a score of 5 on the Composing variable associated with a 2.76 point drop in Classical ratings, breaking the pattern established by the other levels of Composing.  Once again, this may be mostly a small-sample bias issue since only one listener is in category 5 on the Composing variable.  \bf(See calculations below.)\it
\item The pattern of coefficients for GuitarPlay is also hard to understand, but the coefficients for categories 2 and 4 may be in part a consequence of small sample bias.  \bf(See calculations below.)\it
\end{itemize}
Here are the calculations that show that categories of some of the predictors with anomalous coefficients also involve only one or two examinees, and so the anomalies may be simply small-sample bias issues:
```{r}
table(nomissdata.cat$PachListen)
with(nomissdata.cat,table(Subject[PachListen=="4"]))

table(nomissdata.cat$ClsListen)
with(nomissdata.cat,table(Subject[ClsListen=="4"]))

table(nomissdata.cat$KnowAxis)
with(nomissdata.cat,table(Subject[KnowAxis=="1"]))

table(nomissdata.cat$Composing)
with(nomissdata.cat,table(Subject[Composing=="5"]))

table(nomissdata.cat$GuitarPlay)
with(nomissdata.cat,table(Subject[GuitarPlay=="2"]))
with(nomissdata.cat,table(Subject[GuitarPlay=="4"]))

```

\rm




# Problem 4.

\ul{\textit{Musicians vs. Non-musicians.}}  One of the secondary 
hypotheses of the researchers is that people who self-identify as musicians 
may be influenced by things that do not influence non-musicians.  Dichotomize
 "Selfdeclare" ("are you a musician?") so 
that about half the participants are categorized as self-declared musicians, 
and half not. Examine and report on any interactions between the dichotomized 
musician variable and other predictors in the model.  Check to see if the results are sensitive to where you dichotomize.  Provide suitable evidence 
for, and comment on, your results.
\label{p:4}

\it

Interstingly, `Selfdeclare` was not one of the variables in the final model selected by `fitLMER.fnc` above.  Here is the result of dichotomizing is and trying to include it in the model.

Since there are 1541 observations in the `nomissdata.cat` data frame and, as shown in the table below, about half are accounted for in the first two categories of `Selfdeclare`, we'll start by dichotomizing between 2 and 3.
```{r}
cumsum(table(nomissdata.cat$Selfdeclare))

table(Musician <- with(nomissdata.cat,ifelse(Selfdeclare<=2,0,1)))

big.fla.with.Musician <- paste("Classical ~ (", paste(big.predictors,collapse=" + "),
                 "- Selfdeclare)*Musician + (Instrument + Voice + Harmony || Subject)")

lmer.big.with.Musician <- lmer(big.fla.with.Musician,data=nomissdata.cat)

display(lmer.big.with.Musician)

round(summary(lmer.big.with.Musician)$coef,2)

large.fixed.effects <- function(M) {
  x <- as.data.frame(summary(M)$coef)
  return(round(x[abs(x$"t value")>=2.00,],2))
}

large.fixed.effects(lmer.big.with.Musician)


```
We see that 
\begin{itemize}
\item Some coefficients are very poorly estimated (e.g. estimated coefficients in the 100's or 1000's even though the response only ranges from 0 to 10).  This is most often caused by severe collinearity and/or small sample size for estimating those particular coefficients.
\item Many possible interactions with Musician are dropped because of collinearities in the model; and 
\item Musician does seem to have a significant interaction the design variables Instrument, Voice and Harmony.  The interactions really only seem to involve one or two levels of each design variable.  
\end{itemize}

Next, we see what happens if we try to dichotomize Musician at different places:
```{r}
table(nomissdata.cat$Selfdeclare)

table(Musician <- with(nomissdata.cat,ifelse(Selfdeclare<=1,0,1)))
lmer.Musician.1 <- lmer(big.fla.with.Musician,data=nomissdata.cat)
large.fixed.effects(lmer.Musician.1)

## This is our original dichotomization:
table(Musician <- with(nomissdata.cat,ifelse(Selfdeclare<=2,0,1)))
lmer.Musician.2 <- lmer(big.fla.with.Musician,data=nomissdata.cat)
large.fixed.effects(lmer.Musician.2)

table(Musician <- with(nomissdata.cat,ifelse(Selfdeclare<=3,0,1)))
lmer.Musician.3 <- lmer(big.fla.with.Musician,data=nomissdata.cat)
large.fixed.effects(lmer.Musician.3)

table(Musician <- with(nomissdata.cat,ifelse(Selfdeclare<=4,0,1)))
lmer.Musician.4 <- lmer(big.fla.with.Musician,data=nomissdata.cat)
large.fixed.effects(lmer.Musician.4)

table(Musician <- with(nomissdata.cat,ifelse(Selfdeclare<=5,0,1)))
lmer.Musician.5 <- lmer(big.fla.with.Musician,data=nomissdata.cat)
large.fixed.effects(lmer.Musician.5)

```
It's clear that the interactions with Musician vary greatly, depending on where we dichotomize Selfdeclare to get the Musician variable.  Our original dichotomization produces the most parsimonious result (only interactions with the three design variables), but exactly where to dichotomize would depend on a conversation with Dr. Jimenez.

\rm



# Problem 5.

\ul{\textit{Classical vs. Popular.}}  It would be totally appropriate to repeat problems 2, 3 and 4 using "Popular" ratings as the response, rather than classical, and then compare the results for "Classical" ratings, vs. "Popular" ratings.  \textit{In the interest of time, \textbf{you do not have to do this} for this hw/project.  It would make good "future work" for the paper; and/or it is something you could pursue on your own, if you are curious about it.}
\label{p:5}

\it

This is already a very long set of solutions, so I won't repeat any of the work using `Popular` as a response variable! This can just be "future work"!

\rm



# Problem 6.

\ul{\textit{Brief writeup.}} Write an approximately one or two page professional-quality summary\footnote{This should form raw materials for the \textbf{Discussion} section of your final IDMRAD paper, so please review what should be in the \textbf{Discussion} section of an IMDRAD or IMRAD paper.  Also, try to keep this close to one or two pages.} 
of your findings, suitable for Dr. Jimenez.  
Be sure to address: 
\begin{itemize}
\item The influence of the three main experimental 
factors (Instrument, Harmony \& Voice);
\item A brief discussion of variance components---is this a standard repeated
measures model, or did we need to include other variance components?  
\item A discussion of other indivdual covariates in the model.
\end{itemize}
You may
refer to your earlier work (e.g. "As I showed in my answer to part 2(b), 
blah-blah-blah\ldots").  Don't be sloppy about the statistical findings, but 
try to highlight things that will be of substantive interest to Dr. Jimenez.  Make
your summary very readable and clear.
\label{p:6}

\it

The raw data provided by Dr. Jimenez consisted of 2520 observations of 28 variable each, fully balanced across 70 listeners/subjects and 36 experimental conditions.  Observations are groups of 36 audio music stimuli, grouped by listener.  So, as a multi-level data set, the level-1 observations are ratings of how Classical the music in each stimulus sounds, and the level-2 groups are the set of stimuli that each listener rated . 

From the raw data, we created two data sets:
\begin{itemize}
\item `fulldata`, a data set consisting of all of the original variables except for `X1stInstr` and `X1stInstr`, because they had so much missing data (60\% and 87\%, respectively) that they were nearly useless.  We also removed 27 rows from `fulldata` that had missing Classical and Popular response values, and recode two responses of 19 to 10, assuming that the 19's were intended to be 10's but were keying errors (typo's) when entering the data.   This data set has 2493 observations of 26 variables, and is nearly balanced across the 70 subjects and 36 experimental conditions.
\item `nomissdata`, the largest subset of `fulldata` with no missing values.  This data set has 1541 complete observations of 36 variables.  It is reasonably well balanced across experimental conditions and subjects, but represents only 43 of the original 70 subjects.
\end{itemize}
The distributions of the variables in `fulldata` and `nomissdata` are similar.  Although two or three quantitative variables exhibit noticable skewing, we decided not to do any transformations, to make the results as immediately interpretable for Dr. Jimenez as possible.

With the `fulldata` data set, we explored the influence of the three experimental factors, both fixed and random effects, on Classical ratings.

With the `nomissdata` data set, we did variable selection on the person covariate fixed effects and re-considered random effects in light of a final person covariates model.

EDA and preparation of the data sets is shown in the work for Problem 1 above, pp. \pageref{p:1}ff. 

### The influence of the three main experimental factors (Instrument, Harmony \& Voice).

We considered the design factors Instrument, Harmony and Voice using standard linear models with no random effects (Problem 2(a), pp. \pageref{p:2a}ff.), and using random-intercept multi-level models (Problem 2(b)(ii), pp. \pageref{p:2b.ii}ff.).  Overall, the results were
\begin{itemize}
\item In both the standard linear models and in the multi-level models,
\begin{itemize}
\item The model selection criterion AIC (which approximately minimizes prediction error) strongly prefers a model with main effects for Instrument, Voice and Harmony, and an interaction between Voice and Harmony. 
\item The model selection criterion BIC (which approximately minmimizes mis-specification of the model) strongly prefers the model with just the main effects for Instrument, Voice and Harmony.
\end{itemize}
\item The model selection criterion DIC (which is designed to correctly count degrees of freedom when comparing models with different random effects) very strongly preferred models with a random intercept, rather than standard linear models without a random intercept (see Problem 2(b)(iii), pp. \pageref{p:2b.iii} above).
\end{itemize}

The final model from this exploration has main effects for Instrument, Voice and Harmony, possibly an interaction between Voice and Harmony, and a random intercept.  This is a standard repeated-measures model, where the value of the random intercept for each listener represents that listener's overall tendency, over and above the features of the stimulus, to rate any stimulus as sounding `Classical'.   Table~\ref{t:initial-models} summarizes the two models, with and without the interaction of Voice and Harmony:

```{r results='asis',echo=F}
texreg(list(lmer.no.interactions,lmer.vh),stars=numeric(0),
       caption="Initial Random-Intercept Models",
       single.row=T,
       booktabs = T,
       dcolumn=T,
       longtable=T,
       label="t:initial-models",
       custom.model.names=c("Additive Fixed Effects","Additive Fixed Effects + Voice*Harmony"),
       use.packages=F)
```



### Variance components

We then explored other random effects.  Instead of a single random intercept that varies from listener to listener, we considered an expanded set of random effects allowing the coefficients on each level of each design variable to vary from listener to listener as random effects.  This would allow us to accommodate variations in personal tendencies to rate someting as classical, where those personal tendencies vary from from one instrument to another, from one voice to another, or from one harmony to another.

Exploration and comparison of models with various combinations of random slopes on levels of the design variables (problem 2(c)(i), pp. \pageref{p:2c.i}ff.) led to the conclusion that random slopes on all the design variables fit the data better.  We continued to work with a main effects only model, and a model that adds the Voice by Harmony interaction.  Because of a high degree of correlation between the random effects, these models become somehat difficult to interpret. As an alternative we also considered models in which the random effects for different design variables are assumed to be uncorrelated.  Models with correlations forced to be zero, summarized in Table~\ref{t:sel-on-design-REs}, do not fit as well as models that leave all the correlations free, but they are more interpretable and simpler to work with, so we will continue to assume zero correlations between random effects for different desigbn factors.

```{r results='asis',echo=F}
texreg(list(lmer.ni.ivh.c0,lmer.vh.ivh.c0),stars=numeric(0),
       caption="Models with expanded variance components", 
       single.row=T,
       booktabs = T,
       dcolumn=T,
       longtable=T,
       label="t:sel-on-design-REs",
       custom.model.names=c("Additive Fixed Effects","Additive Fixed Effects + Voice*Harmony"),
       fontsize = "small",
       use.packages=F)
```
These models generalize the simple repeated measures model by allowing listeners' individual tendencies to rate music as Classical to depend in part on the features of the simulus music. 
In the remainder of the work we used the simpler model on the left as a basis for further exploration. 

### Person-level covariates.

Backwards variable elimination using either AIC or BIC (Problem 3(a), pp. \pageref{p:3a} ff.), followed by a check of the random effects (Problem 3(b), pp. \pageref{p:3b}), leads to the same "final" model, displayed in Table~\ref{t:sel-on-pers-covs}.

```{r results='asis',echo=F}
texreg(lmer.AIC.ivh,stars=numeric(0),
       caption="Model after selection of person covariates",
       single.row=T,
       booktabs = T,
       dcolumn=T,
       longtable=T,
       label="t:sel-on-pers-covs",
       custom.model.names=c("Additive Fixed Effects"),
       use.packages=F)
```
This model fits the data well, relative to others considered, and has good residual diagnostic plots. Once again, correlations between random effects on different design variables are set to be zero here.   Below, we interpret the fixed effects in the model (see also Problem 3(c), pp. \pageref{p:3c}ff.):
\begin{itemize}
\item The effects with clearest interpretation:
\begin{itemize}
\item The harmony I-V-IV adds nearly a point on average to the Classical rating, vs. the "baseline" harmony I-IV-V.  The other harmonies do not really change the Classical rating much. \bf Listeners generally thing I-V-IV is a more classical sounding harmony structure.\it
\item The instrument "piano" adds over 1.5 points to the Classical rating, vs. the baseline instrument "guitar".  The instrument "string" adds over 3.5 points to the Classical rating, vs. "guitar".  \bf Listeners think piano, and especially strings, sound more classical than guitar.\it
\item The voice-leading patterns "parallel 3rds" and "parallel 5ths" reduce the Classical rating by roughly 0.4 and 0.3 points, respectively, relative to the baseline "contrary motion". \bf Listeners think contrary motion sounds mildly more classical than the other two voice-leading patterns.\it
\item We don't know what `X16.minus17` is, but it has a small negative effect on Classical ratings.
\item People who concentrate on the notes when listening have classical ratings which are on average about 0.3 lower ($\hat\beta_{ConsNotes}=-0.3$), than people who do not concentrate on the notes.
\end{itemize}
\item The effects with more ambiguous interpretation: 
\begin{itemize}
\item The pattern with PachListen (familiarity with Pachelbel's Canon) and ClsListen (how much do you listen to classical music) is not straightforward.  For example, listeners who scored 4 out of 5 on these two measures added nearly 2, and nearly 7 points, respectively, to their Classical ratings, on average, relative to the baseline category of 0 ( = `` not at all'') on each variable.  But listeners who scored 5 out of 5 on these measures actually reduced their Classical ratings somewhat, on average, relative to category 0.  However, these may be due in part to a small-sample bias issue, since only one listener scored 4 on either PachListen or ClsListen. 
\item The patterns with KnowRob (have you heard Rob Paravonian’s Pachelbel Rant) and KnowAxis (have you heard Axis of Evil’s Comedy bit on the 4 Pachelbel chords in popular music?) is more straightforward.  Relative to baseline category 0 (= ``not at all''),  listeners who express modest familiarity with the Rob Paravonian bit have lower Classical ratings by about 0.36 points, whereas those with modest familiarity with the Axis of Awesome bit, reduce their classical ratings by more than six points. The giant jump ($\hat\beta_{KNowAxis=1}-6.30$) may be mainly a small-sample bias issue, since only one listener scored in category 1 on this variable.    
\item Each additional point of familiarity with music from the 1990's conveys a small reduction in Classical ratings on average, but each additional point in the contrast between familiarity with music from the 1990's vs familiarity withmusic from the 1960's conveys a similar-sized increase in Classical ratings; it seems likely that these two effects basically cancel each other out.
\item Each additional Music class that a listener has taken increases their Classical scores by about 0.74; on the other hand, having taken AP Theory conveys a 2.64 increase in classical ratings. 
\item The pattern of estimated coefficients for Composing is again strange, with a score of 5 on the Composing variable associated with a 2.76 point drop in Classical ratings, breaking the pattern established by the other levels of Composing.  Once again, this may be mostly a small-sample bias issue since only one listener is in category 5 on the Composing variable.  
\item The pattern of coefficients for GuitarPlay is also hard to understand, but the coefficients for categories 2 and 4 may be in part a consequence of small sample bias.  
\end{itemize}
\end{itemize}

We also looked at whether a dichotomized version of the `Selfdeclare` variable had interesting interactions with other variables in the data set (Problem 4, pp. \pageref{p:4}ff.).  The dichotomized variable we used was
\[
     Musician_i = \left\{\begin{array}{rl}
                  1, & \mbox{\it if $Selfdeclare\leq c$}\\
                  0, & \mbox{\it otherwise}
                  \end{array}\right.
\]
We found that if $c$ is chosen so that about half the sample in `nomissdata` consider themselves to be Musicians, then we get interesting interactions with the three design variables (only).  That is, Musicians tend to rate certain instruments, voice leading or harmony structures as classsical differently from how non-musicians rate them.  

However, the presence or absence of interactions of the dichotomized 
"Musician" variable with design variables or with other person covariates depends strongly on the value of $c$ in the above equation.  Therefore, a final determination of the question of interactions of a dichotomized version of the `Selfdeclare` with design factors and other terms in the model should really wait for a careful conversation with Dr. Jimenez. 


### Future work.

Here, briefly,  are two items to consider for future work:
\begin{itemize}
\item I did not have time to go back and incoprorate the interaction between Voice and Harmony in the later work above.  It's interesting to note that in the earlier work AIC always wants to include it, and BIC does not, and so ultimately including this interaction probably depends on whether the model is being built to produce good predictions, or to produce and accurately specified model that can help explain listeners' reactions to specific musical features.

\item I also did not try interactions among the person covariates in a model like that of Table~\ref{t:sel-on-pers-covs}, other than the ``Musician'' exploration above.  This would also be interesting to pursue further.

\item All listeners in this study also rated all musical stimuli on how Popular they sounded.  It would be very interesting to repeat all of the above analyses using Popular as the response variable, and compare those analyses with the analyses above.
\end{itemize}

\rm

