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