[R-sig-ME] Correlations among random variables
Avraham Kluger
@vik @ending from @@vion@huji@@c@il
Sun Jan 13 04:42:01 CET 2019
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
> VarCorr(lme4Mlm)
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 using buckeyemail.osu.edu>
Sent: Saturday, January 12, 2019 5:38 PM
To: Avraham Kluger<mailto:avik using 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 using savion.huji.ac.il<mailto:avik using 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 using mail.huji.ac.il<mailto:limor.borut using mail.huji.ac.il>
# written by Avi Kluger: avik using savion.huji.ac.il<mailto:avik using 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 using buckeyemail.osu.edu<mailto:uanhoro.1 using 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 using buckeyemail.osu.edu]
Sent: Saturday, January 12, 2019 4:41 PM
To: Avraham Kluger <avik using savion.huji.ac.il<mailto:avik using 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 using savion.huji.ac.il<mailto:avik using savion.huji.ac.il>> wrote:
Dear James,
As you might have seen in my second message to r-sig-mixed-models using r-project.org<mailto:r-sig-mixed-models using 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 using buckeyemail.osu.edu]
Sent: Saturday, January 12, 2019 3:17 PM
To: Avraham Kluger <avik using savion.huji.ac.il<mailto:avik using savion.huji.ac.il>>
Cc: r-sig-mixed-models using r-project.org<mailto:r-sig-mixed-models using 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 using savion.huji.ac.il<mailto:avik using 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]]
_______________________________________________
R-sig-mixed-models using r-project.org<mailto:R-sig-mixed-models using r-project.org> mailing list
https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
[[alternative HTML version deleted]]
More information about the R-sig-mixed-models
mailing list