---
title: "350 Week 7b  Multilevel modeling"
author: "William Revelle"
date: "05/09/2022"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
options(width=100)
```




```{r}
#install.packages("psych",repos="http://personality-project.org/r", type="source")
library(psych) #make it active
library(psychTools)
sessionInfo()   #psych should be 2.2.5 or higher
```
# Multilevel modeling : An Introduction

One of the more exciting changes in the way that we conduct psychological research is the introduction of intensive longitudinal modeling using *Experience Sampling Methododology* (*ESM*) which allows us to collect *state* data over time.  Alhough originally done with various diary procedures, it is now typically done with smart phones and apps.  At the same time, the procedures for evaluating longitudinal data have seen a computional revolution as computational approaches of *multilevel modeling* have been introduced.  

Here we go through the *R* procedures for doing this analysis.

See the accompanying slides at 
<a href="http:\\personality-project.org/courses/350/350.week7.dynamics.pdf"> 
http:\\personality-project.org/courses/350/350.dynamics </a> and the accompanying article at 
<a href = "http://personality-project.org/revelle/publications/rw.tutorial.19.pdf"> http://personality-project.org/revelle/publications/rw.tutorial.19.pdf </a>.


## A toy problem

It is convenient, when testing new procedures, to *simulate* data.  We do this using the `sim.multi` function.

We set the *random* seed to a specific number \citep[e.g., 42, ][]{adams:hitch} in order to produce replicable results.  If the seed is not specified, a different pseudo random sequence will be generated for each simulation. 

The simulation is set to create data for four subjects with two measures (V1 and V2 and V3 and V4) of each of two factors. 

Items are assumed to have loadings of .6 on one factor, 0 on the other.  The factor intercorrelations at the overall level are randomly varying around 0, but differ for each subject, with factor intercorrelations of -.7, 0, 0, and .7 respectively. 

Although normally we would not round the results of the simulation we do so to create the data in Tables~\ref{tab:fat} and ~\ref{tab:wide}. In addition, the normal output of the \pfun{sim.multi} is already in wide form.  We convert it to Fat form and then back again as demonstrations of the \fun{reshape} and \pfun{dfOrder} functions.  

```{r}
library(psych)   #activate the psych package
#create the data
set.seed(42)
x <- sim.multi(n.obs=4,nvar=4,nfact=2,days=6,ntrials=6,plot=TRUE, 
          phi.i=c(-.7,0,0,.7),loading=.6)
raw <- round(x[3:8])
raw[1:4] <- raw[1:4] + 6
#make a 'Fat' version
XFat <- reshape(raw,idvar="id",timevar="time",times=1:4,direction="wide")
#show it
XFat

#now make it wide
XWide <- reshape(XFat,idvar="id",varying=2:25,direction="long")
Xwide <- dfOrder(XWide,"id")
Xwide   #show the data in typical wide format
#note that there are 4 subjects, with 6 observations per subject on each of 4 variables

#add in the trait information
traits <- data.frame(id = 1:4,extraversion =c(5,10,15,20),
                       neuroticism =c(10,5, 15,10))
traits   #this would be some data collected separately
Xwide.traits <- merge(Xwide,traits, by ="id")  #this merges two data.frames using the id
Xwide.traits   #show the merged data
```

## First, find the descriptive statistics

Basic descriptive statistics may be found by using the appropriate functions in the \Rpkg{psych} package.  We find overall descriptives (\pfun{describe}), descriptives by group (individual), and then  a variety of multilevel statistics including the pooled within person correlations, the between person correlations, and the individual within individual correlations using the \pfun{statsBy} function.  We plot the data using the \pfun{mlPlot} function.

```{r}
help(describe)    #get help on the input commands for the describe function.
#first the overall descriptives
describe(Xwide.traits)
#do it again, but by subjects
describeBy(Xwide.traits, group="id")

#now find the within and between individual correlations
sb <- statsBy(Xwide.traits, group="id",cors=TRUE)
#show them
lowerMat(sb$rwg)   #top part of Table 6
lowerMat(sb$rbg)   #lower part of Table 6
round(sb$within,2)   #Table 7
#
#plot the data
mlPlot(Xwide.traits,grp="id",items =3:6) 
```

# Multilevel Reliability

Reliability is the question of how much variance in a measure is not error. When examining multi-level data, the question is a bit harder, because we can think about several different reliabilities.  Within subjects, across subjects, across time. ...

Traditional reliability measures decompose between person  test variance ($\sigma^2_x$)  into reliable and un-reliable or error variance ($\sigma^2_e$). 

\begin{equation}
\rho_{xx} = \frac{1 - \sigma^2_e}{\sigma^2_x}
\end{equation} 

The problem is then how to find $\sigma^2_e$.  Solutions include test retest correlations and estimates based upon the internal structure of the test \citep{guttman:45,rz:09}.  Generalizability theory \citep{gleser:65} was developed to solve the problem of multiple sources of variance for each score.  This is the situation for multi-level data.  For when subjects have scores on multiple items and multiple time points, we are interested in several different indices of reliability.   The technique developed by \cite{gleser:65} was to estimate the variance components using standard ANOVA procedures to find the appropriate Mean Squares for subjects, time, items, etc.  and then convert these to variance components based upon the expected values of the MS (Table~\ref{tab:anova} top part).  Taking advantage of the power of \R{} to integrate the output of one function into another, we combine `aov` `lme` and
`lmer` into 
 one function (`multilevel.reliability` or `mlr`) 
 which can take wide data, transform it into the long format needed for `aov`etc., do the analyses, and then find the reliabilities based upon these components and the formulae given by \cite{shrout:12a}. Thus the command to find multilevel reliability for a set of variables is just one line of code rather than the complex expressions necessary in SPSS or SAS \citep{shrout:12a}.


An excellent tutorial on this is found in a chapter by Pat Shrout and Sean Lane.  We use their toy example for the analysis.

## Get the Shrout and Lane data

A convenent way to store and transmit small data sets is using the `dput` function which *puts* the data into a readable and compact form.  We can read it directly back into R.

```{r}
shrout <- structure(list(Person = c(1L, 2L, 3L, 4L, 5L, 1L, 2L, 3L, 4L, 
5L, 1L, 2L, 3L, 4L, 5L, 1L, 2L, 3L, 4L, 5L), Time = c(1L, 1L, 
1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 
4L, 4L), Item1 = c(2L, 3L, 6L, 3L, 7L, 3L, 5L, 6L, 3L, 8L, 4L, 
4L, 7L, 5L, 6L, 1L, 5L, 8L, 8L, 6L), Item2 = c(3L, 4L, 6L, 4L, 
8L, 3L, 7L, 7L, 5L, 8L, 2L, 6L, 8L, 6L, 7L, 3L, 9L, 9L, 7L, 8L
), Item3 = c(6L, 4L, 5L, 3L, 7L, 4L, 7L, 8L, 9L, 9L, 5L, 7L, 
9L, 7L, 8L, 4L, 7L, 9L, 9L, 6L)), .Names = c("Person", "Time", 
"Item1", "Item2", "Item3"), class = "data.frame", row.names = c(NA, 
-20L))

headTail(shrout)  #first and last 4 rows of the data frame in normal wide format

#or just show the entire data set
shrout
```


We use the `multilevel.reliability` function (from `psych`) to analyze the data.
Various reliability (or generalizability) coefficients may be found from these variance components.  For instance, the reliability of individual differences over k fixed time points and m multiple items is 
\begin{equation}
R_{kF} = \frac{\sigma^2_{id} + (\sigma^2_{id x items}/m)}{\sigma^2_{id} + (\sigma^2_{id x items}/m) + \sigma^2_{error}/(k m)}
\label{eq:Rkf}
\end{equation}
Equation~\ref{eq:Rkf} is just one (\#6)  of the six generalizability coefficients discussed by \cite{shrout:12a}.   


```{r}
mg <- multilevel.reliability(shrout,grp="Person",Time="Time",items=
        c("Item1","Item2","Item3"),plot=TRUE)
mg
```

As is true of most R functions, there is a lot more there than meets the eye:

```{r}
names(mg) #list the various objects included in the mg object
        
headTail(mg$long)  #show the first and last rows of the long object
```

One purpose of studying multi-level data is to explore individual differences in changes over time.  People can differ in their affect as a function of the situation \citep{wilt:rev:situations}, their scores can increase or decrease over time, and they can show diurnal variation in their moods or their performance.  To find the slope over time, we simply apply a within subject regresssion, and to examine the phase and fit of diurnal variation, we use circular statistics \citep{Jammalamadaka:2006,pewsey} using the `cosinor` function. Later, 
we discuss how to generate and analyze  simulated data with a variety of data structures, particularly growth and decay over time, and diurnal variation.


# A real data set (taken from  Aaron Fisher at UCB)

In a study of 10 participants diagnosed with clinically generalized anxiety disorder, \cite{fisher:15} collected 28 items for at least 60 days per participant.  In an impressive demonstration of how different people are, \cite{fisher:15} examined the dynamic factor structure of each person using procedures discussed by \cite{molenaar:85,molenaar:09}.  The purpose of the study was to encourage the use of personalized care for clinical psychology.  Most importantly, for our purposes, the data set he analyzed is available at Aaron Fisher's website <a href = "https://www.dynamicpsychlab.com"> https://www.dynamicpsychlab.com </a>
for easy download and subsequent analysis. (He seems to have moved the data from there.  But a copy of his data is available at <a href = "http://personality-project.org/courses/350/Fisher_2015_Data"> http://personality-project.org/courses/350/Fisher_2015_data"</a> )

As discussed in the Appendix to our tutorial, the 10 data files available may be merged into one file and we can examine both the reliability of scales made up of subsets of items, as well as the correlational pattern of these scales.  It is important to note that the original paper goes far beyond what we report here, and indeed, the analyses that follow are independent of the main thrust of Fisher's paper.

Of the 28 items, we focus on just eight, four measuring general positive affect (happy, content, relaxed, and positive), and four measuring tension and negative affect (angry, afraid, sad, and lonely).  We see that the rating suggest a clearly reliable  separation between individuals for both positive and negative affects (Table~\ref{tab:fisher:reliability}).  Scoring each subject for these two moods at all time points may done using `scoreItems`

First, we create a little function to our work for us.  This is partially tailored for the Fisher data set, but the concept of concatinating files remains the same for any such arrangement.

We can examine the file structure outside of R by using our browser to see the structure.

### First create a small function to get the data

We use a few R commands such as `getwd`, `setwd` to adjust our working directory from the default location. The function *loops* over all the subjects.
It dynamically calculates how many subjects there are by looking at the length of the `names` parameter.  With a bit of effort, we could even make this automatic.
```{r}

"combine.data" <- function(dir=NULL,names,filename=NULL) {
  new <- NULL
  n <- length(names)
  old.dir <- getwd()   #save the current working directory
for (subject in 1:n) {  #repeat n times, once for each subject
if(is.null(filename)) {setwd(dir)}  else {dir <- filename}     #set the working directory to where the files are
  #this is specific to this particular data structure
  x <- read.file(f=paste0(dir,"/P",names[subject],"/pre",names[subject],".csv")) 
  nx <- nrow(x)
   #add id and time to this data frame
  temp <- data.frame(id=names[subject],time=1:nx,x) 
  new <- rbind(new,temp)  #combine with prior data.frames to make a longer object
   }  #end of the subject loop
 setwd(old.dir)   #set the working directory back to the original
 return(new)}   #end the function by returning the data
```
### We use this general function for the specific case of Fisher's data

We do a little bit of work to get the files from Fisher's server to ours (not shown) and then run our function. The function uses the `read.file` function from `psychTools`.

When we run the function, it echoes out the name of each file read. 

Note that Fisher used a missing value of 999.  We use the `scrub` function t replace all 999 with NA.

```{r}
#now use this function to read in data from a set of files (downloaded from Fisher's web page) and
#  combine into one data.frame
names <- c("002","007","009","010","011","013","022","023","030","065")  #hard coded from his file names
filename="http://personality-project.org/courses/350/Fisher_2015_Data"     #specify where the data are
new <- combine.data(dir=NULL, names=names,filename=filename)    #use the function we just created for this problem
describe(new) #note how missing values were coded as 999 
fisher <- scrub(new,2:29,isvalue=999)   #change 999 to NA 
dim(fisher)  #show the number of rows and columns
describe(fisher) #to see what is there and to check that we recoded everything correctly.
```
The data set contains 28 different mood words.  For the purpose of the demonstration, we want to find the reliability of four positive affect terms (happy, content, relaxed, and positive) as well as the four negaitve affect terms (Table~\ref{tab:fisher:reliability}).   We then find scores for all subjects at all time periods by using the
`scoreItems`
function.  We combine these scores with the id and time information from the original file, and then plot it with the 
`mlPlot` function.   Finally, to examine the within subject correlations we use the
`statsBy` function with the cors=TRUE option (Table~\ref{tab:fisher:alpha}).

## Preliminary analyses

Before we do the multilevel analysis, we can treat the data as a simple data set.  Lets find the correlations of the positive and negative affect terms.
This shows the standard pattern of independence between positive and negative affect.  But this is an overall effect, what happens within subjects?
```{r}
positive <- cs(happy,content,relaxed, positive)
negative <- cs(angry,afraid, sad, lonely)
pana <- c(positive,negative)
R <- lowerCor(fisher[pana])
corPlot(R,numbers=TRUE)
```

Now lets find the simple reliabilities of these two sets of items.  We use the *scoreItems* function.

```{r}
pana.scores <- scoreItems(keys=list(positive=positive,negative=negative), fisher)
pana.scores
```
Looking at the correlations between the scales, we see that there is a slight negative correlation between positive and negative affect for the  pooled data. Lets show that using `pairs.panels` on the scores.

```{r pairs}
pairs.panels(pana.scores$scores)
```


# Now, do  the multilevel analyses

The `mlr` function (aka `multilevel.reliabilty`) function takes advantage of two very powerful packages for doing mixed effects `lmer` and `nlme'.  These need to be installed before you can run `mlr`.  They are loaded automatically when `mlr` starts up.

We also make use of the `lattice` graphic functions which are more powerful than simple core-R graphics but not as powerful as `ggplot`.  

```{r}
#find the multilevel reliabilities
pos.fisher.rel <- mlr(fisher,"id","time",items = c("happy","content", 
                                "relaxed", "positive"), aov=FALSE,lmer=TRUE)
neg.fisher.rel <- mlr(fisher,"id" ,"time",items = c("angry","afraid",
                         "sad","lonely" ),aov= FALSE,lmer=TRUE)
#organize the alpha's for each subject

alpha.df <-  data.frame(positive = pos.fisher.rel$alpha[,c(1,3)], 
                     negative = neg.fisher.rel$alpha[,c(1,3)])


select <-c("happy","content", "relaxed", "positive" ,
               "angry","afraid","sad","lonely" )
affect.keys <- list(positive =  c("happy","content", "relaxed", "positive"),
                    negative =  c("angry","afraid","sad","lonely" ) )
affect.scores <- scoreItems(keys= affect.keys, items = fisher[select], min=0, 
                             max=100)
affect.df <-cbind(fisher[1:2], affect.scores$score)

affect.df <- char2numeric(affect.df,flag=FALSE) #something is breaking statsBy
#by making everything numeric, we seem to solve this problem
#but we need to set flag to FALSE

mlPlot(affect.df, grp = "id",time="time", items= 3:4, typ="p", pch=c(20,21))

describe(affect.df) #debbugging purposes

stats.affect <- statsBy(affect.df,group="id",cors=TRUE)
stats.affect$within
# combine these with the alphas 
 ar <- autoR(affect.df[2:4],group= affect.df["id"])
 alpha.df <- cbind(alpha.df,stats.affect$within,ar$autoR[,2:3])
 alpha.df  #show the results
```




# Simulations as a way of creating test data sets

A very powerful tool in learning how to analyze data and to test how well various models work is to create simulated data set.  The `sim.multi` function allows for creating arbitrarily large data sets with a variety of complex data structures. The basic generator allows one to define a factor model with 1 ... f factors with 1.. j  items per factor and for 1 ... i subjects and 1 ... k time points. It is possible to define factor loadings globally, or individually for each subject.  Factors (and their corresponding items) can increase, remain stable, or decrease over time, and can show diurnal variation with different phase relationships.  Diurnal variation is simulated as a sinusoidal curve varying over 24 hours with a peak (phase angle) at different times of day.  

The default for `sim.multi` is for four subjects over 6 days on two variables:

```{r}
sim.multi()  #just use the default values
```

A more complicate 
example of 16 such subjects is seen in Figure~\ref{fig:diurnal} and various summary statistics are given in Table~\ref{tab:diurnal}. 
16 subjects were simulated with two measures taken 8 times a day for six days.  The scores for some subjects decreased, while for others they increased over the 6 days.  People differ in the strength (amplitude) and phase of their diurnal rhythms. 

We create 16 subjects, with two factors and three items per factor. We simulate data collected 6 times/day over 16 days.  We set the random seed to a fixed number for a reproducible example. After generating the data, we apply the `cosinor` function to estimate the diurnal nature of the signal, `statsBy` to find the within individual correlations, and `autoR` to find the auto correlations of the data.  We combine the various outputs into one data.frame to display in Table~\ref{tab:diurnal}.  


```{r}
set.seed(17)
 x <- sim.multi(n.obs=16,nvar=2,nfact=2,ntrials=48,days=6,
 sin.i =c(1,0,0,1),cos.i=c(1,.5,0,-.5), sigma.i=c(0,0,0,0), 
 sigma=1,phi.i=c(0,0,.5,.5),beta.i=c(-1.,0,1,0,-.5,0,.5,0))
co <-  cosinor(x$time,x[3:6],code="id")
sb <- statsBy(x,group="id",cors=TRUE)
aR <-autoR(x,group="id")
sim.df <- data.frame(co[,c(1,2,4,5)],sb$within[,c(8,9,10)],aR$autoR[,3:4])
sim.df  #show it

#to find the pooled data, we do the same analyses, 
#  but without specifying the group
cos <-  cosinor(x$time,x[3:6])  #pooled results
aR <-autoR(x)     #pooled results
rs <- cor(x[3:5])
rs <- rs[lower.tri(rs)]   #just use the relevant ones
pooled <- c(cos[1:2,1:2],rs,aR$autoR[3:4])
#combine the pooled data with the individual data
sim.df <- rbind(sim.df,pooled)   #rbind adds rows of data frames together
                                #see also cbind to add columns together
round(sim.df,2)  #print it more neatly
df2latex(sim.df)  #if you have LaTex, this makes a table for you.
```

# Random Coefficients on the Fisher data set

The `R` packages `nlme` \citep{nlme:16} and `lme4` \citep{lme4:15} handle a variety of multilevel modeling procedures and can be used to conduct random coefficient modeling (RCM), which is the formal term for models that vary at more than one level.

RCM is done in `nlme` with the `lme` function 
 and in `lme4` with the `lmer` function. These
   functions allow for the specification of fixed effects and random effects within and across multiple levels.  Therefore, main effects of traits and states can be modeled using these functions, as can random effects of states across traits (or any higher level variables).  
 
 To do this analysis we first found the mean positive and negative affect for each subject, and then centered these data around the individual's overall mean `negative.cent`). We see these effects (`lme`) and the random coefficients for each individual (extracted using the `coef` function) for the two variables (positive and negative affect) derived from the Fisher
data set in Table~\ref{tab:ml:effects}.  We see that  negative emotional *states* lead to lower positive emotions, but that the effect of *trait* negative emotion does not affect state positive affect.

We prepare the Fisher data for random coefficient modeling by computing aggregate means (level 2 data) for each participant, merging the level 2 data with the level 1 state data, group-mean centering the state predictor variable around each participant's mean, and merging the centered data.  Then we conduct random coefficient models are conducted in the `nlme`package with the `lme` fuction and in the `lme4`package with the `lmer` function. 
(We might need to install these packages first. )

Random coefficients for each participant are extracted with the `coef`function (Table~\ref{tab:ml:effects}).

### First, do some data *munging* or organization

Functions that work on nested data typically require using the `long` format, instead of the more traditional `wide` format.  Thus, we need to do some work to convert the data into the form we want.  There are several ways of doing this.

The `aggregate` function will apply an arbitrary function to each subset of the data defined by a list, the `statsBy` function will report statistics for each value of a grouping variable.


```{r}

#compute aggregate means for each participant
#affect.mean <- aggregate(fisher,list(affect.df$id),mean, na.rm = TRUE)  #one way

affect.mean <- statsBy(affect.df,group="id")$mean 

# another way
affect.mean.df <- data.frame(group= rownames(affect.mean),affect.mean)

#rename columns to prepare for merge
colnames(affect.mean.df) <- c("id","id.1","time.mean","positive.mean",
                        "negative.mean") 
#merge participant means with individual responses 
affect.mean.df <- merge(affect.df,affect.mean.df,by="id")
#group-mean center positive and negative affect
affect.centered <- affect.mean.df[,c(3,4)] - affect.mean.df[,c(7,8)]  
#rename columns to prepare for merge
colnames(affect.centered) <- c("positive.cent","negative.cent") 
#add centered variables to data frame
affect.mean.centered <- cbind(affect.mean.df,affect.centered)  

```

### Use the nlme package

(perhaps need to install it first)

```{r}
#using the nlme package
library(nlme)
#this model predicts positive affect from time and negative affect (centered), 
#and it allows the slopes relating positive affect to time and negative affect
# to  vary across participants

pa.na.time.nlme  <- lme(positive ~ time + negative.cent + negative.mean 
					+ negative.cent:negative.mean, 
					random= ~time + negative.cent|id, 
					data=affect.mean.centered, 
					na.action = na.omit)
					
summary(pa.na.time.nlme) #shows fixed and random effects
#extract the coefficients for each participant
coef.pa.time.nlme <- coef(pa.na.time.nlme)  
round(coef.pa.time.nlme,2)
describe(coef.pa.time.nlme)  #describe the coefficients

```

### Now use the lme4 package 

(might need to install)

```{r}
#using the lme4 package
library(lme4)
pa.na.time.lme4 <- lmer(positive  ~ time + negative.cent +  negative.mean + 
	 	negative.cent:negative.mean + (time|id) + (negative.cent|id), 
		data = affect.mean.centered, 
		na.action = na.omit)

summary(pa.na.time.lme4)  #the summary function gives the important results
coef.pa.na.time.lme4 <- coef(pa.na.time.lme4)
coef.pa.na.time.lme4
```

## Yet another way to get within subject patterns of results

The `statsBy` function from `psych` will give some useful statistics on the within subject correlational structure of the data.

```{r statsBy}
sb.affect <- statsBy(affect.df,group="id",cors=TRUE)
sb.affect #how much do the subjects differ?
names(sb.affect)  #what are the possible objects to explore?
round(sb.affect$within,2)  # what are the correlations within subjects
```
Compare these correlations with the graphics we drew earlier.

```{r lattice2}
mlPlot(affect.df, grp = "id",time="time", items= 3:4, typ="p", pch=c(20,21))
```


# Conclusions

Modern data collection techniques allow for intensive measurement within subjects.  Analyzing this type of data requires analyzing data at the within subject as well as between subject level.  Although sometimes conclusions will be the same at both levels, it is frequently the case that examining within subject data will show much more complex patterns of results than when they are simply aggregated.  This tutorial is  a simple introduction to the kind of data analytic strategies that are possible.