[R-sig-ME] MCMCglmm with multinomial models
jessica comley
je@@|ecom|ey44 @end|ng |rom gm@||@com
Wed Jul 27 14:36:48 CEST 2022
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
[[alternative HTML version deleted]]
More information about the R-sig-mixed-models
mailing list