Dear Malcom,
Yes, they are covariates that I didn't control as they are a result of my
visual surveys as well. I thought I couldn't consider them as fixed because
of that.
But I follow your suggestions and include them as fixed effects as I had in
the first place:
glm1 <- glm (Biomass ~ fReserve * Roughness + fReserve * DivBoulders ,
family= Gamma, data= myData)
Many many thanks
Barbara
On 18 September 2012 11:14, Malcolm Fairbrother wrote:
> Hi Barbara,
>
> I'm not really understanding what your random classification is, and thus
> why you need a mixed/multilevel model. From the sounds of it, yes, all
> fixed effects should work. What you're currently treating as random
> classifications--DivBoulders and Roughness--just sound like typical
> covariates/predictors, in which case a single-level regression (using R's
> "lm" or "glm") will do.
>
> Good luck,
> Malcolm
>
>
>
> On 17 Sep 2012, at 11:15, barbara costa wrote:
>
> Dear Malcolm,
> thanks a lot for your answer.
>
> My random variables are indices from habitat features I surveyed inside
> and outside one marine reserve that I would like to account for when
> assessing protection effects. Thus, I thought they should be random effects
> as I did not chose them and they are only a sample from variables with a
> wider range of values.
> On the other hand, this way I'm having a doubt on how to interpret my
> results. When I was using all effects as fixed and assuming I was only
> interested for this range of roughness and boulders diversity, I could
> easily test the interaction among those preditors and the protection
> (fReserve).
> a) Do you think I could handle this problem assuming all as fixed effects?
> b) If not (which I think is the most correct option), can you help me
> interpreting these results?
>
> lmer(Biomass ~ fReserve + (1 | DivBoulders) + (1 | Roughness), myData)
>
> Linear mixed model fit by REML
> Formula: Biomass ~ fReserve + (1 | DivBoulders) + (1 | Roughness)
> Data: myData
> AIC BIC logLik deviance REMLdev
> -2973 -2957 1492 -2981 -2983
> Random effects:
> Groups Name Variance Std.Dev.
> Roughness (Intercept) 1.8348e+00 1.3546e+00
> DivBoulders (Intercept) 3.2292e-11 5.6826e-06
> Residual 6.4554e-09 8.0346e-05
> Number of obs: 200, groups: Roughness, 10; DivBoulders, 6
>
> Fixed effects:
> Estimate Std. Error t value
> (Intercept) 10.4599 0.6058 17.267
> fReserveOut -3.4399 0.8567 -4.015
>
> Correlation of Fixed Effects:
> (Intr)
> fReserveOut -0.707
>
>
> i) How can I test the influence or interaction of those random variables
> in the obtained reserve effect?
> ii) And should I test the significance of t-value by myself? In lme I had
> the p-values with summary (lme object).
>
> thanks a lot again.
> Barbara
>
>
>
>
> On 15 September 2012 11:57, Malcolm Fairbrother <
> m.fairbrother@bristol.ac.uk> wrote:
>
>> Hi Barbara,
>>
>> The lme4 package may be a bit easier for you to use, and is more current
>> than nlme in various ways (though nlme does do a few things lme4 can't do).
>> I think the problem is the right hand side of the random part of your
>> model: I don't believe you can have a "+" there, as you have in your first
>> two models. Your third model definitely doesn't make sense, because you're
>> treating fReserve as both a covariate and a random classification.
>>
>> It's a bit hard to understand your query generally. Variables treated as
>> fixed effects can be either categorical or continuous--not a problem. But I
>> don't understand the random classification in your case, which is the part
>> for which you probably want to calculate a variance term (and why you need
>> a mixed model). A random classification is usually coded as 1 to 20, or A
>> to Z, or something like that, whereas your Roughness and DivBoulders look
>> like continuous variables.
>>
>> However, if you are trying to model observations cross-classified in both
>> Roughness and DivBoulders, and treating fReserve as the only
>> covariate/predictor, then try lme4's lmer, along the lines of:
>>
>> lmer(Biomass ~ fReserve + (1 | DivBoulders) + (1 | Roughness), myData)
>>
>> Hope that helps. Follow up with further clarification if not.
>>
>> Cheers,
>> Malcolm
>>
>>
>>
>> > Date: Fri, 14 Sep 2012 18:34:28 +0100
>> > From: barbara costa
>> > To: r-sig-mixed-models@r-project.org
>> > Subject: [R-sig-ME] mixed models with two random variables?
>> >
>> > Dear R users,
>> > Does anyone knows how to run a glmm with one fixed factor and 2 random
>> > numeric variables (indices)? Is there any way to force in the model a
>> > separate interaction of those random variables with the fixed one?
>> > I hope you can help me.
>> >
>> > #eg.
>> > Reserve <- rep(c("In","Out"), 100)
>> > fReserve <- factor(Reserve)
>> > DivBoulders <- rep (c(1.23,2.4,1.26,1.78,1.97,1.35,1.23,2.4,1.26,1.78),
>> 20)
>> > Roughness <-
>> rep(c(3.45,2.56,1.32,5.67,3.73,3.57,2.66,1.52,7.67,2.73),20)
>> > Biomass <- rep(c(8,5.3,3.5,12,25.4,10.1,9.8,2.4,5.6,5.3),20)
>> >
>> > myData <- data.frame (fReserve ,DivBoulders ,Roughness ,Biomass )
>> >
>> > #glm
>> > glm1 <- glm (Biomass ~ fReserve * Roughness + fReserve * DivBoulders ,
>> > family= Gamma, data= myData)
>> >
>> > #glmm:
>> > library (nlme)
>> >
>> > lme1 <- lme (Biomass ~ fReserve , random= ~1 + fReserve | Roughness
>> > + DivBoulders , data= myData ) # random intercept and slope - I
>> suspect
>> > it's my case
>> > #Error in getGroups.data.frame(dataMix, groups) :
>> > # Invalid formula for groups
>> > # if I only use one random variable I have:
>> > #Error in chol.default((value + t(value))/2) :
>> > #the leading minor of order 2 is not positive definite
>> >
>> >
>> > lme2 <- lme( Biomass ~ fReserve , random = ~1 | Roughness
>> > + DivBoulders ,data=myData) #random intercept
>> > #Error in getGroups.data.frame(dataMix, groups) :
>> > # Invalid formula for groups
>> > # if I only use one random variable my result is fine!
>> >
>> >
>> > lme3 <- lme (Biomass ~ fReserve , random= ~ Roughness +
>> DivBoulders
>> > | fReserve , data= myData ) # from help (lme)
>> > summary (lme3)
>> > # I have a result. Is the model correct for what I want?
>> > # BUT my fixed effect (reserve) does not have a p-value due to zero
>> degrees
>> > of freedom. However in glm1 or lm2 with one random variable it has.
>> >
>> > Can you help me please?
>> >
>> > Thanks a lot in advance,
>> > Barbara
>>
>>
>
>
[[alternative HTML version deleted]]