5
$\begingroup$

Simplified version of my model:

glm(cbind(young, adults) ~ as.factor(month) + effort, family = "binomial")

i.e., I study proportion of young as a dependent variable on month (or season), taking into account observer effort. However, the observer effort is dependent on the month:

enter image description here

How to solve this problem? I want to take into account both variables.

I was looking into literature but haven't found any solution. My naive idea is to compute mean effort for each month and instead of taking effort, take difference between effort and this mean. But this is just a naive idea, I would like to hear your advice. Thanks!

EDIT - response to Scortchi question - no, not actually:

> m = glm(cbind(young, adults) ~ as.factor(month) + effort, family = "quasibinomial")
> summary(m)

Call:
glm(formula = cbind(young, adults) ~ as.factor(month) + effort, family = "quasibinomial")

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-2.6829  -1.1138   0.0000   0.9717   4.0090  

Coefficients:
                     Estimate Std. Error t value Pr(>|t|)    
(Intercept)        -0.7868764  0.2170738  -3.625 0.000306 ***
as.factor(month)2   0.8857561  0.2606780   3.398 0.000710 ***
as.factor(month)3   0.7055741  0.2918895   2.417 0.015843 *  
as.factor(month)4   0.3943665  0.3269973   1.206 0.228138    
as.factor(month)5   0.4831113  0.3730987   1.295 0.195713    
as.factor(month)6  -0.5217349  0.5027560  -1.038 0.299676    
as.factor(month)7   0.1612901  0.4333682   0.372 0.709851    
as.factor(month)8   0.5114890  0.3545159   1.443 0.149444    
as.factor(month)9   0.7741060  0.3126087   2.476 0.013466 *  
as.factor(month)10  0.6601093  0.2609937   2.529 0.011608 *  
as.factor(month)11  0.4891778  0.2647303   1.848 0.064967 .  
as.factor(month)12  0.4743091  0.2565709   1.849 0.064849 .  
effort              0.0032506  0.0007976   4.075 5.02e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

(Dispersion parameter for quasibinomial family taken to be 1.270962)

    Null deviance: 1518.0  on 878  degrees of freedom
Residual deviance: 1447.8  on 866  degrees of freedom
  (750 observations deleted due to missingness)
AIC: NA

Number of Fisher Scoring iterations: 4

Warning message:
In summary.glm(m) :
  observations with zero weight not used for calculating dispersion
> 
> 
> m = glm(cbind(young, adults) ~ as.factor(month), family = "quasibinomial")
> summary(m)

Call:
glm(formula = cbind(young, adults) ~ as.factor(month), family = "quasibinomial")

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-3.1142  -1.1266   0.0000   0.9235   3.6484  

Coefficients:
                   Estimate Std. Error t value Pr(>|t|)    
(Intercept)         -0.2296     0.1590  -1.444  0.14890    
as.factor(month)2    0.8610     0.2059   4.181 3.09e-05 ***
as.factor(month)3    1.0887     0.2062   5.280 1.51e-07 ***
as.factor(month)4    0.4184     0.2374   1.762  0.07822 .  
as.factor(month)5    0.1495     0.3086   0.485  0.62802    
as.factor(month)6   -0.6177     0.3872  -1.595  0.11091    
as.factor(month)7   -0.4636     0.3666  -1.265  0.20622    
as.factor(month)8    0.1089     0.2976   0.366  0.71440    
as.factor(month)9    0.4932     0.2490   1.980  0.04787 *  
as.factor(month)10   0.6322     0.2096   3.016  0.00261 ** 
as.factor(month)11   0.4919     0.2152   2.286  0.02243 *  
as.factor(month)12   0.2296     0.2127   1.079  0.28071    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

(Dispersion parameter for quasibinomial family taken to be 1.309215)

    Null deviance: 2400.1  on 1345  degrees of freedom
Residual deviance: 2310.4  on 1334  degrees of freedom
AIC: NA

Number of Fisher Scoring iterations: 4

Warning message:
In summary.glm(m) :
  observations with zero weight not used for calculating dispersion
> 
> 
> 
> m = glm(cbind(young, adults) ~ effort, family = "quasibinomial")
> summary(m)

Call:
glm(formula = cbind(young, adults) ~ effort, family = "quasibinomial")

Deviance Residuals: 
   Min      1Q  Median      3Q     Max  
-2.718  -1.119   0.000   1.011   4.236  

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) -0.326899   0.102473  -3.190  0.00147 ** 
effort       0.003827   0.000688   5.563 3.52e-08 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

(Dispersion parameter for quasibinomial family taken to be 1.268467)

    Null deviance: 1518.0  on 878  degrees of freedom
Residual deviance: 1475.4  on 877  degrees of freedom
  (750 observations deleted due to missingness)
AIC: NA

Number of Fisher Scoring iterations: 4

Warning message:
In summary.glm(m) :
  observations with zero weight not used for calculating dispersion
$\endgroup$
  • $\begingroup$ Are you getting very high variances for coefficient estimates when both month & effort are in the model - compared to what you see when you put each in separately? If not, why do you think there's a problem? $\endgroup$ – Scortchi Jan 28 '13 at 13:59
  • $\begingroup$ I would naively suggest trying package lme4's lmer, where something like lmer (cbind (young, adults) ~ (effort | month), family="binomial") might be what you want. Scortchi's comment sounds much more informed than I can be, though. $\endgroup$ – Wayne Jan 28 '13 at 14:08
  • 2
    $\begingroup$ Both might be important predictors. There's no requirement for predictors to be perfectly orthogonal in generalized linear models, any more than there is in linear regression or ANOVA. If the predictors are too collinear the model fit will still be valid, but useless because the variances for individual model coefficients become enormous. When that doesn't happen you have no cause to worry. $\endgroup$ – Scortchi Jan 28 '13 at 14:33
  • 1
    $\begingroup$ @Scortchi, I'm afraid you are wrong in this. The collinear variables will compete for the fraction of variability they explain. And I don't care just about how the whole model fits, but about the importance of each predictor variable. $\endgroup$ – Curious Jan 28 '13 at 15:29
  • 1
    $\begingroup$ @Tomas: Well worry if you like, but you'll have some degree of collinearity in almost any observational data-set, so you'll be doing a lot of worrying. To my eyes this data-set doesn't look so bad, but if you want a quantitative assessment of the degree of collinearity calculate variance inflation factors & condition indices: stats.stackexchange.com/questions/16692/… stats.stackexchange.com/questions/4099/… $\endgroup$ – Scortchi Jan 28 '13 at 17:05
2
$\begingroup$

I'm converting a previous comment to an answer, expanding a bit based on a follow-up comment from the OP. The original, unedited comment was:

There is no silver bullet for decomposing variation in that situation. One thing you can do with two collinear predictors, $x_1,x_2$, is fit a model $x_1 \sim x_2$, take the residuals from that model, $η$, and replace $x_1$ with $η$ in the model $y \sim x_1+x_2$. This way, you will, definitionally, have uncorrelated predictors and the contribution of η is thought of as the variance explained by $x_1$ that is not subsumed by $x_2$. Of course, which variable is $x_1$ and which is $x_2$ is a judgment call (though the overall model fit will be identical).

In response to the OP's comment:

@Macro, this is a nice thing... maybe worth posting an answer, so we can discuss it with more detail? This is very interesting, because then $x_1=x_2+η$, and if you replace the x1 with η in the original model, you get $y \sim η+x_2=x_1$, which means you loose $x_2$ for the overall fit of the model! And this is strange, paradox! Please post your comment as an answer to discuss it in more detail.

Be careful here, because $x_1 \sim x_2$ is R pseudo-code for the model $x_1 = \beta_0 + \beta_1 x_2 + \eta$, not $x_1 = x_2 + \eta$. So, by my back-of-the-envelope calculation, this means that the model $y \sim \eta + x_2$, which is short hand for $y = \alpha_0 + \alpha_1 \eta + \alpha_2 x_2 + \varepsilon$, can be written as

$$ y = (\alpha_0 - \alpha_1 \beta_0) + \alpha_1 x_1 + (\alpha_2 - \alpha_1 \beta_1) x_2 + \varepsilon $$

So $x_2$ does not drop out of the model. Indeed the model $y \sim \eta + x_2$ can be seen to have identical degrees of freedom, fit statistics, etc. to the model $y \sim x_1 + x_2$, but the predictors are now uncorrelated.

$\endgroup$
  • 1
    $\begingroup$ Thank you Macro! BTW, this means that if I compute the model $y \sim x_1 + x_2$, I can simply transform its coefficients to those from the model $y \sim \eta + x_2$ without even computing this model (but I need coefficients from $x_1 \sim x_2$). Interesting! $\endgroup$ – Curious Jan 28 '13 at 18:07
  • $\begingroup$ You're very welcome, @Tomas! $\endgroup$ – Macro Jan 28 '13 at 23:39

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.