[R-sig-ME] Correlations among random variables

Uanhoro, James u@nhoro@1 @ending from buckeyem@il@o@u@edu
Tue Jan 15 14:39:07 CET 2019



***********************************************************************************
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

        [[alternative HTML version deleted]]

_______________________________________________
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