[R-sig-ME] How to code proportional data with gamm binomial model using gamm4 package?

[Simon Wood]

Julien,

Apologies for the cbind problem: I'll see what I can do about that.

The different results problem is probably because neither of your model 
specifications is quite right. For the `cbind' form of the response you 
should provide something ot the form cbind(success,failure), not 
cbind(success, success+failure). For the proportions form then you must 
supply a `weights' argument containing your `total'.

best,
Simon

On Friday 19 March 2010 19:30, Julien Beguin wrote:
> Dear user,
>  
> I am trying to fit an additive mixed model with proportional data as
> response variable using gamm4 package. My proportional data are composed of
> 1) the number of time an individual was observed at location i (=presence),
> and 2) the number of time each location i was visited (=total). In my data
> set, I also have one random variable (=bloc) and a set of continuous
> independent variables that could explain part of the response variable
> variability. 
> Unfortunately, I have trouble to specify the response variable in the
> model. When I try to code proportional dataas "cbind(presence,total)", I
> don't have access to the "gam" object in the summary() function, while I
> can get the "mer" object which is, fortunately, similar to that from
> glmer() model... When I try to code proportional data directly as
> "proportion" ([0,1]) where proportion = presence/total, I can get all what
> I want except that model parameters and AIC values does not match the
> previous one anymore. I guess that I am doing something wrong here... 
> My question is how to code proportional response variables to get the "gam"
> object in the summary() function and that gives consistant results of AIC,
> parameter estimates, etc...? 
> My R code is below. It is a simple version where I included only one
> independent variable and no smooth term yet but in my real data set, I will
> use several independent variables, potentially with smoothers. 
> Any help would be really appreciated,
>  
> Many thanks,
>  
> Julien Beguin
>  Laval University
>  Ph.D. candidate
>
>
>  ----------------------------------------------- R code
> ---------------------------------------- #################
>  # simulated data#
>  #################
>  
> set.seed(1001)
>  presence <- round(runif(50,0,15))
>  total <- round(runif(50,50,70))
>  prop <- presence/total
>  bloc <- rep(c(1,2,3,4,5), each=10)
>  x1 <- rnorm(50, 20, 5)
>  mydata <- as.data.frame(cbind(presence, total, prop, bloc, x1))
>  mydata$bloc <- as.factor(mydata$bloc)
>  
> #############################################################
>  # gamm model with response variable as cbind(presence,total)#
>  #############################################################
>
> mod1 <- gamm4(cbind(presence,total)~ x1,
> random=~(1|bloc),family=binomial(link = "logit"), data=mydata)
> summary(mod1$gam)
>
> >Error in object$y - object$fitted.values : non-conformable arrays
>
>  
> summary(mod1$mer)
>
> > AIC   BIC logLik deviance
> >
>  > 150.7 156.4 -72.33    144.7
>  >
>  >Fixed effects:
>  >             Estimate Std. Error z value Pr(>|z|)   
>  >X(Intercept) -1.54827    0.21700  -7.135  9.7e-13 ***
>  >Xx1          -0.02927    0.01110  -2.637  0.00837 **
>
>  
> mod2 <- gamm4(prop~ x1, random=~(1|bloc),family=binomial(link = "logit"),
> data=mydata)
>
>  summary(mod2$gam)
>
> >Parametric coefficients:
>  >            Estimate Std. Error z value Pr(>|z|)
>  >(Intercept) -1.35978    1.68946  -0.805    0.421
>  >x1          -0.03236    0.08600  -0.376    0.707
>
>  
> summary(mod2$mer)
>
> > AIC   BIC logLik deviance
> >
>  > 8.983 14.72 -1.492    2.983
> >
> >Fixed effects:
>  >             Estimate Std. Error z value Pr(>|z|)
>  >X(Intercept) -1.35978    1.68946  -0.805    0.421
>  >Xx1          -0.03236    0.08600  -0.376    0.707
>
>  
> mod3 <- glmer(cbind(presence,total)~ x1 + (1|bloc), family=binomial(link =
> "logit"), data=mydata, REML=FALSE) summary(mod3)
>
> > AIC   BIC logLik deviance
> >
>  > 150.7 156.4 -72.33    144.7
> >
> >Fixed effects:
>  >            Estimate Std. Error z value Pr(>|z|)   
>  >(Intercept) -1.54827    0.21700  -7.135  9.7e-13 ***
>  >x1          -0.02927    0.01110  -2.637  0.00837 **
>
>  
> mod4 <- glmer(prop~ x1 + (1|bloc), family=binomial(link = "logit"),
> data=mydata, REML=FALSE) summary(mod4)
>
> > AIC   BIC logLik deviance
> >
>  > 8.983 14.72 -1.492    2.983
> >
> >Fixed effects:
>  >            Estimate Std. Error z value Pr(>|z|)
>  >(Intercept) -1.35978    1.68946  -0.805    0.421
>  >x1          -0.03236    0.08600  -0.376    0.707

-- 
> Simon Wood, Mathematical Sciences, University of Bath, Bath, BA2 7AY UK
> +44 1225 386603  www.maths.bath.ac.uk/~sw283