Chapter 3.1 - Bivariate Intraindividual Covariability Metrics (Part 1)
Overview
This tutorial extends the calculation of intraindividual variability metrics into examination of intraindividual covariation. Specifically, this tutorial demonstrates calculation of intraindividual correlation (and covariance) – as derived metrics that may be useful for articulating a variety of dynamic characteristics from experience sampling data. We begin the move towards the standard multilevel modeling analysis approach by first considering individual-level regressions. At the end of the session we begin setting up the mutlilevel model.
Outline
In this session we cover …
A. Preliminary description of two variables
B. Calculating intraindividual (within-person) correlations (iCorr, iCov)
C. Calculating intraindividual (within-person) regressions (iReg)
D. Setting up the Multilevel Model (MLM)
Prelim - Loading libraries used in this script.
library(psych)
library(ggplot2)
library(data.table)
library(entropy)
library(plyr)
library(nlme)Prelim - Reading in Repeated Measures Data
Note that we are working from a long file. For your own data, there may be some steps to get to this point.
#Setting the working directory
setwd("~/Desktop/Fall_2017") #Person 1 Computer
#setwd("~/Desktop/Fall_2017") #Person 2 Computer
#set filepath for data file
filepath <- "https://quantdev.ssri.psu.edu/sites/qdev/files/AMIBbrief_raw_daily1.csv"
#read in the .csv file using the url() function
daily <- read.csv(file=url(filepath),header=TRUE)
#Little bit of clean-up
var.names.daily <- tolower(colnames(daily))
colnames(daily)<-var.names.dailyEverything we do today uses a long file.
A: Looking at distributions and within-person time-series
Lets look at some descriptions of today’s variables.
#getting a list of the variable names
names(daily)## [1] "id" "day" "date" "slphrs" "weath" "lteq" "pss"
## [8] "se" "swls" "evalday" "posaff" "negaff" "temp" "hum"
## [15] "wind" "bar" "prec"#sample descriptives
describe(daily$posaff)## vars n mean sd median trimmed mad min max range skew kurtosis
## X1 1 1441 4.12 1.1 4.2 4.15 1.19 1 7 6 -0.25 -0.33
## se
## X1 0.03describe(daily$negaff)## vars n mean sd median trimmed mad min max range skew kurtosis
## X1 1 1441 2.45 1.04 2.2 2.34 1.04 1 6.9 5.9 0.96 0.77
## se
## X1 0.03#histograms
ggplot(data=daily, aes(x=posaff)) +
geom_histogram(fill="white", color="black") +
labs(x = "Positive Affect")## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.## Warning: Removed 17 rows containing non-finite values (stat_bin).
ggplot(data=daily, aes(x=negaff)) +
geom_histogram(fill="white", color="black") +
labs(x = "Negative Affect")## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.## Warning: Removed 17 rows containing non-finite values (stat_bin).
And lets look at the bivariate longitudinal time-series plots …
1. For Person A
#ggplot version .. see also http://ggplot.yhathq.com/docs/index.html
ggplot(data = subset(daily, id ==101), aes(x=day), legend=FALSE) +
#geom_rect(mapping=aes(xmin=day, xmax=day+1, ymin=0, ymax=100, fill=posaff), alpha=0.8) +
geom_point(aes(x=day,y = posaff), shape=16, size=3,colour="blue") +
geom_line(aes(x=day,y = posaff), lty=1, size=1,colour="blue") +
geom_point(aes(x=day,y = negaff), shape=17, size=3,colour="red") +
geom_line(aes(x=day,y = negaff), lty=1, size=1,colour="red") +
xlab("Day") +
ylab("Positive (blue) and Negative (red) Affect") + ylim(1,7) +
scale_x_continuous(breaks=seq(0,7,by=1)) 
- For Person B
ggplot(data = subset(daily, id ==102), aes(x=day), legend=FALSE) +
#geom_rect(mapping=aes(xmin=day, xmax=day+1, ymin=0, ymax=100, fill=posaff), alpha=0.8) +
geom_point(aes(x=day,y = posaff), shape=16, size=3,colour="blue") +
geom_line(aes(x=day,y = posaff), lty=1, size=1,colour="blue") +
geom_point(aes(x=day,y = negaff), shape=17, size=3,colour="red") +
geom_line(aes(x=day,y = negaff), lty=1, size=1,colour="red") +
xlab("Day") +
ylab("Positive (blue) and Negative (red) Affect") + ylim(1,7) +
scale_x_continuous(breaks=seq(0,7,by=1)) 
For both individuals we see fluctuations in both positive and negative affect. Within-person, those changes may or may not align well.
B: Calculating intraindividual (within-person) correlations (iCorr, iCov)
As with the intraindvidual variability metrics, calculation of intraindividual correlation (and covariation) as a form of feature extraction. The features (i.e., new variables) are extracted from each individuals’ ensemble of scores.
Lets look at the within-person relation of posaff and negaff.
#Scatter Plots by id
i <- 101
ggplot(data=daily[which(daily$id ==i),], aes(x=negaff, y=posaff)) +
geom_point() +
geom_smooth(method=lm, se=FALSE, fullrange=TRUE) +
stat_ellipse() +
xlab("Negative Afect") + xlim(1,7) +
ylab("Positive Afect") + ylim(1,7) +
ggtitle(paste("ID =",i))
i <- 102
ggplot(data=daily[which(daily$id ==i),], aes(x=negaff, y=posaff)) +
geom_point() +
geom_smooth(method=lm, se=FALSE, fullrange=TRUE) +
stat_ellipse() +
xlab("Negative Afect") + xlim(1,7) +
ylab("Positive Afect") + ylim(1,7) +
ggtitle(paste("ID =",i)) ## Warning: Removed 3 rows containing missing values (geom_smooth).
The within-person associations between postive and negative affect are depicted by the 95% ellipse and the slope of the line.
Intraindividual Correlation (Covariance)
One way to quantify the relation between x and y is the intraindividual correlation or covariance.
We can do this in the same way we did the other “IIV” summaries, using the data.table package.
We simply add the cor() and cov() functions to our list.
#Computing Individual parameters
#converting to a data.table object
daily.dt <- data.table(daily)
#collecting summary info by id
#note: posaff is a continuous variable; evalday is a binary variable
indiv.statbasic <- daily.dt[,list(imean.posaff=mean(posaff, na.rm=TRUE), #continuous variable
isd.posaff=sd(posaff, na.rm=TRUE),
imean.negaff=mean(negaff, na.rm=TRUE), #continuous variable
isd.negaff=sd(negaff, na.rm=TRUE),
icount.posaff=sum(!is.na(posaff)), #data coverage
icount.negaff=sum(!is.na(negaff)),
prop.evalday=mean(evalday, na.rm=TRUE), #binary variable
icount.evalday=sum(!is.na(evalday)),
ientropy.se=entropy(table(se, useNA="no")), #categorical variable
icorr.pana=cor(x=negaff,y=posaff, #correlation
use="pairwise.complete.obs", method="pearson"),
icov.pana=cov(x=negaff,y=posaff, #correlation
use="pairwise.complete.obs", method="pearson")
),
by=id]
#An alternative
#library(plyr)
#count.statbasic <- ddply(daily, "id", summarize, icount=sum(!is.na(posaff)))And lets look at the sample-level descriptives of the individual-level summaries.
describe(indiv.statbasic)## vars n mean sd median trimmed mad min max
## id 1 190 318.29 130.44 321.50 318.99 151.23 101.00 532.00
## imean.posaff 2 190 4.10 0.73 4.10 4.10 0.63 2.23 6.04
## isd.posaff 3 187 0.81 0.35 0.78 0.79 0.33 0.18 2.17
## imean.negaff 4 190 2.48 0.73 2.41 2.43 0.72 1.11 5.09
## isd.negaff 5 187 0.74 0.36 0.68 0.70 0.30 0.12 2.83
## icount.posaff 6 190 7.58 1.22 8.00 7.89 0.00 1.00 8.00
## icount.negaff 7 190 7.58 1.22 8.00 7.89 0.00 1.00 8.00
## prop.evalday 8 190 0.69 0.19 0.71 0.69 0.21 0.00 1.00
## icount.evalday 9 190 7.48 1.33 8.00 7.83 0.00 1.00 8.00
## ientropy.se 10 190 0.75 0.35 0.69 0.76 0.42 0.00 1.56
## icorr.pana 11 187 -0.41 0.43 -0.51 -0.45 0.47 -0.99 1.00
## icov.pana 12 187 -0.38 0.61 -0.22 -0.28 0.32 -4.86 0.54
## range skew kurtosis se
## id 431.00 -0.04 -1.09 9.46
## imean.posaff 3.81 0.04 0.08 0.05
## isd.posaff 1.98 0.84 1.08 0.03
## imean.negaff 3.98 0.68 0.45 0.05
## isd.negaff 2.70 1.51 5.23 0.03
## icount.posaff 7.00 -3.99 16.55 0.09
## icount.negaff 7.00 -3.99 16.55 0.09
## prop.evalday 1.00 -0.42 0.33 0.01
## icount.evalday 7.00 -3.45 12.32 0.10
## ientropy.se 1.56 -0.16 -0.11 0.03
## icorr.pana 1.99 0.85 0.21 0.03
## icov.pana 5.40 -3.40 17.79 0.04There are clearly individual differences, even if the variables are calculated on “small ensembles”. But, a correlation on 8 bivariate observations is highly questionable.
Lets look at the distribution of the correlations.
#histogram
ggplot(data=indiv.statbasic, aes(x=icorr.pana)) +
geom_histogram(fill="white", color="black") +
labs(x = "iCorr PA~NA")## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.## Warning: Removed 3 rows containing non-finite values (stat_bin).
Note that, by definition correlations are bounded, -1.0 to +1.0. Thus, if we want to use the iCorr as an outcome variable, we MUST do a transformation: Fisher’s r-to-z transformation.
\[z_{r} = \frac{1}{2} [log_{e}(1+r) - log_{e}(1-r)]\]
\[z_{r} = \frac{1}{2} ln\left(\frac{1+r}{1-r}\right)\] We can either calculate ourselves or trust the fisherz() function in the psych package.
indiv.statbasic$icorr.pana.z <- fisherz(indiv.statbasic$icorr.pana)
describe(indiv.statbasic$icorr.pana.z)## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 187 Inf NaN -0.56 -0.57 0.68 -2.58 Inf Inf NaN NaN NaNNote that there is a problem. We have an Inf score for id == 128. That would need to be fixed.
#histogram
ggplot(data=indiv.statbasic, aes(x=icorr.pana.z)) +
geom_histogram(fill="white", color="black") +
labs(x = "iCorr PA~NA: z-transform")## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.## Warning: Removed 5 rows containing non-finite values (stat_bin).
See how the distribution is better.
C: Calculating intraindividual (within-person) regressions (iReg)
An alternative measure of association is to run individual regressions and use the \(\beta_{i}\) coefficients as an individual differences variable.
We use the plyr library to break up the data by id and fit a linear regression to each piece of data. Then collect the coefficients for each individual.
#using the plyr library
# Break up data by "id", then fit the specified model to each piece and return a list
models <- dlply(daily, "id", function(df)
lm(posaff ~ negaff,
data = df,
na.action=na.exclude))
# Apply coef() to each model and return a data frame
icoeffs <- ldply(models, coef)
#rename the variables
names(icoeffs) <- c("id","iintercept.pana","islope.pana")
#Alternative
#collecting regression output by id
#indiv.reg <- daily.dt[,c(reg.1 = as.list(coef(lm(posaff ~ negaff))),
# reg.1.sigma = as.list(summary(lm(posaff ~ negaff))$sigma)),
# by=id]And lets look how the individual regression varaibles are correlated (across-persons) with the icorr.pana and icov.pana variables.
#merging together
indiv.statbasic <- merge(indiv.statbasic,icoeffs,by="id")
#correlation
cor(x=indiv.statbasic$icorr.pana, y=indiv.statbasic$islope.pana, use="pairwise.complete.obs", method="pearson")## [1] 0.8050574#sclatter plot
ggplot(data=indiv.statbasic, aes(x=icorr.pana, y=islope.pana)) +
geom_point() +
geom_smooth(method=lm, se=FALSE, fullrange=TRUE, color="red") +
stat_ellipse()## Warning: Removed 3 rows containing non-finite values (stat_smooth).## Warning: Removed 3 rows containing non-finite values (stat_ellipse).## Warning: Removed 3 rows containing missing values (geom_point).
We see pretty good concordance, with some outlier cases probably related to number of observations.
D. Setting up the Multilevel Model (MLM)
Both of the above approaches are prone to inflation of Type-1 error rate, becasue the measures of IntraCov (e.g., icorr, icov, ireg) are imperfect estimates. They have a standard error that is ignored in the Stage 2 analysis of individual differences.
The MultiLevel Model (MLM) was developed to accommodate the interdependence in nested data.
It is particularly useful for examining individual differences in bivariate associations - thus providing for more robust statistical inference about within-person associations (intraindividual covariation) and the between-person differences therein.
The nlme package contains tools for fitting mixed effects models, also known as multilevel (MLM) or hierarchical linear models (HLM).
The lme() function fits the MLMs
The fixed argument takes the fixed model
The random argument takes the random model
The data argument specifies the data sources
The na.action argument specifies what to do with missing data
We set up the basic model. The requires choosing which variable will the “outcome” and which variable will be the “predictor”. We choose posaff as the outcome variable because it is more normally distributed.
- First we fit the Unconditional Means model to see how much between-person and how much within-person variance there is in our outcome variable.
#Use the constant '1' to designate that only intercepts are being modeled
um.fit <- lme(fixed= posaff ~ 1,
random= ~ 1|id,
data=daily,
na.action=na.exclude)
summary(um.fit)## Linear mixed-effects model fit by REML
## Data: daily
## AIC BIC logLik
## 4043.118 4058.935 -2018.559
##
## Random effects:
## Formula: ~1 | id
## (Intercept) Residual
## StdDev: 0.6565496 0.8806322
##
## Fixed effects: posaff ~ 1
## Value Std.Error DF t-value p-value
## (Intercept) 4.108614 0.05324996 1251 77.15714 0
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -3.38303960 -0.56263785 0.06669917 0.61455439 3.65052987
##
## Number of Observations: 1441
## Number of Groups: 190Lets extract the random effects with the VarCorr() function
VarCorr(um.fit)## id = pdLogChol(1)
## Variance StdDev
## (Intercept) 0.4310574 0.6565496
## Residual 0.7755131 0.8806322We can then compute the intra-class correlation (ICC) as the ratio of the random intercept variance (between-person) to the total variance (between + within).
First lets store the random effect variances, which will be the first column of the VarCorr object (see above).
RandomEffects <- as.numeric(VarCorr(um.fit)[,1])
RandomEffects## [1] 0.4310574 0.7755131Next lets compute the ICC. It is the ratio of the random intercept variance (between-person var) over the total variance (between + within var)
ICC_between <- RandomEffects[1]/(RandomEffects[1]+RandomEffects[2])
ICC_between## [1] 0.3572584So, between-person variance = 35.73%, and within-person variance = 100 - 35.73 = 64.27% This means there is lots of within-person variance to model!
- Lets do the same breakdown for our predictor variable, using the Unconditional Means model to see how much between-person and how much within-person variance there is.
#Use the constant '1' to designate that only intercepts are being modeled
um.fit <- lme(fixed= negaff ~ 1, random= ~ 1|id,
data=daily,
na.action=na.exclude)
summary(um.fit)## Linear mixed-effects model fit by REML
## Data: daily
## AIC BIC logLik
## 3839.46 3855.277 -1916.73
##
## Random effects:
## Formula: ~1 | id
## (Intercept) Residual
## StdDev: 0.6534748 0.8140798
##
## Fixed effects: negaff ~ 1
## Value Std.Error DF t-value p-value
## (Intercept) 2.463675 0.0522867 1251 47.11858 0
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -3.8739086 -0.6122847 -0.1608485 0.4657840 3.9394042
##
## Number of Observations: 1441
## Number of Groups: 190#Extract the random effects
VarCorr(um.fit)## id = pdLogChol(1)
## Variance StdDev
## (Intercept) 0.4270294 0.6534748
## Residual 0.6627260 0.8140798#Compute the ICC
RandomEffects <- as.numeric(VarCorr(um.fit)[,1])
ICC_between <- RandomEffects[1]/(RandomEffects[1]+RandomEffects[2])
ICC_between## [1] 0.391858OK - that looks good as well. Generally, experience sampling data has quite a bit of within-person variance because researchers tend to measure variables that vary over time. They do not waste valuable design dollars on variables that do not need to be measured repeatedly.
- We need to split our predictor variable into “trait” and “state” components.
We calculate the imeans.
daily.imeans <- ddply(daily, "id", summarize, imean.negaff=mean(negaff, na.rm=TRUE),
imean.posaff=mean(posaff, na.rm=TRUE))
describe(daily.imeans)## vars n mean sd median trimmed mad min max
## id 1 190 318.29 130.44 321.50 318.99 151.23 101.00 532.00
## imean.negaff 2 190 2.48 0.73 2.41 2.43 0.72 1.11 5.09
## imean.posaff 3 190 4.10 0.73 4.10 4.10 0.63 2.23 6.04
## range skew kurtosis se
## id 431.00 -0.04 -1.09 9.46
## imean.negaff 3.98 0.68 0.45 0.05
## imean.posaff 3.81 0.04 0.08 0.05Note that we really only need to do all of this for the predictor variable. I am also doing it for the outcome because of OCD, and because of some later checking and plotting.
We calculate a sample-centered versions of these person-level variables for later use.
daily.imeans$imean.negaff.c <- scale(daily.imeans$imean.negaff,center=TRUE,scale=FALSE)
daily.imeans$imean.posaff.c <- scale(daily.imeans$imean.posaff,center=TRUE,scale=FALSE)
describe(daily.imeans)## vars n mean sd median trimmed mad min max
## id 1 190 318.29 130.44 321.50 318.99 151.23 101.00 532.00
## imean.negaff 2 190 2.48 0.73 2.41 2.43 0.72 1.11 5.09
## imean.posaff 3 190 4.10 0.73 4.10 4.10 0.63 2.23 6.04
## imean.negaff.c 4 190 0.00 0.73 -0.07 -0.05 0.72 -1.37 2.61
## imean.posaff.c 5 190 0.00 0.73 0.00 0.00 0.63 -1.88 1.94
## range skew kurtosis se
## id 431.00 -0.04 -1.09 9.46
## imean.negaff 3.98 0.68 0.45 0.05
## imean.posaff 3.81 0.04 0.08 0.05
## imean.negaff.c 3.98 0.68 0.45 0.05
## imean.posaff.c 3.81 0.04 0.08 0.05We then merge these “trait” scores back into the long data file and calculate the “state” scores.
daily <- merge(daily,daily.imeans,by="id")
daily$negaff.state <- daily$negaff - daily$imean.negaff
daily$posaff.state <- daily$posaff - daily$imean.posaffLets plot the state and trait scores to see if it looks right.
#Trait Scores
ggplot(data = daily[which(daily$id <=105),], aes(x = day, y = imean.negaff, group = id, color=factor(id))) +
geom_line(data=daily[which(daily$id <= 105 & daily$imean.negaff !="NA"),]) +
xlab("Day") +
ylab("Trait Negative Affect") + ylim(1,7) +
scale_x_continuous(breaks=seq(0,7,by=1)) 
#State Scores (note change of the ylim)
ggplot(data = daily[which(daily$id <=105),], aes(x = day, y = negaff.state, group = id, color=factor(id))) +
geom_line(data=daily[which(daily$id <= 105 & daily$negaff.state !="NA"),]) +
xlab("Day") +
ylab("State Negative Affect") + ylim(-3,3) +
scale_x_continuous(breaks=seq(0,7,by=1)) 
OK - all looking good.
- Lets look at the between-person differences in within-person associations.
Since we split our predictor varaible into two pieces, we probably need to include both components into the regression. The state piece gives us the within-person association, and the trait piece gives us the between-person association.
model1.fit <- lme(fixed= posaff ~ 1 + negaff.state + imean.negaff,
random= ~ 1 + negaff.state|id,
data=daily,
na.action=na.exclude)
summary(model1.fit)## Linear mixed-effects model fit by REML
## Data: daily
## AIC BIC logLik
## 3570.623 3607.52 -1778.312
##
## Random effects:
## Formula: ~1 + negaff.state | id
## Structure: General positive-definite, Log-Cholesky parametrization
## StdDev Corr
## (Intercept) 0.6181958 (Intr)
## negaff.state 0.2782322 -0.165
## Residual 0.7103830
##
## Fixed effects: posaff ~ 1 + negaff.state + imean.negaff
## Value Std.Error DF t-value p-value
## (Intercept) 5.117124 0.17292438 1250 29.591684 0
## negaff.state -0.558281 0.03509220 1250 -15.908970 0
## imean.negaff -0.409654 0.06730408 188 -6.086621 0
## Correlation:
## (Intr) ngff.s
## negaff.state -0.003
## imean.negaff -0.959 -0.022
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -4.3004541 -0.5791329 0.0218156 0.5863607 3.2437095
##
## Number of Observations: 1441
## Number of Groups: 190Lets do some interpretation.
The Intercept parameter is the expected level of posaff when everything else is zero, which in this case for a person with trait negaff = 0 (i.e., the average person) on a day where state negaff = 0 (the average day). For us, that is 5.12.
The imean.negaff parameter is the expected difference in posaff for a person with 1 unit higher trait negaff. This is the between-person association \(\gamma_{01} = -0.41\), p < .05.
Individuals with higher trait negative affect have lower trait positive affect.
The negaff.state parameter is the expected difference in posaff on a day with 1 unit higher state negaff for the prototypical person. This is the average within-person association, \(\gamma_{10} = -0.56\), p < .05.
On days with higher state negative affect, the protopyical individual will have lower state positive affect.
- Lets plot the different kinds of associations.
For the between-person association, we have 190 persons with 1 estimated trait score (= imean) per person. That is 190 dots.
ggplot(data=daily.imeans, aes(x=imean.negaff, y=imean.posaff)) +
geom_point() +
geom_smooth(method=lm, se=FALSE, fullrange=TRUE, lty=1, size=2, color="red") +
xlab("Trait Negative Affect") + xlim(1,7) +
ylab("Trait Positive Affect") + ylim(1,7) +
ggtitle("Between-person Association")
Note that I have “fudged” this plot. It is not derived from the model, but from the background data files (and a sample-level regression outside the MLM).
For the within-person associations, we have 190 persons. That is 1 estimated regression line for each person = 190 lines (each derived from up to 8 bivariate observations).
ggplot(data=daily, aes(x=negaff.state, y=posaff.state, group=id, color="gray")) +
#geom_point() +
geom_smooth(method=lm, se=FALSE, fullrange=FALSE, lty=1, size=.5, color="gray40") +
geom_smooth(aes(group=1), method=lm, se=FALSE, fullrange=FALSE, lty=1, size=2, color="hotpink") +
xlab("State Negative Affect") + xlim(-4,4) +
ylab("State Positive Affect") + ylim(-4,4) +
ggtitle("Within-person Associations")## Warning: Removed 17 rows containing non-finite values (stat_smooth).
## Warning: Removed 17 rows containing non-finite values (stat_smooth).
Note that I have “fudged” this plot. It is not derived from the model, but from the background data files (and individual regression outside the MLM).
- Lets check how the coefficients from the MLM match up with the individual-level regressions.
We obtain estimates of the indivual level esimates from MLM model
#Fixed effect estimates
FE <- fixef(model1.fit)
FE## (Intercept) negaff.state imean.negaff
## 5.1171237 -0.5582807 -0.4096544#Individual random effect estimates
RE <- ranef(model1.fit)
head(RE)## (Intercept) negaff.state
## 101 0.05007425 0.02077294
## 102 0.83375669 -0.13558575
## 103 -0.49398945 -0.25869199
## 104 -0.49180712 -0.11149607
## 105 0.02171479 0.04352667
## 106 0.34837596 -0.41678040#Individual intercepts (MLM model based)
MLM.intercept <- FE[1]+RE[,1]
#Individual slopes (MLM model based)
MLM.negaff.state <- FE[2]+RE[,2]
#Combining the individual regression intercepts and slopes and the model based intercepts and slopes
indiv.parm.combined <- cbind(icoeffs,MLM.intercept,MLM.negaff.state)
#Just note that this is dangerous becasue the lengths may differ, alignment may be off.
head(indiv.parm.combined)## id iintercept.pana islope.pana MLM.intercept MLM.negaff.state
## 1 101 5.304688 -0.4947917 5.167198 -0.5375078
## 2 102 7.185714 -0.9062374 5.950880 -0.6938665
## 3 103 6.294681 -1.1564198 4.623134 -0.8169727
## 4 104 5.328046 -0.9213198 4.625317 -0.6697768
## 5 105 5.091801 -0.3892060 5.138838 -0.5147541
## 6 106 8.306612 -1.7648805 5.465500 -0.9750611Note the “shrinkage”.
Also apparent in the descriptives
describe(indiv.parm.combined)## vars n mean sd median trimmed mad min
## id 1 190 318.29 130.44 321.50 318.99 151.23 101.00
## iintercept.pana 2 190 5.30 1.61 5.34 5.34 1.43 -0.30
## islope.pana 3 187 -0.51 0.63 -0.54 -0.50 0.52 -3.46
## MLM.intercept 4 190 5.12 0.57 5.13 5.11 0.52 3.69
## MLM.negaff.state 5 190 -0.56 0.16 -0.57 -0.55 0.14 -1.08
## max range skew kurtosis se
## id 532.00 431.00 -0.04 -1.09 9.46
## iintercept.pana 9.93 10.23 -0.24 0.74 0.12
## islope.pana 1.80 5.25 -0.42 3.26 0.05
## MLM.intercept 6.56 2.87 0.05 -0.08 0.04
## MLM.negaff.state -0.16 0.92 -0.34 0.44 0.01Looking at the correlations between MLM model based parameters and individual-level parameters
corr.test(indiv.parm.combined[,-1])## Call:corr.test(x = indiv.parm.combined[, -1])
## Correlation matrix
## iintercept.pana islope.pana MLM.intercept
## iintercept.pana 1.00 -0.85 0.50
## islope.pana -0.85 1.00 -0.08
## MLM.intercept 0.50 -0.08 1.00
## MLM.negaff.state -0.81 0.76 -0.27
## MLM.negaff.state
## iintercept.pana -0.81
## islope.pana 0.76
## MLM.intercept -0.27
## MLM.negaff.state 1.00
## Sample Size
## iintercept.pana islope.pana MLM.intercept
## iintercept.pana 190 187 190
## islope.pana 187 187 187
## MLM.intercept 190 187 190
## MLM.negaff.state 190 187 190
## MLM.negaff.state
## iintercept.pana 190
## islope.pana 187
## MLM.intercept 190
## MLM.negaff.state 190
## Probability values (Entries above the diagonal are adjusted for multiple tests.)
## iintercept.pana islope.pana MLM.intercept
## iintercept.pana 0 0.00 0.00
## islope.pana 0 0.00 0.25
## MLM.intercept 0 0.25 0.00
## MLM.negaff.state 0 0.00 0.00
## MLM.negaff.state
## iintercept.pana 0
## islope.pana 0
## MLM.intercept 0
## MLM.negaff.state 0
##
## To see confidence intervals of the correlations, print with the short=FALSE option#plots
pairs.panels(indiv.parm.combined[,-1])
WARNING: It is not always clear which is the “correct” model - the collection of individual-level regressions, or the MLM model. In both cases, post-hoc estimated “random effect scores” are not exact scores. They are estimated scores.
OK - this gets us started with the MLM. Next time we shall go back and take care of some important details.