[R-sig-ME] MCMCglmm with multinomial models

Jarrod Hadfield j@h@d||e|d @end|ng |rom ed@@c@uk
Wed Jul 27 10:35:19 CEST 2022


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 <mailto: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 <mailto: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 <mailto: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 <mailto: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 <mailto: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 <mailto: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 <mailto: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 <mailto: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 <mailto: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 <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 <mailto: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 <mailto: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 <mailto: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 <mailto: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 <mailto:R-sig-mixed-models using r-project.org> mailing list
>>>> >>>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models <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 <mailto:R-sig-mixed-models using r-project.org> mailing list
>>>> > https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models <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 <mailto: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 <mailto: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 <mailto: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 <mailto: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 <mailto:jessica.comley using ubd.edu.bn>
> 


	[[alternative HTML version deleted]]



More information about the R-sig-mixed-models mailing list