[R-sig-ME] Correlations among random variables
Avraham Kluger
@v|k @end|ng |rom @@v|on@huj|@@c@||
Thu Jan 24 05:39:43 CET 2019
Hi again,
A reminder, I am trying to get a stable estimate of the correlation between two random errors in a mode like:
lme(outcome ~ 0 + focalcode + partcode, random = ~ 0 + focalcode + partcode|focalid/dyadid, data = x)
I figured out that five different methods yield identical covariance estimate, as I show below. My new question is CAN YOU CONSTRAIN RESIDUAL TO ZERO in either lme4 or nlme? This may allow convergence and hence estimate of standard error for the covariance.
-------------------------------------------------------------------------------------------------------------------------------------------------
Method Var Focal Var Partner Residual Cov Corr
--------------------------------------------------------------------------------------------------------------------------------------------------
SPSS .549 .423 .239
lavaan (SEM) .552 .425 .116 .239
lme4 .326 .200 .223 .115 .451
nlme .376 .502 .046 .265
lavaan (SEM) with latent residuals .331 .457 .094 .116 .298
-------------------------------------------------------------------------------------------------------------------------------------------------
Some methods printout the covariance and some do not (or fail to converge properly), but all suggest that the covariance is .115 or .116. For MLM in R:
lme4, cov = .451 * sqrt(.326 * .200) = .115
nlme, cov = .265 * sqrt(.376 * .502) = .115
for SPSS that runs without any warning, cov = .239 * sqrt(.549 * .423) = .115
I appreciate your patience with a novice,
Best,
Avi
-----Original Message-----
From: Avraham Kluger
Sent: Wednesday, January 16, 2019 12:27 PM
To: Uanhoro, James <uanhoro.1 using buckeyemail.osu.edu>; R-sig-mixed-models using r-project.org
Cc: Kenny, David <david.kenny using uconn.edu>
Subject: RE: [R-sig-ME] Correlations among random variables
Hi,
I thank James, Wolfgang , Thierry , and David for educating me. To summarize all responses, you seem to suggest: Do not trust your data when it comes to the correlation between the error terms. I still have some questions about intepreration. Note that I do not care particularity about these data, but I am trying to grasp the principle. First, I repeated the calculations with SPSS and with SEM (in the later using wide-data formal, and placing equality constraints within Focal and within Partner on intercepts, variances, and covariances). The summary of all the results can be found below.
Error variances, covariances, and correlations by software/analytic approach
---------------------------------------------------------------------------------------------------------------------------------
Var Focal Var Partner Residual Total Var Focal Total Var Partner Cov (SE/CI) Corr (SE/CI)
---------------------------------------------------------------------------------------------------------------------------------
SPSS .549 .423 NA .549 .423 NA .239 (.046)
lme4 .326 .200 .223 .549 .423 0.115(*) .451
nlme .376 .502 .046 .549 .423 (*) .265
lavaan (SEM) .331 .457 .094 .551 .425 0.116 (.024) .298
---------------------------------------------------------------------------------------------------------------------------------
* cannot get confidence intervals on var-cov components: Non-positive definite approximate variance-covariance
1. How are the correlations calculated. Assuming that the covariance is .115 (or .116), what are the variances that their product serve as a denominator in transforming the covariance into the correlation?
2. David Dupphy wrote:
FWIW, the raw correlations (that are being analysed) when reordered as a a bivariate setup:
partner-focal r=0.31 [dyads]
focal intraclass r=0.10 (jack SE=0.06) [clustered on focal] partner intraclass r=0.33 (jSE=0.07) [clustered on focal]
Given that all analyses suggest a POSITIVE covariance, would it be reasonable to conclude that there is a positive correlation in the population, although its magnitude is uncertain? Or would you still believe these signals are noise?
I with deep gratitude to this form and all the help received,
Avi
From: Uanhoro, James [uanhoro.1 using buckeyemail.osu.edu]
Sent: Tuesday, January 15, 2019 3:39 PM
To: Avraham Kluger
Cc: R-sig-mixed-models using r-project.org
Subject: Re: [R-sig-ME] Correlations among random variables
***********************************************************************************
This email contains links to a hosting company that is frequented by scammers.
Please think before you click any links.
The Ohio State University will NEVER ask you for your account information by email.
If you receive such a message, please report it to report-phish using osu.edu
NEVER reply to any email asking you for your account information
or other personal details. For more information or to get help,
contact the IT Service Desk by calling 614-688-HELP (4357).
***********************************************************************************
A correlation always exists in the population, if we assume the realistic position that exactly nil effects are just not true.
The points Wolfgang Viechtbauer makes are very important. That is also the crux of this discussion: https://github.com/lme4/lme4/issues/175#issuecomment-33580591.
If your question pertains to the variance covariance matrix that changes from lme4 to nlme, then you're probably asking too much of the data. And the only information you should trust about this particular variance covariance matrix is that you should trust nothing else about it.
James.
On Jan 15, 2019 00:58, Avraham Kluger <avik using savion.huji.ac.il> wrote:
Dear Thierry,
I thank you for the reference to your blog. It does raise questions abo ut noise in the data. However, if the correlations are significant in any method, would you conclude that a correlation exists in the population?
Avi
From: Thierry Onkelinx<mailto:thierry.onkelinx using inbo.be>
Sent: Monday, January 14, 2019 2:42 PM
To: Avraham Kluger<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>; Kenny, David<mailto:david.kenny using uconn.edu>
Subject: Re: [R-sig-ME] Correlations among random variables
Dear Avraham,
Do you have a huge amount of random effects? If not, the variance estimates have a large uncertainty. So that you precieve as a strong diverence is actually just noise from the model uncertainty. I wrote a small blog post on the number of random effect levels and the resulting uncertainty on the variance estimates: https://www.muscardinus.be/2018/09/number-random-effect-levels/
Best regards,
ir. Thierry Onkelinx
Statisticus / Statistician
Vlaamse Overheid / Government of Flanders
INSTITUUT VOOR NATUUR- EN BOSONDERZOEK / RESEARCH INSTITUTE FOR NATURE AND FOREST
Team Biometrie & Kwaliteitszorg / Team Biometrics & Quality Assurance
thierry.onkelinx using inbo.be<mailto:thierry.onkelinx using inbo.be>
Havenlaan 88 bus 73, 1000 Brussel
www.inbo.be<http://www.inbo.be>
///////////////////////////////////////////////////////////////////////////////////////////
To call in the statistician after the experiment is done may be no more than asking him to perform a post-mortem examination: he may be able to say what the experiment died of. ~ Sir Ronald Aylmer Fisher
The plural of anecdote is not data. ~ Roger Brinner
The combination of some data and an aching desire for an answer does not ensure that a reasonable answer can be extracted from a given body of data. ~ John Tukey
///////////////////////////////////////////////////////////////////////////////////////////
[https://inbo-website-prd-532750756126.s3-eu-west-1.amazonaws.com/inbologoleeuw_nl.png]<https://www.inbo.be>
Op zo 13 jan. 2019 om 04:42 schreef Avraham Kluger <avik using savion.huji.ac.il<mailto:avik using savion.huji.ac.il>>:
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<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<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><mailto: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><mailto: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><mailto: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><mailto: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<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><mailto: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><mailto: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><mailto: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<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><mailto: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><mailto: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><mailto: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><mailto: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]]
_______________________________________________
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
d
[[alternative HTML version deleted]]
_______________________________________________
R-sig-mixed-models using r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
More information about the R-sig-mixed-models
mailing list