[R-sig-ME] Correlations among random variables
Thierry Onkelinx
thierry@onkelinx @ending from inbo@be
Mon Jan 14 13:42:24 CET 2019
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
Havenlaan 88 bus 73, 1000 Brussel
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://www.inbo.be>
Op zo 13 jan. 2019 om 04:42 schreef Avraham Kluger <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>
> 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]]
>
> _______________________________________________
> 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