[R-sig-ME] MCMCglmm with multinomial models

jessica comley je@@|ecom|ey44 @end|ng |rom gm@||@com
Wed Jul 27 06:24:43 CEST 2022


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

	[[alternative HTML version deleted]]



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