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