Hi,
Following help from James Uanhoro, I produced models with correlated
random variables both with nlms and lme4. Curiously, the estimates of the
variances and their correlation are identical, but the error variances and
their correlation are not. Yet, the sum of the error variances are
identical. For example, in the nlme code below the error for focalcode +
residual is .376 + .046 = .422, and in the lme4 it is .2 + .222 = .422.
Can anyone guide me to read about these different decompositions? This may
explain the different correlations among the error terms .265 vs. .451.
MLM with nlme
Variance StdDev Corr
focalid = pdLogChol(0 + focalcode + partcode)
focalcode 0.20840953 0.4565189 foclcd
partcode 0.06089854 0.2467763 0.699
dyadid = pdLogChol(0 + focalcode + partcode)
focalcode 0.37674500 0.6137956 foclcd
partcode 0.50282362 0.7091006 0.265
Residual 0.04641050 0.2154310
MLM with lme4
as.data.frame(VarCorr(lme4Mlm))
grp var1 var2 vcov sdcor
1 dyadid:focalid focalcode <NA> 0.20018050 0.4474153
2 dyadid:focalid partcode <NA> 0.32625940 0.5711912
3 dyadid:focalid focalcode partcode 0.11523355 0.4509065
4 focalid focalcode <NA> 0.20840930 0.4565187
5 focalid partcode <NA> 0.06089846 0.2467761
6 focalid focalcode partcode 0.07872740 0.6988182
7 Residual <NA> <NA> 0.22297497 0.4722023
Groups Name Std.Dev. Corr
dyadid:focalid focalcode 0.44742
partcode 0.57119 0.451
focalid focalcode 0.45652
partcode 0.24678 0.699
Residual 0.47220
Yours,
Avi Kluger<http://avikluger.wix.com/avi-kluger>
From: Uanhoro, James<mailto:uanhoro.1 at buckeyemail.osu.edu>
Sent: Saturday, January 12, 2019 5:38 PM
To: Avraham Kluger<mailto:avik at savion.huji.ac.il>
Subject: Re: [R-sig-ME] Correlations among random variables
That is the exact lme4 syntax that replicates the nlme model. I did not
attempt to run it but when I did, I noticed that the problem is you have
two random effects, each having 624 values, 624 by 2 equals the sample
size. lme4 will not run when this happens hence the error message. You can
tell lme4 to run nevertheless
summary(mlms <- lmer(
outcome ~ 0 + focalcode + partcode +
(0 + focalcode + partcode | focalid/dyadid),
Chapter10_df,
control = lmerControl(check.nobs.vs.nRE = "warning")))
This will force the program to run the code, and will print out warnings.
The log likelihood was the same with that from nlme indicating that the
models are the same. But the variance-covariance matrices for the random
effects by the interaction between focalid and dyadid - the inner cluster
variable - are different. Which one do you trust? Probably neither - there
is just not enough information to estimate this varCov matrix. See this
thread for commentary on the issue:
https://github.com/lme4/lme4/issues/175
Hope this helps, -James.
On Jan 12 2019, at 9:49 am, Avraham Kluger <avik at savion.huji.ac.il<mailto:
avik at savion.huji.ac.il>> wrote:
Oops,
Actually the code produces error
Error: number of observations (=1248) <= number of random effects (=1248)
for term (0 + focalcode + partcode | dyadid:focalid); the random-effects
parameters and the residual variance (or scale parameter) are probably
unidentifiable
The code below should run on any machine.
Best
Avi
################################################################################
# **************************** R companion for
**************************
#
# Kenny, D. A., Kashy, D. A., & Cook, W. L. (2006). Dyadic data analysis.
# New York: Guilford Press.
#
# lme code developed by Limor Borut: limor.borut at mail.huji.ac.il<mailto:
limor.borut at mail.huji.ac.il>
# written by Avi Kluger: avik at savion.huji.ac.il<mailto:
avik at savion.huji.ac.il>
#
# CHAPTER 10 -- one with many SRM
#
################################################################################
rm(list = ls()) # Clean the Global
Environment
cat ("\014") # Clean the R console
if (is.null(dev.list()) == FALSE) dev.off() # Clean Plots
# Read (in SPSS format) from Kenny's book site and replicate Table 9.1
if (!require('foreign')) install.packages('foreign'); library('foreign')
Chapter10_df <- read.spss("http://davidakenny.net/kkc/c10/c10_recip.sav",
to.data.frame = TRUE, use.value.labels = FALSE)
head(Chapter10_df)
if (!require("nlme")) install.packages("nlme");
suppressMessages(library(nlme))
mlm <- lme(outcome ~ 0 + focalcode + 0 + partcode,
random = ~ 0 + focalcode + partcode|focalid/dyadid,
data = Chapter10_df)
summary(mlm)
intervals(mlm)
mlmOutput <- VarCorr(mlm)
VarCorr(mlm)
cat(
"Actor variance = ", round(as.numeric(VarCorr(mlm)[, "Variance"][3]),
3),
"\nPartner variance = ", round(as.numeric(VarCorr(mlm)[, "Variance"][2]),
3),
"\nGeneralized Reciprocity = ", round(as.numeric(VarCorr(mlm)[,
"Corr"][3]), 3),
"\nDyadic Reciprocity = ", round(as.numeric(VarCorr(mlm)[, "Corr"][6]),
3), "\n"
)
# Very Important Note. The original data coded with 0 the focal person.
# Therefore the first random variable above is partner variance. Reversing
# the codes below make the results more intuitive. I thank David Kenny for
# Clarifying this issue.
Chapter10_df$focalcode <- 1- Chapter10_df$focalcode
Chapter10_df$partcode <- 1- Chapter10_df$partcode
mlm <- lme(outcome ~ 0 + focalcode + 0 + partcode,
random = ~ 0 + focalcode + partcode|focalid/dyadid,
data = Chapter10_df)
VarCorr(mlm)
# An alternative suggsted by James Uanhoro <uanhoro.1 at buckeyemail.osu.edu
<mailto:uanhoro.1 at buckeyemail.osu.edu>>
if (!require("lme4")) install.packages("lme4");
suppressMessages(library(lme4))
mlm <- lmer(outcome ~ 0 + focalcode + partcode + role +
(0 + focalcode + partcode | focalid/ dyadid),
data = Chapter10_df)
summary(mlm)
VarCorr(mlm)
From: Uanhoro, James [mailto:uanhoro.1 at buckeyemail.osu.edu]
Sent: Saturday, January 12, 2019 4:41 PM
To: Avraham Kluger <avik at savion.huji.ac.il<mailto:avik at savion.huji.ac.il>>
Subject: RE: [R-sig-ME] Correlations among random variables
My reply addressed to issues: covariance between error terms; and between
random effects.
The only change to lme4 for random effects is to switch the double pipe to
a single pipe in the random effects specification of the model, as I have
done below:
mlm <- lmer(outcome ~ 0 + focalcode + partcode + role +
(0 + focalcode + partcode | focalid/ dyadid),
data = df)
James.
On Jan 12, 2019 09:05, Avraham Kluger <avik at savion.huji.ac.il<mailto:
avik at savion.huji.ac.il>> wrote:
Dear James,
As you might have seen in my second message to
r-sig-mixed-models at r-project.org<mailto:r-sig-mixed-models at r-project.org>,
my student solved this problem with nlme. Would you know how to write it
in lme4?
Here is the working nlme code
mlm <- lme(outcome ~ 0 + focalcode + 0 + partcode,
random = ~ 0 + focalcode + partcode|focalid/dyadid,
data = Chapter10_df)
Best,
Avi
From: Uanhoro, James [mailto:uanhoro.1 at buckeyemail.osu.edu]
Sent: Saturday, January 12, 2019 3:17 PM
To: Avraham Kluger <avik at savion.huji.ac.il<mailto:avik at savion.huji.ac.il>>
Cc: r-sig-mixed-models at r-project.org<mailto:
r-sig-mixed-models at r-project.org>
Subject: Re: [R-sig-ME] Correlations among random variables
In the lme4 syntax, you'd have to change the double pipe, ||, when
specifying the random effects to a single pipe, |, to permit a correlation
between random effects. lme4 is faster than nlme.
Assuming lme4 and nlme are the only options ... If you want to specify an
error covariance structure beyond the covariance structure implied by
standard multilevel models, you will have to use nlme. nlme has a
`correlation =` argument that allows different covariance structures,
corSymm (general/unstructured), corCompSymm (exchangeable), ...
On Jan 12, 2019 02:01, Avraham Kluger <avik at savion.huji.ac.il<mailto:
avik at savion.huji.ac.il>> wrote:
Hi,
I am struggling to analyze, in R, MLM models that specify correlations
among random variables, as can be done with SPSS, SAS, or MlWin.
Consider the following code in SPSS
-----------------------------
MIXED
Outcome BY role WITH focalcode partcode
/FIXED = focalcode partcode | NOINT
/PRINT = SOLUTION TESTCOV
/RANDOM focalcode partcode | SUBJECT(focalid) COVTYPE(UNR)
/REPEATED = role | SUBJECT(focalid*dyadid) COVTYPE(UNR).
-----------------------------
And a minimal code (with data) in R
-----------------------------
df <- read.csv("
https://raw.githubusercontent.com/avi-kluger/RCompanion4DDABook/master/Chapter%2010/Chapter10_df.csv
")
head(df)
library(lme4)
mlm <- lmer(outcome ~ 0 + focalcode + partcode + role +
(0 + focalcode + partcode|| focalid/ dyadid),
data = df)
summary(mlm)
-----------------------------
These SPSS and R codes produce the same variance estimates. However, SPSS
also produces a correlation among "focalcode" and "partcode." How can this
be done in R? Is it also possible to produce the correlation among the
respective error variances (as in SPSS)?
Additional information
1. MOTIVATION. The question arises from David Kenny's work on
one-with-many reciprocal designs (e.g., a manager rate all subordinates,
and all subordinates rate the same manager). These models estimate the
variance stemming from the one (e.g., managers) and the many (e.g.,
subordinates), and the correlation among them (termed generalized
reciprocity). The data and codes for SAS etc. are available at
http://davidakenny.net/kkc/c10/c10.htm.
2. SPSS OUTPUT (download HTML file):
https://www.dropbox.com/s/eqch0kq6djtbsfx/One%20with%20many%20SPSS%20output.htm?dl=1
Sincerely,
Avi Kluger
https://www.avi-kluger.com/
[[alternative HTML version deleted]]