[R-sig-ME] MCMCglmm with multinomial models
jessica comley
je@@|ecom|ey44 @end|ng |rom gm@||@com
Tue Jul 25 04:06:09 CEST 2023
Hi Jarrod,
Thank you for all your help on this last year. We submitted our paper and
have gotten comments back from a reviewer asking if we could test for the
interaction between culling activity and predator presence with the data we
have?
Below is the model I ran previously.
*prior1=list(R=list(V=1, nu=0.002))*
*m1<-MCMCglmm(cbind(dawn, diurnal, dusk, nocturnal)~trait
-1+trait:culling+trait:predator, rcov=~idv(units+trait:units),
prior=prior1, data=bbj, family="multinomial4", nitt= 150000)*
All the best,
Jess
On Wed, Jul 27, 2022 at 8:36 PM jessica comley <jessiecomley44 using gmail.com>
wrote:
> Okay great, thank you. I have run Wald tests and they show that culling is
> significant but predator presence is not.
>
> On Wed, Jul 27, 2022 at 4:35 PM Jarrod Hadfield <j.hadfield using ed.ac.uk>
> wrote:
>
>> Not quite - it means the odds of diurnal activity compared to nocturnal
>> activity is greater when there is no culling compared to lethal culling.
>> However, I would do the omnibus Wald test before interpreting the
>> significance of individual tests.
>>
>>
>> On 27 Jul 2022, at 06:31, jessica comley <jessiecomley44 using gmail.com>
>> wrote:
>>
>> This email was sent to you by someone outside the University.
>> You should only click on links or attachments if you are certain that the
>> email is genuine and the content is safe.
>> Sorry, let me try rephrasing my question. From my results, how do I
>> interpret the base activity nocturnal. Does the significant result of *traitdiurnal:cullingnone
>> *mean that there is more diurnal activity than nocturnal activity when
>> culling is none compared to lethal?
>>
>> On Wed, Jul 27, 2022 at 12:54 PM Jarrod Hadfield <j.hadfield using ed.ac.uk>
>> wrote:
>>
>>> The trait:culling interactions are the degree to which activity is
>>> higher when culling is none rather than lethal. If you exponentiate the
>>> effects you get the promotional change in the odds ratio none:lethal.
>>>
>>>
>>>
>>> On 27 Jul 2022, at 05:48, jessica comley <jessiecomley44 using gmail.com>
>>> wrote:
>>>
>>> This email was sent to you by someone outside the University.
>>> You should only click on links or attachments if you are certain that
>>> the email is genuine and the content is safe.
>>> Sorry Jarrod, one last question I swear.
>>>
>>> From my corrected outcome:
>>>
>>> *Location effects: cbind(dawn, diurnal, dusk, nocturnal) ~ trait - 1 +
>>> trait:culling + trait:predator*
>>>
>>> * post.mean l-95% CI u-95% CI eff.samp
>>> pMCMC *
>>> *traitdawn -1.49696 -1.91443 -1.03494 449.8 < 7e-05
>>> ****
>>> *traitdiurnal -1.26128 -1.65569 -0.87573 173.7 < 7e-05
>>> ****
>>> *traitdusk -1.69165 -2.14582 -1.25074 304.0 < 7e-05
>>> ****
>>> *traitdawn:cullingnone -0.35018 -0.75539 0.03488 1453.5 0.06068
>>> . *
>>> *traitdiurnal:cullingnone 0.61075 0.22682 0.99247 646.2 0.00653
>>> ***
>>> *traitdusk:cullingnone 0.32833 -0.07876 0.72035 1218.7
>>> 0.10136 *
>>> *traitdawn:predatorhigh 0.12826 -0.48062 0.81929 696.5
>>> 0.70884 *
>>> *traitdiurnal:predatorhigh -0.45382 -1.03340 0.15183 293.7
>>> 0.12313 *
>>> *traitdusk:predatorhigh -0.14870 -0.77156 0.49585 475.2
>>> 0.62830 *
>>> *traitdawn:predatorlow 0.41706 -0.11004 0.93341 579.5
>>> 0.10626 *
>>> *traitdiurnal:predatorlow -0.01279 -0.47660 0.49818 206.6
>>> 0.95741 *
>>> *traitdusk:predatorlow 0.16600 -0.36760 0.71419 323.0 0.55102
>>> *
>>>
>>> How would I work out whether traitnocturnal:cullingnone occured more or
>>> less than cullinglethal? Is this something that can be done?
>>>
>>> Cheers,
>>> Jess
>>>
>>> On Wed, Jul 27, 2022 at 12:24 PM jessica comley <
>>> jessiecomley44 using gmail.com> wrote:
>>>
>>>> Hi Jarrod,
>>>>
>>>> Yes I did have a typo in my data, which I have corrected and now all is
>>>> working well.
>>>>
>>>> Thank you very much for your help, I really do appreciate your
>>>> responses!
>>>>
>>>> Cheers,
>>>> Jess
>>>>
>>>> On Wed, Jul 27, 2022 at 12:14 PM Jarrod Hadfield <j.hadfield using ed.ac.uk>
>>>> wrote:
>>>>
>>>>> Hi,
>>>>>
>>>>> 1/ My guess is that there is a mistake/typo in your data.frame:
>>>>> culling has 3 levels not 2. Does table(bbj$culling) return what you expect?
>>>>>
>>>>> 2/ They mean i) there is more diurnal activity under culling compared
>>>>> to whatever the mystery level of culling is. ii) there is more dawn
>>>>> activity when predators is low compared to absent.
>>>>>
>>>>> 3/ The indices should be for all terms involving the thing to be
>>>>> tested. So in the current model they should be 4:9 and 10:15 (not 3:5 and
>>>>> 6:8). This will change when you sort out your culling column (probably to
>>>>> 4:6 and 7:12). In my previous email I said you should be testing 3 effects
>>>>> for predator, but in fact there should be 6 (I thought predator had 3
>>>>> levels not 2).
>>>>>
>>>>> You might want a -1 in your model formula (i.e
>>>>> trait-1+trait:culling+trait:predator) to make the interpretation of
>>>>> the first 3 terms a little easier, but up to you.
>>>>>
>>>>> Cheers,
>>>>>
>>>>> Jarrod
>>>>>
>>>>>
>>>>>
>>>>>
>>>>>
>>>>>
>>>>> On 27 Jul 2022, at 04:42, jessica comley <jessiecomley44 using gmail.com>
>>>>> wrote:
>>>>>
>>>>> This email was sent to you by someone outside the University.
>>>>> You should only click on links or attachments if you are certain that
>>>>> the email is genuine and the content is safe.
>>>>> Dear Jarrod,
>>>>>
>>>>> Sorry to bother you again, I just want to make sure I am doing this
>>>>> correctly and understanding my results.
>>>>>
>>>>> I used the model you suggested:
>>>>>
>>>>> *prior1=list(R=list(V=1, nu=0.002)) m1<-MCMCglmm(cbind(dawn, diurnal,
>>>>> dusk, nocturnal)~trait+trait:culling+trait:predator,
>>>>> rcov=~idv(units+trait:units), prior=prior1, data=bbj,
>>>>> family="multinomial4", nitt= 150000)*
>>>>>
>>>>> And this is my outcome:
>>>>>
>>>>>
>>>>>
>>>>>
>>>>>
>>>>>
>>>>>
>>>>>
>>>>>
>>>>>
>>>>>
>>>>>
>>>>>
>>>>>
>>>>>
>>>>>
>>>>>
>>>>>
>>>>>
>>>>>
>>>>>
>>>>>
>>>>>
>>>>>
>>>>>
>>>>>
>>>>>
>>>>>
>>>>>
>>>>>
>>>>> *Iterations = 3001:149991 Thinning interval = 10 Sample size =
>>>>> 14700 DIC: 10312.95 R-structure: ~idv(units + trait:units)
>>>>> post.mean l-95% CI u-95% CI eff.samp trait:units 0.01932 0.0002102
>>>>> 0.07042 441 Location effects: cbind(dawn, diurnal, dusk, nocturnal)
>>>>> ~ trait + trait:culling + trait:predator
>>>>> post.mean l-95% CI u-95% CI eff.samp pMCMC (Intercept)
>>>>> -2.018903 -2.677830 -1.335305 609.1 < 7e-05 *** traitdiurnal
>>>>> 0.542636 -0.363766 1.405230 644.6 0.22068 traitdusk
>>>>> -0.047952 -0.984923 0.917710 374.6 0.91850
>>>>> traitdawn:cullingLethal 0.534474 0.003076 1.058211 596.8 0.05524 .
>>>>> traitdiurnal:cullingLethal 0.232597 -0.369959 0.834347 404.7 0.42789
>>>>> traitdusk:cullingLethal 0.376191 -0.217707 0.964636 408.2
>>>>> 0.19782 traitdawn:cullingnone -0.163961 -0.573477 0.298182
>>>>> 2945.5 0.38245 traitdiurnal:cullingnone 0.674041 0.247724 1.101917
>>>>> 2825.0 0.00952 ** traitdusk:cullingnone 0.449710 -0.014263
>>>>> 0.874925 1683.9 0.05102 . traitdawn:predatorhigh 0.456119
>>>>> -0.206435 1.151283 561.7 0.18531 traitdiurnal:predatorhigh
>>>>> -0.303114 -0.976842 0.377294 479.9 0.36939 traitdusk:predatorhigh
>>>>> 0.108262 -0.674552 0.895819 256.7 0.76122
>>>>> traitdawn:predatorlow 0.756262 0.160407 1.323520 418.1 0.01279 *
>>>>> traitdiurnal:predatorlow 0.136875 -0.446984 0.750518 305.7 0.65619
>>>>> traitdusk:predatorlow 0.422497 -0.303586 1.145038 194.9
>>>>> 0.22857 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘
>>>>> ’ 1*
>>>>>
>>>>> 1) Why do the 2 categories for culling both show up but then only 2 of
>>>>> the three categories for predator show up? i.e. predatorabsent is missing?
>>>>>
>>>>> 2) Do these results mean that i) diurnal activity and lethal culling
>>>>> is sig different from nocturnal activity and lethal culling; ii) dawn
>>>>> activity and predator low is sig different from nocturnal activity and
>>>>> predator low?
>>>>>
>>>>> 3) Is this the correct way and interpretation of the within group
>>>>> effects?
>>>>> *##culling lethal*
>>>>>
>>>>>
>>>>> * aod::wald.test(cov(m1$Sol[,3:5]),
>>>>> colMeans(m1$Sol[,3:5]),Terms=1:3)$result$chi2["P"] P 0.1638938*
>>>>>
>>>>>
>>>>>
>>>>>
>>>>> * ##culling none aod::wald.test(cov(m1$Sol[,6:8]),
>>>>> colMeans(m1$Sol[,6:8]),Terms=1:3)$result$chi2["P"] P 0.006497424*
>>>>>
>>>>>
>>>>> So these results show us that culling none has an effect on activity?
>>>>>
>>>>> Thank you in advance,
>>>>> Jess
>>>>>
>>>>> On Wed, Jul 27, 2022 at 8:20 AM jessica comley <
>>>>> jessiecomley44 using gmail.com> wrote:
>>>>>
>>>>>> Dear Jarrod,
>>>>>>
>>>>>> Thank you so much for your help, I greatly appreciate it!
>>>>>>
>>>>>> All the best,
>>>>>> Jess
>>>>>>
>>>>>> On Wed, Jul 27, 2022 at 3:32 AM Jarrod Hadfield <j.hadfield using ed.ac.uk>
>>>>>> wrote:
>>>>>>
>>>>>>> Hi Jess,
>>>>>>>
>>>>>>> Section should definitely not be left out, but I would imagine it is
>>>>>>> going to be very difficult to separate culling, predator and Section
>>>>>>> effects - I would expect the credible intervals to be large.
>>>>>>>
>>>>>>> As mentioned in my previous post you can test for an effect of
>>>>>>> culling by fitting the model
>>>>>>>
>>>>>>> ~trait+trait:culling+trait:predator
>>>>>>>
>>>>>>> And then fitting a Wald test to the three terms with 'culling' in.
>>>>>>> The effect of predator can be tested similarly but with the 3 terms with
>>>>>>> 'predator' in.
>>>>>>>
>>>>>>> Since your covariates do not vary within Section it will be much
>>>>>>> easier to aggregate the counts at the Section level (i.e have a data frame
>>>>>>> with 14 rows and 1 column for each activity with the number observed for
>>>>>>> each activity) and fit family="multinomial". You can then get rid of the
>>>>>>> random formula as the Section effects are now effectively the residuals.
>>>>>>> Given the lack of replication I would advise using the idv formula that I
>>>>>>> suggested previously and hope the model isn't too misspecified:
>>>>>>>
>>>>>>> prior=list(R=list(V=1, nu=0.002))
>>>>>>>
>>>>>>> m1<-MCMCglmm(cbind(dawn, diurnal, dusk,
>>>>>>> nocturnal)~trait+trait:culling+trait:predator,
>>>>>>> rcov=~idv(units+trait:units), prior=prior, ...)
>>>>>>>
>>>>>>> Note this models is identical to the original model, it's just
>>>>>>> parameterised in a more efficient way.
>>>>>>>
>>>>>>> Cheers,
>>>>>>>
>>>>>>> Jarrod
>>>>>>>
>>>>>>>
>>>>>>>
>>>>>>>
>>>>>>> On 25 Jul 2022, at 03:52, jessica comley <jessiecomley44 using gmail.com>
>>>>>>> wrote:
>>>>>>>
>>>>>>> This email was sent to you by someone outside the University.
>>>>>>> You should only click on links or attachments if you are certain
>>>>>>> that the email is genuine and the content is safe.
>>>>>>> Dear Jarrod and Walid,
>>>>>>>
>>>>>>> Thank you for your replies, it is greatly appreciated.
>>>>>>>
>>>>>>> The predator and culling factors do not vary within sites. As shown
>>>>>>> in the example data in one of my previous emails, Bucklands only has
>>>>>>> culling as lethal and predator as low, whereas Colchester only has predator
>>>>>>> as high and culling as none.
>>>>>>>
>>>>>>> We are trying to submit a paper on black-backed jackal and caracal
>>>>>>> activity in the presence of different culling practices and
>>>>>>> predator presence. The reviewers want us to try a GLMM approach to
>>>>>>> determine whether culling or predators have an effect on black-backed
>>>>>>> jackal or caracal activity.
>>>>>>>
>>>>>>> Therefore, in your opinion how could be go about this given our
>>>>>>> data? Would it be advisable to leave out the random effect of Section?
>>>>>>>
>>>>>>> All the best,
>>>>>>> Jess
>>>>>>>
>>>>>>> On Wed, Jul 20, 2022 at 3:06 PM Jarrod Hadfield <j.hadfield using ed.ac.uk>
>>>>>>> wrote:
>>>>>>>
>>>>>>>> Hi Jess
>>>>>>>>
>>>>>>>> In multinomial models the linear model is set up as a (logit)
>>>>>>>> difference in probability between an outcome and some base-line outcome.
>>>>>>>> Often, as here, the base-line outcome is arbitrary, and so the idh
>>>>>>>> structure is a little odd. For example, if A is the base line category, idh
>>>>>>>> assumes COV(B-A, C-A) = 0 which therefore assumes
>>>>>>>> COV(B,C)+VAR(A) =COV(A,B)+COV(C,A). It's not clear why this would
>>>>>>>> be the case. Perhaps a more reasonable, but less parameter rich, option
>>>>>>>> would be to have:
>>>>>>>>
>>>>>>>> ~idv(Section+trait:Section)
>>>>>>>>
>>>>>>>> which parameterises the Section covariance matrix by a single
>>>>>>>> parameter (rather than 6). The term idv(Section+trait:Section) fits a 3x3
>>>>>>>> covariance matrix of the form v*(I+J) where v is the estimated variance.
>>>>>>>> This assumes i) Sections are repeatable in outcome, but knowing that a
>>>>>>>> Section has an increased 'preference' for A doesn’t tell you whether it
>>>>>>>> also has an increased preference for one of the other categories and ii)
>>>>>>>> the repeatability for each outcome within sites is the same (on the latent
>>>>>>>> scale).
>>>>>>>>
>>>>>>>> To test groups of effects (in your case the 3 culling:trait
>>>>>>>> effects), I usually use a Wald test and the posterior covariances (see here
>>>>>>>>
>>>>>>>> https://stat.ethz.ch/pipermail/r-sig-mixed-models/2017q3/025930.html).
>>>>>>>> It's far from correct and so Walid's suggestions may be better, but
>>>>>>>> small-scale simulations suggests it has good frequentist properties.
>>>>>>>>
>>>>>>>> To add predator presence you can just add a predator:trait effect
>>>>>>>> into the linear model. If the culling and predator factors do not vary
>>>>>>>> within sites then you probably don't have enough information to reliably
>>>>>>>> estimate these effects.
>>>>>>>>
>>>>>>>> Cheers,
>>>>>>>>
>>>>>>>> Jarrod
>>>>>>>>
>>>>>>>>
>>>>>>>>
>>>>>>>>
>>>>>>>>
>>>>>>>>
>>>>>>>> > On 19 Jul 2022, at 18:17, Walid Mawass <walidmawass10 using gmail.com>
>>>>>>>> wrote:
>>>>>>>> >
>>>>>>>> > This email was sent to you by someone outside the University.
>>>>>>>> > You should only click on links or attachments if you are certain
>>>>>>>> that the email is genuine and the content is safe.
>>>>>>>> >
>>>>>>>> > Hey Jess,
>>>>>>>> >
>>>>>>>> > 1) Yes that is correct
>>>>>>>> >
>>>>>>>> > 2) To my knowledge there is a rule of thumb, where you set the
>>>>>>>> nitt (# of
>>>>>>>> > iterations) to a large number that includes the burnin amount,
>>>>>>>> then you
>>>>>>>> > choose your thinning interval (sampling of the chain). For
>>>>>>>> example, this is
>>>>>>>> > what I would use: nitt= 150000, burnin=50000, thin=100. This will
>>>>>>>> give you
>>>>>>>> > a decent burnin and a final sample of 1000 saved iterations. Note
>>>>>>>> however
>>>>>>>> > that this does not have to increase the effective sample size for
>>>>>>>> certain
>>>>>>>> > variables, but it might do the trick.
>>>>>>>> >
>>>>>>>> > 3) hmm...I think one way to do it is to make predictions using
>>>>>>>> the above
>>>>>>>> > model and interpret the patterns you see for each relationship
>>>>>>>> you are
>>>>>>>> > interested in. Another way to compare effect size would be to use
>>>>>>>> bayesian
>>>>>>>> > posterior indices. I suggest these two papers by Makowski et al.
>>>>>>>> (2019a &
>>>>>>>> > b) that present both interesting posterior indices to use with
>>>>>>>> Bayesian
>>>>>>>> > statistical analysis and an associated R package that does the
>>>>>>>> job of
>>>>>>>> > computing these indices, *bayestestR*.
>>>>>>>> >
>>>>>>>> > Good luck
>>>>>>>> > --
>>>>>>>> > Walid Mawass
>>>>>>>> > Ph.D. candidate in Evolutionary Biology - UQTR
>>>>>>>> > *Currently* Postdoctoral Research Associate
>>>>>>>> > Masel Lab - University of Arizona
>>>>>>>> >
>>>>>>>> >
>>>>>>>> > On Sun, Jul 17, 2022 at 11:32 PM jessica comley <
>>>>>>>> jessiecomley44 using gmail.com>
>>>>>>>> > wrote:
>>>>>>>> >
>>>>>>>> >> Hi Walid,
>>>>>>>> >>
>>>>>>>> >> Thank you for your reply, I greatly appreciate it. I have a few
>>>>>>>> more
>>>>>>>> >> questions and if you could help that would be great.
>>>>>>>> >>
>>>>>>>> >> I tested for correlation between activities and the 14 Sections
>>>>>>>> and the
>>>>>>>> >> correlation comes out as low. Therefore I have changed my code
>>>>>>>> to use idh()
>>>>>>>> >> instead of us as suggested:
>>>>>>>> >>
>>>>>>>> >> test1c.5b <- MCMCglmm(activity ~ -1 + at.level(culling,1):trait +
>>>>>>>> >> at.level(culling, 2):trait, random=~idh(trait):Section, rcov =
>>>>>>>> >> ~idh(trait):units, data = caracal, family = "categorical", prior
>>>>>>>> = prior,
>>>>>>>> >> burnin=5000, nitt=80000)
>>>>>>>> >>
>>>>>>>> >> 1) Is this correct?
>>>>>>>> >>
>>>>>>>> >> 2) Increasing the number of interactions increases the effective
>>>>>>>> sample
>>>>>>>> >> size, therefore is there a general rule of thumb as to how large
>>>>>>>> your
>>>>>>>> >> effective sample size should be?
>>>>>>>> >>
>>>>>>>> >> 3) I understand how to use and interpret the results of
>>>>>>>> HPDinterval (i.e.
>>>>>>>> >> if intervals do not overlap 0 then relationship is strong), but
>>>>>>>> how am I
>>>>>>>> >> able to test the relationship between all four activities and
>>>>>>>> fixed effects
>>>>>>>> >> and not just have the three categories (i.e. diurnal, dusk,
>>>>>>>> nocturnal)
>>>>>>>> >> compared to the base category (dawn)? For example, I am also
>>>>>>>> interested in
>>>>>>>> >> whether there is a significant/strong relationship between
>>>>>>>> activities of
>>>>>>>> >> caracal at dusk with culling(Lethal)/no culling(none) compared to
>>>>>>>> >> activities of caracal at diurnal with culling(Lethal)/no
>>>>>>>> culling(none).
>>>>>>>> >>
>>>>>>>> >> Below is an example of our dataset:
>>>>>>>> >> Camera Section CameraID Animal predator culling activity
>>>>>>>> >> 1a Bucklands Bucklands1a Caracal low Lethal diurnal
>>>>>>>> >> 1a Bucklands Bucklands1a Caracal low Lethal dawn
>>>>>>>> >> 2a Bucklands Bucklands2a Caracal low Lethal dusk
>>>>>>>> >> 2a Bucklands Bucklands2a Caracal low Lethal diurnal
>>>>>>>> >> 3a Bucklands Bucklands3a Caracal low Lethal dawn
>>>>>>>> >> Cam 1 Colchester ColchesterCam 1 Caracal high none diurnal
>>>>>>>> >> Cam 1 Colchester ColchesterCam 1 Caracal high none diurnal
>>>>>>>> >> Cam 1 Colchester ColchesterCam 1 Caracal high none diurnal
>>>>>>>> >> Cam 1 Colchester ColchesterCam 1 Caracal high none diurnal
>>>>>>>> >> Cam 2 Colchester ColchesterCam 2 Caracal high none diurnal
>>>>>>>> >> Cam 2 Colchester ColchesterCam 2 Caracal high none diurnal
>>>>>>>> >> Cam 3 Colchester ColchesterCam 3 Caracal high none diurnal
>>>>>>>> >> Cam 3 Colchester ColchesterCam 3 Caracal high none diurnal
>>>>>>>> >> Cam 3 Colchester ColchesterCam 3 Caracal high none diurnal
>>>>>>>> >> Cam 4 Colchester ColchesterCam 4 Caracal high none diurnal
>>>>>>>> >> Cam 4 Colchester ColchesterCam 4 Caracal high none diurnal
>>>>>>>> >> Cam 4 Colchester ColchesterCam 4 Caracal high none nocturnal
>>>>>>>> >> 1a Connaught Connaught1a Caracal low Lethal nocturnal
>>>>>>>> >> 1a Connaught Connaught1a Caracal low Lethal nocturnal
>>>>>>>> >> 1d Connaught Connaught1d Caracal low Lethal diurnal
>>>>>>>> >> 3B Connaught Connaught3B Caracal low Lethal diurnal
>>>>>>>> >> 3B Connaught Connaught3B Caracal low Lethal diurnal
>>>>>>>> >> 4a Connaught Connaught4a Caracal low Lethal nocturnal
>>>>>>>> >> 4a Connaught Connaught4a Caracal low Lethal nocturnal
>>>>>>>> >> 4b Connaught Connaught4b Caracal low Lethal diurnal
>>>>>>>> >> 6a Connaught Connaught6a Caracal low Lethal nocturnal
>>>>>>>> >> 6b Connaught Connaught6b Caracal low Lethal diurnal
>>>>>>>> >> 7a Connaught Connaught7a Caracal low Lethal nocturnal
>>>>>>>> >> 9a Connaught Connaught9a Caracal low Lethal nocturnal
>>>>>>>> >> 9d Connaught Connaught9d Caracal low Lethal nocturnal
>>>>>>>> >> 9d Connaught Connaught9d Caracal low Lethal dusk
>>>>>>>> >> 7d Diepdam Diepdam7d Caracal absent Lethal dusk
>>>>>>>> >> 8d Diepdam Diepdam8d Caracal absent Lethal diurnal
>>>>>>>> >> 9c Diepdam Diepdam9c Caracal absent Lethal nocturnal
>>>>>>>> >>
>>>>>>>> >> All the best,
>>>>>>>> >> Jess
>>>>>>>> >>
>>>>>>>> >>
>>>>>>>> >> On Fri, Jul 15, 2022 at 11:37 PM Walid Mawass <
>>>>>>>> walidmawass10 using gmail.com>
>>>>>>>> >> wrote:
>>>>>>>> >>
>>>>>>>> >>> Hello,
>>>>>>>> >>>
>>>>>>>> >>> I don't think I can specifically help you with some of your
>>>>>>>> inquiries.
>>>>>>>> >>> However, I do want to comment on a few things that might need
>>>>>>>> some
>>>>>>>> >>> attention.
>>>>>>>> >>>
>>>>>>>> >>> First, MCMCglmm is based on a Bayesian implementation and does
>>>>>>>> not
>>>>>>>> >>> compute p-values to compare. What you need to compare are the
>>>>>>>> posterior
>>>>>>>> >>> distributions of your effect sizes. This can be done visually
>>>>>>>> using the
>>>>>>>> >>> base plot function in R. Or by comparing the HPD intervals and
>>>>>>>> the mode (or
>>>>>>>> >>> mean) of the posterior distributions.
>>>>>>>> >>>
>>>>>>>> >>> Second, I have no idea what your data structure looks like
>>>>>>>> (which makes
>>>>>>>> >>> it hard to interpret model results), but the effective sample
>>>>>>>> size (from
>>>>>>>> >>> the 5500 saved iterations sample) for your random variable
>>>>>>>> Section is very
>>>>>>>> >>> low (the same applies for your fixed effects). You should
>>>>>>>> consider this
>>>>>>>> >>> issue and look again at your assumption of correlation between
>>>>>>>> >>> activities for the 14 sections you have in your dataset. If you
>>>>>>>> do not
>>>>>>>> >>> expect among activity correlations then you can use the idh()
>>>>>>>> function
>>>>>>>> >>> instead of us().
>>>>>>>> >>>
>>>>>>>> >>> Hopefully this helps and in hope that people on this list with
>>>>>>>> more
>>>>>>>> >>> knowledge of these models will help out.
>>>>>>>> >>>
>>>>>>>> >>> Best,
>>>>>>>> >>> --
>>>>>>>> >>> Walid Mawass
>>>>>>>> >>> Ph.D. candidate in Evolutionary Biology - UQTR
>>>>>>>> >>> *Currently* Postdoctoral Research Associate
>>>>>>>> >>> Masel Lab - University of Arizona
>>>>>>>> >>>
>>>>>>>> >>>
>>>>>>>> >>> On Fri, Jul 15, 2022 at 8:49 AM jessica comley <
>>>>>>>> jessiecomley44 using gmail.com>
>>>>>>>> >>> wrote:
>>>>>>>> >>>
>>>>>>>> >>>> Dear all,
>>>>>>>> >>>>
>>>>>>>> >>>> I am hoping that someone will be able to help me with
>>>>>>>> conducting MCMCglmm
>>>>>>>> >>>> multinomial models.
>>>>>>>> >>>>
>>>>>>>> >>>> The data I am working with is for black-backed jackal (bbj)
>>>>>>>> and carcal.
>>>>>>>> >>>> For
>>>>>>>> >>>> each species we have a multinomial response variable called
>>>>>>>> activity
>>>>>>>> >>>> which
>>>>>>>> >>>> has four categories (dawn, diurnal, dusk, nocturnal). We have
>>>>>>>> two
>>>>>>>> >>>> categorical fixed effects which are 1) culling (none, lethal)
>>>>>>>> and 2)
>>>>>>>> >>>> predator presence (absent, high, low). We also have a
>>>>>>>> categorical
>>>>>>>> >>>> variable
>>>>>>>> >>>> called Section (made up of 14 different reserves/ farms where
>>>>>>>> the
>>>>>>>> >>>> activity
>>>>>>>> >>>> of caracal and bbj were recorded). There are 273 observations
>>>>>>>> for caracal
>>>>>>>> >>>> and 4399 for bbj. We are wanting to test the effects of
>>>>>>>> culling and
>>>>>>>> >>>> predators on caracal and bbj activity separately.
>>>>>>>> >>>>
>>>>>>>> >>>> I have been working through Jarrod Hadfields course notes,
>>>>>>>> particularly
>>>>>>>> >>>> with regards to Chapter 5.2. The chi-square analyses reveal
>>>>>>>> that the
>>>>>>>> >>>> frequencies of culling and predators differ as do activities.
>>>>>>>> >>>>
>>>>>>>> >>>> I have managed to work out the specific probabilities for the
>>>>>>>> culling
>>>>>>>> >>>> none
>>>>>>>> >>>> vs culling lethal for each activity (dawn, diurnal, dusk,
>>>>>>>> nocturnal) for
>>>>>>>> >>>> caracal, but I'm confused as to how to determine p-values to
>>>>>>>> determine
>>>>>>>> >>>> which activities culling none vs culling lethal are affecting?
>>>>>>>> >>>>
>>>>>>>> >>>> Myy code and outcomes are pasted below with questions stated
>>>>>>>> in bold.
>>>>>>>> >>>>
>>>>>>>> >>>> caracal2 <- read.csv("caracal_new.csv", header=T)
>>>>>>>> >>>> caracal <- as.data.frame(unclass(caracal2), stringsAsFactors =
>>>>>>>> TRUE)
>>>>>>>> >>>>
>>>>>>>> >>>> #Chi-squared tests
>>>>>>>> >>>> Ctable1 <- table(caracal$activity, caracal$culling)
>>>>>>>> >>>> chisq.test(rowSums(Ctable1)) #strongly suggests activities
>>>>>>>> differ
>>>>>>>> >>>> chisq.test(Ctable1)#strongly suggests culling category differs
>>>>>>>> >>>>
>>>>>>>> >>>> Ctable2 <- table(caracal$activity, caracal$predator)
>>>>>>>> >>>> chisq.test(rowSums(Ctable2))#strongly suggests activities
>>>>>>>> differ
>>>>>>>> >>>> chisq.test(Ctable2)#strongly suggests predator category differs
>>>>>>>> >>>>
>>>>>>>> >>>> prior = list(R = list(fix=1, V=(1/k) * (I + J)), G =
>>>>>>>> list(G1=list(V =
>>>>>>>> >>>> diag(k-1), nu=1)))
>>>>>>>> >>>> test1c.5 <- MCMCglmm(activity ~ -1 + at.level(culling,1):trait
>>>>>>>> +
>>>>>>>> >>>> at.level(culling, 2):trait, random=~us(trait):Section, rcov =
>>>>>>>> >>>> ~us(trait):units, data = caracal, family = "categorical",
>>>>>>>> prior = prior,
>>>>>>>> >>>> burnin=5000, nitt=60000)
>>>>>>>> >>>> *##I'm not sure how to add the three predator levels to this
>>>>>>>> model or if
>>>>>>>> >>>> it
>>>>>>>> >>>> would be appropriate?*
>>>>>>>> >>>>
>>>>>>>> >>>>
>>>>>>>> >>>> k <- length(levels(caracal$activity))
>>>>>>>> >>>> I <- diag(k-1)
>>>>>>>> >>>> J <- matrix(rep(1, (k-1)^2), c(k-1, k-1))
>>>>>>>> >>>> IJ <- (1/k) *(diag(k-1) + matrix(1,k-1, k-1))
>>>>>>>> >>>>
>>>>>>>> >>>> contrasts(caracal$activity)
>>>>>>>> >>>>
>>>>>>>> >>>> #culling lethal
>>>>>>>> >>>> Delta <- cbind(c(0,1,0,0), c(0,0,1,0), c(0,0,0,1))
>>>>>>>> >>>> c2 <- (16 * sqrt(3)/(15 * pi))^2
>>>>>>>> >>>> D <- ginv(Delta %*% t(Delta)) %*% Delta
>>>>>>>> >>>> Int <- t(apply(test1c.5$Sol[,1:3],1, function(x) + D %*%
>>>>>>>> (x/sqrt(1 + c2 *
>>>>>>>> >>>> diag(IJ)))))
>>>>>>>> >>>> summary(mcmc(exp(Int)/rowSums(exp(Int))))
>>>>>>>> >>>>
>>>>>>>> >>>> prop.table(Ctable1[,1])
>>>>>>>> >>>>
>>>>>>>> >>>> #culling none
>>>>>>>> >>>> Delta <- cbind(c(0,1,0,0), c(0,0,1,0), c(0,0,0,1))
>>>>>>>> >>>> c2 <- (16 * sqrt(3)/(15 * pi))^2
>>>>>>>> >>>> D <- ginv(Delta %*% t(Delta)) %*% Delta
>>>>>>>> >>>> Int <- t(apply(test1c.5$Sol[,4:6],1, function(x) + D %*%
>>>>>>>> (x/sqrt(1 + c2 *
>>>>>>>> >>>> diag(IJ)))))
>>>>>>>> >>>> summary(mcmc(exp(Int)/rowSums(exp(Int))))
>>>>>>>> >>>>
>>>>>>>> >>>> prop.table((Ctable1[,2]))
>>>>>>>> >>>>
>>>>>>>> >>>> HPDinterval(test1c.5$Sol)
>>>>>>>> >>>>
>>>>>>>> >>>> #model summary
>>>>>>>> >>>>> summary(test1c.5)
>>>>>>>> >>>>
>>>>>>>> >>>> Iterations = 5001:59991
>>>>>>>> >>>> Thinning interval = 10
>>>>>>>> >>>> Sample size = 5500
>>>>>>>> >>>>
>>>>>>>> >>>> DIC: 699.7014
>>>>>>>> >>>>
>>>>>>>> >>>> G-structure: ~us(trait):Section
>>>>>>>> >>>>
>>>>>>>> >>>>
>>>>>>>> post.mean l-95%
>>>>>>>> >>>> CI
>>>>>>>> >>>> u-95% CI eff.samp
>>>>>>>> >>>> traitactivity.diurnal:traitactivity.diurnal.Section
>>>>>>>> 1.8124
>>>>>>>> >>>> 0.09784
>>>>>>>> >>>> 5.665 77.01
>>>>>>>> >>>> traitactivity.dusk:traitactivity.diurnal.Section
>>>>>>>> 0.8450
>>>>>>>> >>>> -0.83585
>>>>>>>> >>>> 3.856 64.17
>>>>>>>> >>>> traitactivity.nocturnal:traitactivity.diurnal.Section
>>>>>>>> 1.3621
>>>>>>>> >>>> -1.19129
>>>>>>>> >>>> 6.157 58.48
>>>>>>>> >>>> traitactivity.diurnal:traitactivity.dusk.Section
>>>>>>>> 0.8450
>>>>>>>> >>>> -0.83585
>>>>>>>> >>>> 3.856 64.17
>>>>>>>> >>>> traitactivity.dusk:traitactivity.dusk.Section
>>>>>>>> 1.2034
>>>>>>>> >>>> 0.07090
>>>>>>>> >>>> 3.681 102.16
>>>>>>>> >>>> traitactivity.nocturnal:traitactivity.dusk.Section
>>>>>>>> 0.7505
>>>>>>>> >>>> -1.77113
>>>>>>>> >>>> 4.524 43.53
>>>>>>>> >>>> traitactivity.diurnal:traitactivity.nocturnal.Section
>>>>>>>> 1.3621
>>>>>>>> >>>> -1.19129
>>>>>>>> >>>> 6.157 58.48
>>>>>>>> >>>> traitactivity.dusk:traitactivity.nocturnal.Section
>>>>>>>> 0.7505
>>>>>>>> >>>> -1.77113
>>>>>>>> >>>> 4.524 43.53
>>>>>>>> >>>> traitactivity.nocturnal:traitactivity.nocturnal.Section
>>>>>>>> 2.7148
>>>>>>>> >>>> 0.09401
>>>>>>>> >>>> 8.397 76.59
>>>>>>>> >>>>
>>>>>>>> >>>> R-structure: ~us(trait):units
>>>>>>>> >>>>
>>>>>>>> >>>> post.mean
>>>>>>>> l-95% CI
>>>>>>>> >>>> u-95% CI eff.samp
>>>>>>>> >>>> traitactivity.diurnal:traitactivity.diurnal.units
>>>>>>>> 0.50 0.50
>>>>>>>> >>>> 0.50 0
>>>>>>>> >>>> traitactivity.dusk:traitactivity.diurnal.units
>>>>>>>> 0.25 0.25
>>>>>>>> >>>> 0.25 0
>>>>>>>> >>>> traitactivity.nocturnal:traitactivity.diurnal.units
>>>>>>>> 0.25 0.25
>>>>>>>> >>>> 0.25 0
>>>>>>>> >>>> traitactivity.diurnal:traitactivity.dusk.units
>>>>>>>> 0.25 0.25
>>>>>>>> >>>> 0.25 0
>>>>>>>> >>>> traitactivity.dusk:traitactivity.dusk.units
>>>>>>>> 0.50 0.50
>>>>>>>> >>>> 0.50 0
>>>>>>>> >>>> traitactivity.nocturnal:traitactivity.dusk.units
>>>>>>>> 0.25 0.25
>>>>>>>> >>>> 0.25 0
>>>>>>>> >>>> traitactivity.diurnal:traitactivity.nocturnal.units
>>>>>>>> 0.25 0.25
>>>>>>>> >>>> 0.25 0
>>>>>>>> >>>> traitactivity.dusk:traitactivity.nocturnal.units
>>>>>>>> 0.25 0.25
>>>>>>>> >>>> 0.25 0
>>>>>>>> >>>> traitactivity.nocturnal:traitactivity.nocturnal.units
>>>>>>>> 0.50 0.50
>>>>>>>> >>>> 0.50 0
>>>>>>>> >>>>
>>>>>>>> >>>> Location effects: activity ~ -1 + at.level(culling, 1):trait +
>>>>>>>> >>>> at.level(culling, 2):trait
>>>>>>>> >>>>
>>>>>>>> >>>> post.mean l-95% CI
>>>>>>>> u-95% CI
>>>>>>>> >>>> eff.samp pMCMC
>>>>>>>> >>>> at.level(culling, 1):traitactivity.diurnal 1.2306
>>>>>>>> -0.0533 2.6793
>>>>>>>> >>>> 145.29 0.0418 *
>>>>>>>> >>>> at.level(culling, 1):traitactivity.dusk 0.6605
>>>>>>>> -0.6006 2.0761
>>>>>>>> >>>> 92.91 0.2840
>>>>>>>> >>>> at.level(culling, 1):traitactivity.nocturnal 1.6090
>>>>>>>> 0.0914 3.1356
>>>>>>>> >>>> 151.02 0.0265 *
>>>>>>>> >>>> traitactivity.diurnal:at.level(culling, 2) 1.2664
>>>>>>>> -0.1552 2.7750
>>>>>>>> >>>> 226.40 0.0604 .
>>>>>>>> >>>> traitactivity.dusk:at.level(culling, 2) 0.3533
>>>>>>>> -0.9898 1.5218
>>>>>>>> >>>> 148.44 0.5447
>>>>>>>> >>>> traitactivity.nocturnal:at.level(culling, 2) 1.0447
>>>>>>>> -0.6405 2.8354
>>>>>>>> >>>> 346.40 0.1618
>>>>>>>> >>>> ---
>>>>>>>> >>>> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
>>>>>>>> >>>>
>>>>>>>> >>>> *##So for the model summary I get that lethal culling at
>>>>>>>> activity diurnal
>>>>>>>> >>>> is significantly different from lethal culling at dawn (its
>>>>>>>> the base
>>>>>>>> >>>> reference), but I'm also interested in whether lethal culling
>>>>>>>> at activity
>>>>>>>> >>>> diurnal is different from lethal culling at dusk for example.
>>>>>>>> Is this
>>>>>>>> >>>> possible? *
>>>>>>>> >>>>
>>>>>>>> >>>> #outcomes culling lethal
>>>>>>>> >>>>> summary(mcmc(exp(Int)/rowSums(exp(Int))))
>>>>>>>> >>>>
>>>>>>>> >>>> Iterations = 1:5500
>>>>>>>> >>>> Thinning interval = 1
>>>>>>>> >>>> Number of chains = 1
>>>>>>>> >>>> Sample size per chain = 5500
>>>>>>>> >>>>
>>>>>>>> >>>> 1. Empirical mean and standard deviation for each variable,
>>>>>>>> >>>> plus standard error of the mean:
>>>>>>>> >>>>
>>>>>>>> >>>> Mean SD Naive SE Time-series SE
>>>>>>>> >>>> [1,] 0.1253 0.05565 0.0007504 0.002484
>>>>>>>> >>>> [2,] 0.3748 0.10497 0.0014155 0.003204
>>>>>>>> >>>> [3,] 0.1757 0.06640 0.0008954 0.002515
>>>>>>>> >>>> [4,] 0.3242 0.11939 0.0016099 0.003514
>>>>>>>> >>>>
>>>>>>>> >>>> 2. Quantiles for each variable:
>>>>>>>> >>>>
>>>>>>>> >>>> 2.5% 25% 50% 75% 97.5%
>>>>>>>> >>>> var1 0.03641 0.08695 0.1198 0.1554 0.2553
>>>>>>>> >>>> var2 0.17298 0.30580 0.3704 0.4431 0.5896
>>>>>>>> >>>> var3 0.06166 0.12913 0.1705 0.2161 0.3215
>>>>>>>> >>>> var4 0.12610 0.23999 0.3090 0.3901 0.6045
>>>>>>>> >>>>
>>>>>>>> >>>>> prop.table(Ctable1[,1])
>>>>>>>> >>>> dawn diurnal dusk nocturnal
>>>>>>>> >>>> 0.1250000 0.2812500 0.1770833 0.4166667
>>>>>>>> >>>>
>>>>>>>> >>>>
>>>>>>>> >>>> #outcomes culling none
>>>>>>>> >>>>> summary(mcmc(exp(Int)/rowSums(exp(Int))))
>>>>>>>> >>>>
>>>>>>>> >>>> Iterations = 1:5500
>>>>>>>> >>>> Thinning interval = 1
>>>>>>>> >>>> Number of chains = 1
>>>>>>>> >>>> Sample size per chain = 5500
>>>>>>>> >>>>
>>>>>>>> >>>> 1. Empirical mean and standard deviation for each variable,
>>>>>>>> >>>> plus standard error of the mean:
>>>>>>>> >>>>
>>>>>>>> >>>> Mean SD Naive SE Time-series SE
>>>>>>>> >>>> [1,] 0.1288 0.06141 0.0008280 0.002787
>>>>>>>> >>>> [2,] 0.3804 0.10406 0.0014032 0.002662
>>>>>>>> >>>> [3,] 0.1710 0.06844 0.0009228 0.002592
>>>>>>>> >>>> [4,] 0.3198 0.11812 0.0015928 0.002956
>>>>>>>> >>>>
>>>>>>>> >>>> 2. Quantiles for each variable:
>>>>>>>> >>>>
>>>>>>>> >>>> 2.5% 25% 50% 75% 97.5%
>>>>>>>> >>>> var1 0.02891 0.08896 0.1220 0.1594 0.2685
>>>>>>>> >>>> var2 0.18007 0.31094 0.3783 0.4474 0.5965
>>>>>>>> >>>> var3 0.05840 0.12425 0.1634 0.2083 0.3250
>>>>>>>> >>>> var4 0.12430 0.23921 0.3077 0.3862 0.5964
>>>>>>>> >>>>
>>>>>>>> >>>>> prop.table((Ctable1[,2]))
>>>>>>>> >>>> dawn diurnal dusk nocturnal
>>>>>>>> >>>> 0.1306818 0.4375000 0.1875000 0.2443182
>>>>>>>> >>>>
>>>>>>>> >>>> Any help or guidance will be greatly appreciated.
>>>>>>>> >>>>
>>>>>>>> >>>> All the best,
>>>>>>>> >>>> Jess
>>>>>>>> >>>>
>>>>>>>> >>>> --
>>>>>>>> >>>> Jessica Comley (PhD)
>>>>>>>> >>>> Research Scientist
>>>>>>>> >>>>
>>>>>>>> >>>> [[alternative HTML version deleted]]
>>>>>>>> >>>>
>>>>>>>> >>>> _______________________________________________
>>>>>>>> >>>> R-sig-mixed-models using r-project.org mailing list
>>>>>>>> >>>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>>>>>>>> >>>>
>>>>>>>> >>>
>>>>>>>> >>
>>>>>>>> >> --
>>>>>>>> >> Jessica Comley (PhD)
>>>>>>>> >> Research Scientist
>>>>>>>> >>
>>>>>>>> >>
>>>>>>>> >
>>>>>>>> > [[alternative HTML version deleted]]
>>>>>>>> >
>>>>>>>> > _______________________________________________
>>>>>>>> > R-sig-mixed-models using r-project.org mailing list
>>>>>>>> > https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>>>>>>>>
>>>>>>>> The University of Edinburgh is a charitable body, registered in
>>>>>>>> Scotland, with registration number SC005336. Is e buidheann carthannais a
>>>>>>>> th’ ann an Oilthigh Dhùn Èideann, clàraichte an Alba, àireamh clàraidh
>>>>>>>> SC005336.
>>>>>>>>
>>>>>>>
>>>>>>>
>>>>>>> --
>>>>>>> Jessica Comley (PhD)
>>>>>>> Research Scientist
>>>>>>>
>>>>>>>
>>>>>>>
>>>>>>
>>>>>> --
>>>>>> Dr Jessica Comley
>>>>>> Lecturer: Environmental and Life Sciences
>>>>>> Faculty of Science
>>>>>> Universiti Brunei Darussalam
>>>>>>
>>>>>> Email: jessica.comley using ubd.edu.bn
>>>>>>
>>>>>>
>>>>>
>>>>> --
>>>>> Dr Jessica Comley
>>>>> Lecturer: Environmental and Life Sciences
>>>>> Faculty of Science
>>>>> Universiti Brunei Darussalam
>>>>>
>>>>> Email: jessica.comley using ubd.edu.bn
>>>>>
>>>>>
>>>>>
>>>>
>>>> --
>>>> Dr Jessica Comley
>>>> Lecturer: Environmental and Life Sciences
>>>> Faculty of Science
>>>> Universiti Brunei Darussalam
>>>>
>>>> Email: jessica.comley using ubd.edu.bn
>>>>
>>>>
>>>
>>> --
>>> Dr Jessica Comley
>>> Lecturer: Environmental and Life Sciences
>>> Faculty of Science
>>> Universiti Brunei Darussalam
>>>
>>> Email: jessica.comley using ubd.edu.bn
>>>
>>>
>>>
>>
>> --
>> Dr Jessica Comley
>> Lecturer: Environmental and Life Sciences
>> Faculty of Science
>> Universiti Brunei Darussalam
>>
>> Email: jessica.comley using ubd.edu.bn
>>
>>
>>
>
> --
> Dr Jessica Comley
> Lecturer: Environmental and Life Sciences
> Faculty of Science
> Universiti Brunei Darussalam
>
> Email: jessica.comley using ubd.edu.bn
>
>
--
Dr Jessica Comley
Lecturer: Environmental and Life Sciences
Faculty of Science
Universiti Brunei Darussalam
Email: jessica.comley using ubd.edu.bn
[[alternative HTML version deleted]]
More information about the R-sig-mixed-models
mailing list