[R-sig-ME] parameter estimates for all factor levels MCMCglmm

[Conor Michael Goold]

Hi Dena, 

I personally find removing the intercept without pre-defining how the categorical predictors are coded is confusing. 

It might be easiest if you set up the contrasts you want to look at before running the model, then you can keep track of exactly what the output is showing. A great webpage for different types of coding for categorical predictors is http://www.ats.ucla.edu/stat/r/library/contrast_coding.htm. 

If you are running a large model (e.g. many predictors/levels of categorical variables), or a model that you may repeat with different data, I would set up the design matrix of predictors first using the model.matrix() function, within which you can specify a list of contrast arguments (see contrast.arg in the model.matrix() function). Then you can plug that into the regression formula. In the long run, I've found it saves a lot of time doing it this way. As a simple example:

#==========================================================================

require(nlme)
data("Machines")
d <- Machines

# set predictors as factor variables
d[,c("Machine","Worker")] <- apply(d[,c("Machine","Worker")], 2, as.factor)

# define model matrix, coding predictors as deflections from their reference categories
designMatrix <- model.matrix(~ Machine + Worker, data = d, 
                             contrasts.arg = list(Machine = "contr.treatment", Worker = "contr.treatment"))

# glm model, using Worker as a fixed effect rather than a random effect
# remove intercept because it is already included in designMatrix. Alternatively, remove the intercept column from the design matrix 
summary(fit <- glm(d$score ~ -1 + designMatrix))

#==========================================================================

You can then post-process the fitted model object as you like. 

Hope that helps a little.

Best regards
Conor Goold
PhD Student
Phone:        +47 67 23 27 24



Norwegian University of Life Sciences
Campus Ås. www.nmbu.no

________________________________________
From: R-sig-mixed-models <r-sig-mixed-models-bounces at r-project.org> on behalf of Dena Paris <dp at dparis.com>
Sent: Tuesday, December 6, 2016 3:00 AM
To: Jarrod Hadfield
Cc: r-sig-mixed-models at r-project.org
Subject: Re: [R-sig-ME] parameter estimates for all factor levels MCMCglmm

Hi Jarrod,

Thanks for the very prompt reply! I’m still struggling with this a bit. My understanding was that by removing the global intercept (-1) all factor levels would be estimated, and if I have only one fixed effect, say Arthropod, this is exactly what happens. However, when I include the second fixed effect, Size, it returns estimates for all 16 levels of Arthropod, but estimates for only 2 of the 3 Size classes. So, in the case with two fixed effects and the global intercept removed, is there a way to calculate the estimate for the third Size class, or is it the level from which the other levels deviate?

Thanks again for all your help.

Dena

> On 5 Dec. 2016, at 5:49 pm, Jarrod Hadfield <j.hadfield at ed.ac.uk> wrote:
>
> Hi,
>
> You have an intercept which will be the estimate for the first level of Arthropod and Size, the remaining effects are deviations from these levels (15 and 2 contrasts respectively).
>
> Cheers,
>
> Jarrod
>
>
>
> On 05/12/2016 04:29, Dena Paris wrote:
>> Hi all,
>>
>> I'm new to Bayesian stats and the MCMCglmm package. I'm trying to understand how to interpret the results from a fitted model. I've fitted a model with a binomial response (reject/select), two categorical fixed effects (taxon with 16 levels, and size class with 3 levels), and a single random effect (bird ID). My model is mixing well and the results fit the data. My problem is that I'm not getting estimates for all levels of the fixed effects (19). How do I get the post mean and CIs for all levels in order to correctly interpret/write up the results?
>>
>> As I have (near) complete separation in the data, I've used the fixed effect prior structure suggested in the Course Notes, fixed the residuals, and removed the global intercept.
>>
>> prior.1 = list(
>>   B = list(mu = rep(0, 18), V = (diag(18)) * (1 + pi^2/3)),
>>   R = list(fix=1, V=1, n = k - 1),
>>   G = list(G1 = list(V = 1, n = 1))
>> )?
>>
>> m.1 <- MCMCglmm(Selected ~ -1 + Arthropod + Size, random = ~bID, family = "categorical", prior = prior.1, data = type.selected, verbose = FALSE, nitt = 5e+05, burnin = 5000, thin = 100)
>>
>> Thank you for your guidance,
>> Dena
>>
>>
>> Dena Paris
>>
>> -----
>> Dena Paris
>> Ph.D Candidate
>> School of Environmental Sciences
>> Institute for Land, Water and Society
>> Charles Sturt University
>> PO Box 789
>> Albury NSW 2640
>> M: +61 424 451 858?
>>
>>
>>      [[alternative HTML version deleted]]
>>
>> _______________________________________________
>> R-sig-mixed-models at 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.
>

_______________________________________________
R-sig-mixed-models at r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models