Kornbrot, Diana
2018-11-25 17:55:40 UTC
Error means squares in GLMER and LMER (Kornbrot, Diana)
Have now got lmer in R ad MIXED in SPSS to agree for my data in s4 a 2 between, 2 within anova, iwht pno as random subjects and freq as dependent
You are all correct do not nee ls means directly BUT fa values have numerators and denominators and would really like to see both
Rscript
a<- lmer(freq~b1*b2*w1*w2+(w1|pno) + (w2|pno), data=s4)
anova(a)for inferential and ls_mean(a) for descriptive means.
want to go on to ‘correct' analysis which has a binomial proprtion of freq/Nmax
so tried
b<- glmer(cbind(freq, Nmax-freq) ~ b1*b2*w1*w2 +(w1|pno)+(w2|pno), data= s4,family="binomial" )
Warning message:
In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
Model failed to converge with max|grad| = 0.00692938 (tol = 0.001, component 1)
So, what can I change so it does converge. Its by no means a big data set
it was happy with print(b), summary(b) and anova(b) but results do not agree with SPSS, which did converge
BUT would not give me mea
no applicable method for 'ls_means' applied to an object of class "c('glmerMod', 'merMod')"
So how does one get means from objects of class "c('glmerMod', 'merMod’)”?
Also when I try to look up glmer in packages, it says it ca’t be installed in R3.51. although it must be installed to obey the command. Mystified
All help gratefully received
best
Diana
Message: 3
Date: Thu, 22 Nov 2018 17:09:56 +0000
From: "Kornbrot, Diana" <***@herts.ac.uk<mailto:***@herts.ac.uk>>
To: roee maor <***@gmail.com<mailto:***@gmail.com>>
Cc: "r-sig-mixed-***@r-project.org<mailto:r-sig-mixed-***@r-project.org>"
<r-sig-mixed-***@r-project.org<mailto:r-sig-mixed-***@r-project.org>>
Subject: Re: [R-sig-ME] Error means squares in GLMER and LMER
Message-ID: <8D3FBC0E-8BFD-4AC9-89BF-***@herts.ac.uk<mailto:8D3FBC0E-8BFD-4AC9-89BF-***@herts.ac.uk>>
Content-Type: text/plain; charset="utf-8"
So I am comparing standard ANOVA on raw frequencies (or equivalently probabilities) with GLMM for binomial proportion with both logit and probit link
ALL analyses have been completed in SPSS using MIXED for: response = identity , link = normal ; response = proportion=freq/Nmax, link = probit; response = proportion link = logit)
I want to show how to do identical analyses in R using lmer for are freq and glmer for proportions
So I wasn’t SAME results from R and SPSS (and a diamond necklace for Christmas, celebrated as an EU citizen in UK - I am a demanding woman)
Results are NOT quite the same.
I am checking using the raw probabilities as raw response with lmer before moving on to glmer for proportions
Check 1, in SPSS for raw freq or probability REPEATED give same result as MIXED (response = identity, link = normal). where there are differences it is REPEATED WITHIN comparisons not MULTIVARIATE.
Check 2. Compare R, lmer with SPSS mixed
if repeated groups are w1, w2, etc and between groups are b1, b2 etc, I use:
result <- lmer(freq~b1*b2*w1*w2 + (w1|subject) + (w2|subject), data = test)
anova(result)
Fand df from R and SPSS do not always match, even when they do match on sums of squares
I am trying to work out WHY there is a mismatch
Thought that knowing what was in the DENOMINATOR of the F values - which i perhaps wrongly termed error sums of squares, might help
I want F for usual reasons: to test significance and estimate effect size.
I also want all my packages to give me the SAME F and df2 and to UNDERSTAND what is happening
Sorry this is so long, but hope it is now clearer
best
Diana
and as an extra treat would like marginal means from object of type lmer
Dear Diana,
If indeed what you're looking for is what Rolf mentioned, you might find Nakagawa & Schielzeth (2013) helpful.
It's titled "A general and simple method for obtaining R2 from generalized linear mixed-effects models". There is no dispute about the method's generality, but simple is a relative term...
Here's the link: https://besjournals.onlinelibrary.wiley.com/doi/epdf/10.1111/j.2041-210x.2012.00261.x
Hope this helps,
--
Roi Maor
PhD candidate
School of Zoology, Tel Aviv University
Centre for Biodiversity and Environment Research, UCL
_____________________________________
Professor Diana Kornbrot
Mobile
+44 (0) 7403 18 16 12
Work
University of Hertfordshire
College Lane, Hatfield, Hertfordshire AL10 9AB, UK
+44 (0) 170 728 4626
***@herts.ac.uk<mailto:***@herts.ac.uk><mailto:***@herts.ac.uk>
http://dianakornbrot.wordpress.com/
http://go.herts.ac.uk/Diana_Kornbrot
skype: kornbrotme
Home
19 Elmhurst Avenue
London N2 0LT, UK
+44 (0) 208 444 2081
------------------------------------------------------------
[[alternative HTML version deleted]]
------------------------------
Message: 4
Date: Thu, 22 Nov 2018 14:44:27 -0500
From: Ben Bolker <***@gmail.com>
To: "Kornbrot, Diana" <***@herts.ac.uk>
Cc: roee maor <***@gmail.com>, R SIG Mixed Models
<r-sig-mixed-***@r-project.org>
Subject: Re: [R-sig-ME] Error means squares in GLMER and LMER
Message-ID:
<***@mail.gmail.com>
Content-Type: text/plain; charset="utf-8"
For marginal means, use the emmeans package.
If you use the lmerTest package, you can get Satterthwaite or
Kenward-Roger df: you can use lme (from the nlme package), or
<https://github.com/bbolker/mixedmodels-misc/blob/master/R/calcDenDF.R>,
to get df via a simple "parameter-counting" exercise.
The problem is that the "F statistics" are quite poorly defined for
GLMMs. Can you show us the contrasting results you're getting for SPSS
and glmer? Do you know how SPSS is computing the F statistics? (This
<https://www.ibm.com/support/knowledgecenter/en/SS3RA7_15.0.0/com.ibm.spss.modeler.help/idh_glmm_build_options.htm>
makes it seem like it might be using Satterthwaite approximations ...)
On Thu, Nov 22, 2018 at 12:09 PM Kornbrot, Diana
<***@herts.ac.uk> wrote:
So I am comparing standard ANOVA on raw frequencies (or equivalently probabilities) with GLMM for binomial proportion with both logit and probit link
ALL analyses have been completed in SPSS using MIXED for: response = identity , link = normal ; response = proportion=freq/Nmax, link = probit; response = proportion link = logit)
I want to show how to do identical analyses in R using lmer for are freq and glmer for proportions
So I wasn’t SAME results from R and SPSS (and a diamond necklace for Christmas, celebrated as an EU citizen in UK - I am a demanding woman)
Results are NOT quite the same.
I am checking using the raw probabilities as raw response with lmer before moving on to glmer for proportions
Check 1, in SPSS for raw freq or probability REPEATED give same result as MIXED (response = identity, link = normal). where there are differences it is REPEATED WITHIN comparisons not MULTIVARIATE.
Check 2. Compare R, lmer with SPSS mixed
if repeated groups are w1, w2, etc and between groups are b1, b2 etc, I use:
result <- lmer(freq~b1*b2*w1*w2 + (w1|subject) + (w2|subject), data = test)
anova(result)
Fand df from R and SPSS do not always match, even when they do match on sums of squares
I am trying to work out WHY there is a mismatch
Thought that knowing what was in the DENOMINATOR of the F values - which i perhaps wrongly termed error sums of squares, might help
I want F for usual reasons: to test significance and estimate effect size.
I also want all my packages to give me the SAME F and df2 and to UNDERSTAND what is happening
Sorry this is so long, but hope it is now clearer
best
Diana
and as an extra treat would like marginal means from object of type lmer
Dear Diana,
If indeed what you're looking for is what Rolf mentioned, you might find Nakagawa & Schielzeth (2013) helpful.
It's titled "A general and simple method for obtaining R2 from generalized linear mixed-effects models". There is no dispute about the method's generality, but simple is a relative term...
Here's the link: https://besjournals.onlinelibrary.wiley.com/doi/epdf/10.1111/j.2041-210x.2012.00261.x
Hope this helps,
--
Roi Maor
PhD candidate
School of Zoology, Tel Aviv University
Centre for Biodiversity and Environment Research, UCL
_____________________________________
Professor Diana Kornbrot
Mobile
+44 (0) 7403 18 16 12
Work
University of Hertfordshire
College Lane, Hatfield, Hertfordshire AL10 9AB, UK
+44 (0) 170 728 4626
***@herts.ac.uk<mailto:***@herts.ac.uk>
http://dianakornbrot.wordpress.com/
http://go.herts.ac.uk/Diana_Kornbrot
skype: kornbrotme
Home
19 Elmhurst Avenue
London N2 0LT, UK
+44 (0) 208 444 2081
------------------------------------------------------------
[[alternative HTML version deleted]]
_______________________________________________
R-sig-mixed-***@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
------------------------------
Message: 5
Date: Thu, 22 Nov 2018 22:59:33 +0000
From: roee maor <***@gmail.com>
To: ***@ed.ac.uk
Cc: R-sig-mixed-***@r-project.org
Subject: Re: [R-sig-ME] Unclear output from MCMCglmm with categorical
predictors
Message-ID:
<CACxNx6uqb0C3HCEpvOfKDF=pOKTzXjgZStFS-***@mail.gmail.com>
Content-Type: text/plain; charset="utf-8"
Dear Jarrod and Ben,
Thank you both for the advice.
Apologies if I accidentally caused anyone to receive an unexpected email.
Running the model specified above returns the error: "non-reduced nodes do
not appear first". The error persists when using a fully-bifurcating tree
so the issue is not tree polytomies.
I tried to go through the source code line-by-line to locate the problem
but got lost in the many loops and conditional statements. This error
message is conditioned on a match() output that also returns a warning that
the two arguments matched are unequal in length: one corresponds to the
tips while the other includes the tree's internal nodes as well.
As always, any help is much appreciated.
Best,
Roi
On Wed, 21 Nov 2018 at 15:05, HADFIELD Jarrod <***@ed.ac.uk> wrote:
Hi,
You could upload the tarball to winbuilder (
https://win-builder.r-project.org/) and build a Windows source package.
Cheers,
Jarrod
On 21/11/2018 14:09, roee maor wrote:
Hi Jarrod,
Many thanks for your reply.
I couldn't install the tarball on R v3.4.3 or v3.5.1 so sourced the files
directly to the workspace.
I tried to run this model as you suggested:
T1 <- MCMCglmm(Activity ~ -1 + log(Mass) + Max.Temp * Annual.Precip,
random = ~ animal,
prior = list(R = list(fix=1, V=1e-15), G = list(G1 =
list(V=1, nu=0.002))),
pedigree = datatree,
reduced = TRUE,
burnin = 50000, nitt = 750001, thin = 700,
family = "threshold",
data = Rdata,
pl = TRUE, saveX = TRUE, saveZ = TRUE,
verbose = TRUE)
It returns some errors about missing functions "is.positive.definite" and
"Matrix", which I addressed with:
library("corpcor", lib.loc="~/R/win-library/3.5")
library("MatrixModels", lib.loc="~/R/win-library/3.5")
but I can't figure this one out:
'Error in .C("MCMCglmm", as.double(data$MCMC_y),
as.double(data$MCMC_y.additional), :
C symbol name "MCMCglmm" not in load table'
Detaching these packages doesn't necessarily cause that same error to
appear although I execute the exact same code.
Also, several attempts (same code again) caused a fatal error and
automatic session termination (info for a similar session below if
interesting).
I tried to use the fully bifurcating tree as an experiment but that made
no difference.
Any ideas what this last error means?
Thanks!
Roi
sessionInfo()
R version 3.5.1 (2018-07-02)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows >= 8 x64 (build 9200)
Matrix products: default
locale:
[1] LC_COLLATE=English_United Kingdom.1252 LC_CTYPE=English_United
Kingdom.1252 LC_MONETARY=English_United Kingdom.1252
[4] LC_NUMERIC=C LC_TIME=English_United
Kingdom.1252
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] corpcor_1.6.9 Matrix_1.2-14 phytools_0.6-60 maps_3.3.0
ape_5.2
loaded via a namespace (and not attached):
[1] igraph_1.2.2 Rcpp_0.12.19 magrittr_1.5
MASS_7.3-50 mnormt_1.5-5
[6] scatterplot3d_0.3-41 lattice_0.20-35 quadprog_1.5-5
fastmatch_1.1-0 tools_3.5.1
[11] parallel_3.5.1 grid_3.5.1 nlme_3.1-137
clusterGeneration_1.3.4 phangorn_2.4.0
[16] plotrix_3.7-4 coda_0.19-2 yaml_2.2.0
numDeriv_2016.8-1 animation_2.5
[21] compiler_3.5.1 combinat_0.0-8 expm_0.999-3
pkgconfig_2.0.2
On Tue, 20 Nov 2018 at 20:01, HADFIELD Jarrod <***@ed.ac.uk> wrote:
Hi,
Most likely the phylogenetic heritability (the phylogenetic variance /
the phylogenetic +residual variance) is approaching one resulting in
numerical difficulties. Probably the best thing is to assume that the
phylogenetic heritability equals 1 and use the reduced phylogenetic mixed
model implementation. This allows the phylogenetic heritability to be equal
to 1 without causing numerical issues. At some point I will integrate these
models into the main MCMCglmm package, but for now you can download it from
here: http://jarrod.bio.ed.ac.uk/MCMCglmmRAM_2.24.tar.gz.
Change the name of the ‘’Binomial’ column to ‘animal’ and fit:
T1 <- MCMCglmm(Activity ~ -1 + log(Mass) + Max.Temp * Annual.Precip,
random = ~ animal
prior = list(R = list(fix=1, V=1e-15), G =
list(G1 = list(V=1, nu=0.002))),
pedigree = tree,
reduced=TRUE,
burnin = 150000, nitt = 2650001, thin =
2500,
family = "threshold", data = Tdata,
pl = TRUE, pr = TRUE, saveX = TRUE, saveZ
= TRUE,
verbose = FALSE)
You should need fewer iterations.
Cheers,
Jarrod
On 20 Nov 2018, at 18:08, roee maor <***@gmail.com> wrote:
Dear Jarrod (and list),
Following your previous comment I added "random = ~ Binomial" to my model
to allow for a phylogenetic analysis.
This causes convergence problems: the trace plots show increasing
oscillations along each chain (although no directional trends, so it's not
a burn-in issue). Also, the posterior samples are highly correlated,
residual variance estimates are >10^3 and threshold estimates are high (>20
on the latent scale).
Surprisingly (to me), predictors that are strongly significant in the
non-phylogenetic model lose their effect in the phylogenetic model (I tried
several alternative parameter configurations).
It seems that this model attributes the explained variance to phylogeny
alone.
Can anyone explain what is going on here? Am I specifying the model
poorly or just asking my data more than it can answer?
I tried to overcome this issue by using a fully resolved variant of the
phylogeny, which only improved things slightly.
I also changed the random effect to "random=~Family" or "random=~Order",
which reduced the variance and threshold estimates to more acceptable
levels (<10), but still no significant predictors (and I'm not sure how the
algorithm calculates covariance between higher taxa in the phylogeny).
Separately I tried parameter expanded prior: "prior = list(R = list(V=1,
fix=1), G = list(G1 = list(V=1, nu=1, alpha.mu=0, alpha.V=1000)))". That
didn't help, and messing with priors for this reason feels like poor
practice.
This is the model:
T1 <- MCMCglmm(Activity ~ -1 + log(Mass) + Max.Temp * Annual.Precip,
random = ~ Binomial,
prior = list(R = list(fix=1, V=1), G =
list(G1 = list(V=1, nu=0.002))),
ginverse = list(Binomial=INphylo$Ainv),
burnin = 150000, nitt = 2650001, thin =
2500,
family = "threshold", data = Tdata,
pl = TRUE, pr = TRUE, saveX = TRUE, saveZ
= TRUE,
verbose = FALSE)
The data I use looks like this (not all variables appear in each model):
str(Tdata)
'data.frame': 1389 obs. of 10 variables:
$ Binomial : Factor w/ 1421 levels "Abrocoma_bennettii",..: 1
2 3 4 5 6 7 8 9 10 ...
$ Order : Factor w/ 27 levels "Afrosoricida",..: 24 24 24
24 24 3 24 24 2 2 ...
$ Family : Factor w/ 126 levels "Abrocomidae",..: 1 26 26
26 26 46 74 87 10 10 ...
$ Activity : Factor w/ 3 levels "1","2","3": 1 3 2 2 2 3 2
3 1 3 ...
$ Habitat : Factor w/ 6 levels "Aqua","Arbo",..: 5 5 5 5 5
5 5 5 5 5 ...
$ Diet : Factor w/ 3 levels "Faun","Herb",..: 2 3 3 3
3 1 3 2 2 2 ...
$ Mass : num 250.5 24.9 34.5 38.9 24.5 ...
$ Max.Temp : num 22 16.6 19.1 19.8 17.2 ...
$ Annual.Precip : num 166 645 558 903 1665 ...
Any advice would be much appreciated!
Many thanks,
--
Roi Maor
PhD candidate
School of Zoology, Tel Aviv University
Centre for Biodiversity and Environment Research, UCL
The University of Edinburgh is a charitable body, registered in
Scotland, with registration number SC005336.
--
Roi Maor
PhD candidate
School of Zoology, Tel Aviv University
Centre for Biodiversity and Environment Research, UCL
The University of Edinburgh is a charitable body, registered in
Scotland, with registration number SC005336.
--
Roi Maor
PhD candidate
School of Zoology, Tel Aviv University
Centre for Biodiversity and Environment Research, UCL
[[alternative HTML version deleted]]
------------------------------
Subject: Digest Footer
_______________________________________________
R-sig-mixed-models mailing list
R-sig-mixed-***@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
------------------------------
End of R-sig-mixed-models Digest, Vol 143, Issue 33
***************************************************
_____________________________________
Professor Diana Kornbrot
Mobile
+44 (0) 7403 18 16 12
Work
University of Hertfordshire
College Lane, Hatfield, Hertfordshire AL10 9AB, UK
+44 (0) 170 728 4626
***@herts.ac.uk<mailto:***@herts.ac.uk>
http://dianakornbrot.wordpress.com/
http://go.herts.ac.uk/Diana_Kornbrot
skype: kornbrotme
Home
19 Elmhurst Avenue
London N2 0LT, UK
+44 (0) 208 444 2081
------------------------------------------------------------
[[alternative HTML version deleted]]
Have now got lmer in R ad MIXED in SPSS to agree for my data in s4 a 2 between, 2 within anova, iwht pno as random subjects and freq as dependent
You are all correct do not nee ls means directly BUT fa values have numerators and denominators and would really like to see both
Rscript
a<- lmer(freq~b1*b2*w1*w2+(w1|pno) + (w2|pno), data=s4)
anova(a)for inferential and ls_mean(a) for descriptive means.
want to go on to ‘correct' analysis which has a binomial proprtion of freq/Nmax
so tried
b<- glmer(cbind(freq, Nmax-freq) ~ b1*b2*w1*w2 +(w1|pno)+(w2|pno), data= s4,family="binomial" )
Warning message:
In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
Model failed to converge with max|grad| = 0.00692938 (tol = 0.001, component 1)
So, what can I change so it does converge. Its by no means a big data set
it was happy with print(b), summary(b) and anova(b) but results do not agree with SPSS, which did converge
BUT would not give me mea
ls_means(b)
Error in UseMethod("ls_means") :no applicable method for 'ls_means' applied to an object of class "c('glmerMod', 'merMod')"
So how does one get means from objects of class "c('glmerMod', 'merMod’)”?
Also when I try to look up glmer in packages, it says it ca’t be installed in R3.51. although it must be installed to obey the command. Mystified
All help gratefully received
best
Diana
Message: 3
Date: Thu, 22 Nov 2018 17:09:56 +0000
From: "Kornbrot, Diana" <***@herts.ac.uk<mailto:***@herts.ac.uk>>
To: roee maor <***@gmail.com<mailto:***@gmail.com>>
Cc: "r-sig-mixed-***@r-project.org<mailto:r-sig-mixed-***@r-project.org>"
<r-sig-mixed-***@r-project.org<mailto:r-sig-mixed-***@r-project.org>>
Subject: Re: [R-sig-ME] Error means squares in GLMER and LMER
Message-ID: <8D3FBC0E-8BFD-4AC9-89BF-***@herts.ac.uk<mailto:8D3FBC0E-8BFD-4AC9-89BF-***@herts.ac.uk>>
Content-Type: text/plain; charset="utf-8"
So I am comparing standard ANOVA on raw frequencies (or equivalently probabilities) with GLMM for binomial proportion with both logit and probit link
ALL analyses have been completed in SPSS using MIXED for: response = identity , link = normal ; response = proportion=freq/Nmax, link = probit; response = proportion link = logit)
I want to show how to do identical analyses in R using lmer for are freq and glmer for proportions
So I wasn’t SAME results from R and SPSS (and a diamond necklace for Christmas, celebrated as an EU citizen in UK - I am a demanding woman)
Results are NOT quite the same.
I am checking using the raw probabilities as raw response with lmer before moving on to glmer for proportions
Check 1, in SPSS for raw freq or probability REPEATED give same result as MIXED (response = identity, link = normal). where there are differences it is REPEATED WITHIN comparisons not MULTIVARIATE.
Check 2. Compare R, lmer with SPSS mixed
if repeated groups are w1, w2, etc and between groups are b1, b2 etc, I use:
result <- lmer(freq~b1*b2*w1*w2 + (w1|subject) + (w2|subject), data = test)
anova(result)
Fand df from R and SPSS do not always match, even when they do match on sums of squares
I am trying to work out WHY there is a mismatch
Thought that knowing what was in the DENOMINATOR of the F values - which i perhaps wrongly termed error sums of squares, might help
I want F for usual reasons: to test significance and estimate effect size.
I also want all my packages to give me the SAME F and df2 and to UNDERSTAND what is happening
Sorry this is so long, but hope it is now clearer
best
Diana
and as an extra treat would like marginal means from object of type lmer
Dear Diana,
If indeed what you're looking for is what Rolf mentioned, you might find Nakagawa & Schielzeth (2013) helpful.
It's titled "A general and simple method for obtaining R2 from generalized linear mixed-effects models". There is no dispute about the method's generality, but simple is a relative term...
Here's the link: https://besjournals.onlinelibrary.wiley.com/doi/epdf/10.1111/j.2041-210x.2012.00261.x
Hope this helps,
--
Roi Maor
PhD candidate
School of Zoology, Tel Aviv University
Centre for Biodiversity and Environment Research, UCL
_____________________________________
Professor Diana Kornbrot
Mobile
+44 (0) 7403 18 16 12
Work
University of Hertfordshire
College Lane, Hatfield, Hertfordshire AL10 9AB, UK
+44 (0) 170 728 4626
***@herts.ac.uk<mailto:***@herts.ac.uk><mailto:***@herts.ac.uk>
http://dianakornbrot.wordpress.com/
http://go.herts.ac.uk/Diana_Kornbrot
skype: kornbrotme
Home
19 Elmhurst Avenue
London N2 0LT, UK
+44 (0) 208 444 2081
------------------------------------------------------------
[[alternative HTML version deleted]]
------------------------------
Message: 4
Date: Thu, 22 Nov 2018 14:44:27 -0500
From: Ben Bolker <***@gmail.com>
To: "Kornbrot, Diana" <***@herts.ac.uk>
Cc: roee maor <***@gmail.com>, R SIG Mixed Models
<r-sig-mixed-***@r-project.org>
Subject: Re: [R-sig-ME] Error means squares in GLMER and LMER
Message-ID:
<***@mail.gmail.com>
Content-Type: text/plain; charset="utf-8"
For marginal means, use the emmeans package.
If you use the lmerTest package, you can get Satterthwaite or
Kenward-Roger df: you can use lme (from the nlme package), or
<https://github.com/bbolker/mixedmodels-misc/blob/master/R/calcDenDF.R>,
to get df via a simple "parameter-counting" exercise.
The problem is that the "F statistics" are quite poorly defined for
GLMMs. Can you show us the contrasting results you're getting for SPSS
and glmer? Do you know how SPSS is computing the F statistics? (This
<https://www.ibm.com/support/knowledgecenter/en/SS3RA7_15.0.0/com.ibm.spss.modeler.help/idh_glmm_build_options.htm>
makes it seem like it might be using Satterthwaite approximations ...)
On Thu, Nov 22, 2018 at 12:09 PM Kornbrot, Diana
<***@herts.ac.uk> wrote:
So I am comparing standard ANOVA on raw frequencies (or equivalently probabilities) with GLMM for binomial proportion with both logit and probit link
ALL analyses have been completed in SPSS using MIXED for: response = identity , link = normal ; response = proportion=freq/Nmax, link = probit; response = proportion link = logit)
I want to show how to do identical analyses in R using lmer for are freq and glmer for proportions
So I wasn’t SAME results from R and SPSS (and a diamond necklace for Christmas, celebrated as an EU citizen in UK - I am a demanding woman)
Results are NOT quite the same.
I am checking using the raw probabilities as raw response with lmer before moving on to glmer for proportions
Check 1, in SPSS for raw freq or probability REPEATED give same result as MIXED (response = identity, link = normal). where there are differences it is REPEATED WITHIN comparisons not MULTIVARIATE.
Check 2. Compare R, lmer with SPSS mixed
if repeated groups are w1, w2, etc and between groups are b1, b2 etc, I use:
result <- lmer(freq~b1*b2*w1*w2 + (w1|subject) + (w2|subject), data = test)
anova(result)
Fand df from R and SPSS do not always match, even when they do match on sums of squares
I am trying to work out WHY there is a mismatch
Thought that knowing what was in the DENOMINATOR of the F values - which i perhaps wrongly termed error sums of squares, might help
I want F for usual reasons: to test significance and estimate effect size.
I also want all my packages to give me the SAME F and df2 and to UNDERSTAND what is happening
Sorry this is so long, but hope it is now clearer
best
Diana
and as an extra treat would like marginal means from object of type lmer
Dear Diana,
If indeed what you're looking for is what Rolf mentioned, you might find Nakagawa & Schielzeth (2013) helpful.
It's titled "A general and simple method for obtaining R2 from generalized linear mixed-effects models". There is no dispute about the method's generality, but simple is a relative term...
Here's the link: https://besjournals.onlinelibrary.wiley.com/doi/epdf/10.1111/j.2041-210x.2012.00261.x
Hope this helps,
--
Roi Maor
PhD candidate
School of Zoology, Tel Aviv University
Centre for Biodiversity and Environment Research, UCL
_____________________________________
Professor Diana Kornbrot
Mobile
+44 (0) 7403 18 16 12
Work
University of Hertfordshire
College Lane, Hatfield, Hertfordshire AL10 9AB, UK
+44 (0) 170 728 4626
***@herts.ac.uk<mailto:***@herts.ac.uk>
http://dianakornbrot.wordpress.com/
http://go.herts.ac.uk/Diana_Kornbrot
skype: kornbrotme
Home
19 Elmhurst Avenue
London N2 0LT, UK
+44 (0) 208 444 2081
------------------------------------------------------------
[[alternative HTML version deleted]]
_______________________________________________
R-sig-mixed-***@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
------------------------------
Message: 5
Date: Thu, 22 Nov 2018 22:59:33 +0000
From: roee maor <***@gmail.com>
To: ***@ed.ac.uk
Cc: R-sig-mixed-***@r-project.org
Subject: Re: [R-sig-ME] Unclear output from MCMCglmm with categorical
predictors
Message-ID:
<CACxNx6uqb0C3HCEpvOfKDF=pOKTzXjgZStFS-***@mail.gmail.com>
Content-Type: text/plain; charset="utf-8"
Dear Jarrod and Ben,
Thank you both for the advice.
Apologies if I accidentally caused anyone to receive an unexpected email.
Running the model specified above returns the error: "non-reduced nodes do
not appear first". The error persists when using a fully-bifurcating tree
so the issue is not tree polytomies.
I tried to go through the source code line-by-line to locate the problem
but got lost in the many loops and conditional statements. This error
message is conditioned on a match() output that also returns a warning that
the two arguments matched are unequal in length: one corresponds to the
tips while the other includes the tree's internal nodes as well.
As always, any help is much appreciated.
Best,
Roi
On Wed, 21 Nov 2018 at 15:05, HADFIELD Jarrod <***@ed.ac.uk> wrote:
Hi,
You could upload the tarball to winbuilder (
https://win-builder.r-project.org/) and build a Windows source package.
Cheers,
Jarrod
On 21/11/2018 14:09, roee maor wrote:
Hi Jarrod,
Many thanks for your reply.
I couldn't install the tarball on R v3.4.3 or v3.5.1 so sourced the files
directly to the workspace.
I tried to run this model as you suggested:
T1 <- MCMCglmm(Activity ~ -1 + log(Mass) + Max.Temp * Annual.Precip,
random = ~ animal,
prior = list(R = list(fix=1, V=1e-15), G = list(G1 =
list(V=1, nu=0.002))),
pedigree = datatree,
reduced = TRUE,
burnin = 50000, nitt = 750001, thin = 700,
family = "threshold",
data = Rdata,
pl = TRUE, saveX = TRUE, saveZ = TRUE,
verbose = TRUE)
It returns some errors about missing functions "is.positive.definite" and
"Matrix", which I addressed with:
library("corpcor", lib.loc="~/R/win-library/3.5")
library("MatrixModels", lib.loc="~/R/win-library/3.5")
but I can't figure this one out:
'Error in .C("MCMCglmm", as.double(data$MCMC_y),
as.double(data$MCMC_y.additional), :
C symbol name "MCMCglmm" not in load table'
Detaching these packages doesn't necessarily cause that same error to
appear although I execute the exact same code.
Also, several attempts (same code again) caused a fatal error and
automatic session termination (info for a similar session below if
interesting).
I tried to use the fully bifurcating tree as an experiment but that made
no difference.
Any ideas what this last error means?
Thanks!
Roi
sessionInfo()
R version 3.5.1 (2018-07-02)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows >= 8 x64 (build 9200)
Matrix products: default
locale:
[1] LC_COLLATE=English_United Kingdom.1252 LC_CTYPE=English_United
Kingdom.1252 LC_MONETARY=English_United Kingdom.1252
[4] LC_NUMERIC=C LC_TIME=English_United
Kingdom.1252
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] corpcor_1.6.9 Matrix_1.2-14 phytools_0.6-60 maps_3.3.0
ape_5.2
loaded via a namespace (and not attached):
[1] igraph_1.2.2 Rcpp_0.12.19 magrittr_1.5
MASS_7.3-50 mnormt_1.5-5
[6] scatterplot3d_0.3-41 lattice_0.20-35 quadprog_1.5-5
fastmatch_1.1-0 tools_3.5.1
[11] parallel_3.5.1 grid_3.5.1 nlme_3.1-137
clusterGeneration_1.3.4 phangorn_2.4.0
[16] plotrix_3.7-4 coda_0.19-2 yaml_2.2.0
numDeriv_2016.8-1 animation_2.5
[21] compiler_3.5.1 combinat_0.0-8 expm_0.999-3
pkgconfig_2.0.2
On Tue, 20 Nov 2018 at 20:01, HADFIELD Jarrod <***@ed.ac.uk> wrote:
Hi,
Most likely the phylogenetic heritability (the phylogenetic variance /
the phylogenetic +residual variance) is approaching one resulting in
numerical difficulties. Probably the best thing is to assume that the
phylogenetic heritability equals 1 and use the reduced phylogenetic mixed
model implementation. This allows the phylogenetic heritability to be equal
to 1 without causing numerical issues. At some point I will integrate these
models into the main MCMCglmm package, but for now you can download it from
here: http://jarrod.bio.ed.ac.uk/MCMCglmmRAM_2.24.tar.gz.
Change the name of the ‘’Binomial’ column to ‘animal’ and fit:
T1 <- MCMCglmm(Activity ~ -1 + log(Mass) + Max.Temp * Annual.Precip,
random = ~ animal
prior = list(R = list(fix=1, V=1e-15), G =
list(G1 = list(V=1, nu=0.002))),
pedigree = tree,
reduced=TRUE,
burnin = 150000, nitt = 2650001, thin =
2500,
family = "threshold", data = Tdata,
pl = TRUE, pr = TRUE, saveX = TRUE, saveZ
= TRUE,
verbose = FALSE)
You should need fewer iterations.
Cheers,
Jarrod
On 20 Nov 2018, at 18:08, roee maor <***@gmail.com> wrote:
Dear Jarrod (and list),
Following your previous comment I added "random = ~ Binomial" to my model
to allow for a phylogenetic analysis.
This causes convergence problems: the trace plots show increasing
oscillations along each chain (although no directional trends, so it's not
a burn-in issue). Also, the posterior samples are highly correlated,
residual variance estimates are >10^3 and threshold estimates are high (>20
on the latent scale).
Surprisingly (to me), predictors that are strongly significant in the
non-phylogenetic model lose their effect in the phylogenetic model (I tried
several alternative parameter configurations).
It seems that this model attributes the explained variance to phylogeny
alone.
Can anyone explain what is going on here? Am I specifying the model
poorly or just asking my data more than it can answer?
I tried to overcome this issue by using a fully resolved variant of the
phylogeny, which only improved things slightly.
I also changed the random effect to "random=~Family" or "random=~Order",
which reduced the variance and threshold estimates to more acceptable
levels (<10), but still no significant predictors (and I'm not sure how the
algorithm calculates covariance between higher taxa in the phylogeny).
Separately I tried parameter expanded prior: "prior = list(R = list(V=1,
fix=1), G = list(G1 = list(V=1, nu=1, alpha.mu=0, alpha.V=1000)))". That
didn't help, and messing with priors for this reason feels like poor
practice.
This is the model:
T1 <- MCMCglmm(Activity ~ -1 + log(Mass) + Max.Temp * Annual.Precip,
random = ~ Binomial,
prior = list(R = list(fix=1, V=1), G =
list(G1 = list(V=1, nu=0.002))),
ginverse = list(Binomial=INphylo$Ainv),
burnin = 150000, nitt = 2650001, thin =
2500,
family = "threshold", data = Tdata,
pl = TRUE, pr = TRUE, saveX = TRUE, saveZ
= TRUE,
verbose = FALSE)
The data I use looks like this (not all variables appear in each model):
str(Tdata)
'data.frame': 1389 obs. of 10 variables:
$ Binomial : Factor w/ 1421 levels "Abrocoma_bennettii",..: 1
2 3 4 5 6 7 8 9 10 ...
$ Order : Factor w/ 27 levels "Afrosoricida",..: 24 24 24
24 24 3 24 24 2 2 ...
$ Family : Factor w/ 126 levels "Abrocomidae",..: 1 26 26
26 26 46 74 87 10 10 ...
$ Activity : Factor w/ 3 levels "1","2","3": 1 3 2 2 2 3 2
3 1 3 ...
$ Habitat : Factor w/ 6 levels "Aqua","Arbo",..: 5 5 5 5 5
5 5 5 5 5 ...
$ Diet : Factor w/ 3 levels "Faun","Herb",..: 2 3 3 3
3 1 3 2 2 2 ...
$ Mass : num 250.5 24.9 34.5 38.9 24.5 ...
$ Max.Temp : num 22 16.6 19.1 19.8 17.2 ...
$ Annual.Precip : num 166 645 558 903 1665 ...
Any advice would be much appreciated!
Many thanks,
--
Roi Maor
PhD candidate
School of Zoology, Tel Aviv University
Centre for Biodiversity and Environment Research, UCL
The University of Edinburgh is a charitable body, registered in
Scotland, with registration number SC005336.
--
Roi Maor
PhD candidate
School of Zoology, Tel Aviv University
Centre for Biodiversity and Environment Research, UCL
The University of Edinburgh is a charitable body, registered in
Scotland, with registration number SC005336.
--
Roi Maor
PhD candidate
School of Zoology, Tel Aviv University
Centre for Biodiversity and Environment Research, UCL
[[alternative HTML version deleted]]
------------------------------
Subject: Digest Footer
_______________________________________________
R-sig-mixed-models mailing list
R-sig-mixed-***@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
------------------------------
End of R-sig-mixed-models Digest, Vol 143, Issue 33
***************************************************
_____________________________________
Professor Diana Kornbrot
Mobile
+44 (0) 7403 18 16 12
Work
University of Hertfordshire
College Lane, Hatfield, Hertfordshire AL10 9AB, UK
+44 (0) 170 728 4626
***@herts.ac.uk<mailto:***@herts.ac.uk>
http://dianakornbrot.wordpress.com/
http://go.herts.ac.uk/Diana_Kornbrot
skype: kornbrotme
Home
19 Elmhurst Avenue
London N2 0LT, UK
+44 (0) 208 444 2081
------------------------------------------------------------
[[alternative HTML version deleted]]