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