[R-sig-ME] MCMCglmm with multinomial models

Jarrod Hadfield j@h@d||e|d @end|ng |rom ed@@c@uk
Wed Jul 27 06:14:30 CEST 2022


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). 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
>>>>
>>>
>>
>> --
>> 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

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>


	[[alternative HTML version deleted]]



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