Tweedie GLM, what var.power should I chose?

Question

Im trying to do a GLM on my data about bats to see how different variables affect bat activity on 8 species. Orignally my data was count data, but because of hardware difficulties I had to divide the count data with amount of recorded hours so I get comparable data. Data also has many zero-values and together with positive non-intergers and four different explaining variables, I thought a Tweedie GLM with gamma distribution was the right choice.

Here you have my dataset:

    > str(flagermusdata)
tibble [91 x 46] (S3: tbl_df/tbl/data.frame)
 $ Location                 : chr [1:91] "Stenderup" "Stenderup" "Stenderup" "Stenderup" ...
 $ AM-nr                    : num [1:91] 35 36 41 43 40 36 35 52 31 32 ...
 $ Recorded hour            : num [1:91] 14 2.5 10.5 11 3.5 3.5 14.5 14.5 14.5 14.5 ...
 $ lat                      : num [1:91] 55.5 55.5 55.5 55.5 55.4 ...
 $ lon                      : num [1:91] 9.65 9.67 9.62 9.62 9.44 ...
 $ StandAge                 : num [1:91] 103 217 86 27 70 65 46 13 44 27 ...
 $ StandSp                  : chr [1:91] "Beech" "Beech" "Beech" "Other" ...
 $ DistWater                : num [1:91] 600 10 0 100 200 0 100 100 250 200 ...
 $ DistOpen                 : num [1:91] 10 10 20 10 5 30 15 0 10 5 ...
 $ Region                   : chr [1:91] "Sønderjylland/fyn" "Sønderjylland/fyn" "Sønderjylland/fyn" "Sønderjylland/fyn" ...
 $ Plecotus auritus         : num [1:91] 0 0 0 0 0 6 0 1 0 0 ...
 $ Barbastella barbastellus : num [1:91] 0 0 0 0 0 0 0 0 0 0 ...
 $ Nyctalus Notula          : num [1:91] 0 0 0 0 0 6 1 0 0 0 ...
 $ Nyctalus Leisleri        : num [1:91] 0 0 0 0 0 0 0 0 0 0 ...
 $ Eptesicus serotinus      : num [1:91] 0 0 0 0 0 0 0 4 7 7 ...
 $ Eptesicus nilssoni       : num [1:91] 0 0 0 0 0 0 0 0 0 0 ...
 $ Myotis nattereri         : num [1:91] 0 0 0 0 0 0 0 1 0 0 ...
 $ Myotis brandti           : num [1:91] 0 0 0 0 0 0 0 0 0 0 ...
 $ Myotis mystacinus        : num [1:91] 0 0 0 0 0 0 0 0 0 0 ...
 $ Myotis dasycneme         : num [1:91] 0 0 0 0 0 0 0 0 0 0 ...
 $ Myotis daubentoni        : num [1:91] 0 0 0 0 0 33 1 21 0 0 ...
 $ Myotis bechsteini        : num [1:91] 0 0 0 0 0 0 0 0 0 0 ...
 $ Myotis myotis            : num [1:91] 0 0 0 0 0 0 0 0 0 0 ...
 $ Vespertilio murinus      : num [1:91] 0 0 0 0 0 0 0 0 2 0 ...
 $ Pipistrellus nathusii    : num [1:91] 0 0 0 0 1 49 0 9 1 1 ...
 $ Pipistrellus pipistrellus: num [1:91] 0 0 0 1 69 19 5 63 0 4 ...
 $ Pipistrellus pygmaeus    : num [1:91] 0 0 0 1 0 14 0 58 0 0 ...
 $ NoSpec                   : num [1:91] 0 0 0 2 2 6 3 6 3 3 ...
 $ PLEAURactivity           : num [1:91] 0 0 0 0 0 ...
 $ BARBARactivity           : num [1:91] 0 0 0 0 0 0 0 0 0 0 ...
 $ NYCNOCactivity           : num [1:91] 0 0 0 0 0 ...
 $ NYCLEIactivity           : num [1:91] 0 0 0 0 0 0 0 0 0 0 ...
 $ EPTSERactivity           : num [1:91] 0 0 0 0 0 ...
 $ EPTNILactivity           : num [1:91] 0 0 0 0 0 0 0 0 0 0 ...
 $ MYONATactivity           : num [1:91] 0 0 0 0 0 ...
 $ MYOBRAactivity           : num [1:91] 0 0 0 0 0 0 0 0 0 0 ...
 $ MYOMYSactivity           : num [1:91] 0 0 0 0 0 0 0 0 0 0 ...
 $ MYODASactivity           : num [1:91] 0 0 0 0 0 0 0 0 0 0 ...
 $ MYODAUactivity           : num [1:91] 0 0 0 0 0 ...
 $ MYOBECactivity           : num [1:91] 0 0 0 0 0 0 0 0 0 0 ...
 $ MYOMYOactivity           : num [1:91] 0 0 0 0 0 0 0 0 0 0 ...
 $ VESMURactivity           : num [1:91] 0 0 0 0 0 ...
 $ PIPNATactivity           : num [1:91] 0 0 0 0 0.286 ...
 $ PIPPIPactivity           : num [1:91] 0 0 0 0.0909 19.7143 ...
 $ PIPPYGactivity           : num [1:91] 0 0 0 0.0909 0 ...
 $ Totalactivity            : num [1:91] 0 0 0 0.182 20 ...

With some of the species, the tests went fine:

Call:
glm(formula = PIPNATactivity ~ StandAge + DistWater + DistOpen + 
    StandSp, family = tweedie(var.power = 2, link.power = 0), 
    data = flagermusdata, maxit = 100)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-2.8715  -2.1226  -1.1956   0.1082   2.9974  

Coefficients:
                    Estimate Std. Error t value Pr(>|t|)    
(Intercept)        -0.698292   0.651280  -1.072   0.2873    
StandAge            0.007474   0.004581   1.631   0.1073    
DistWater          -0.006833   0.000907  -7.534 1.32e-10 ***
DistOpen            0.093781   0.035616   2.633   0.0104 *  
StandSpBirch       -0.801524   0.880568  -0.910   0.3658    
StandSpForest pine  0.495872   0.742192   0.668   0.5063    
StandSpOak         -0.134608   0.751640  -0.179   0.8584    
StandSpOther        1.429348   0.657085   2.175   0.0330 *  
StandSpSpruce       0.607052   0.666613   0.911   0.3656    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for Tweedie family taken to be 3.595696)

    Null deviance: 301.31  on 78  degrees of freedom
Residual deviance: 207.58  on 70  degrees of freedom
  (12 observations deleted due to missingness)
AIC: NA

Number of Fisher Scoring iterations: 30

But I'm a bit worried about the dispersion parameter, is it normal to be that "high"?

And for other species i get errors and warnings:

> glmnycnoc <- glm(NYCNOCactivity ~ StandAge + DistWater + DistOpen + StandSp, data = flagermusdata, family=tweedie(var.power=2, link.power=0), maxit = 100)
Error in glm.fit(x = c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,  : 
  NA/NaN/Inf in 'x'
In addition: Warning message:
step size truncated due to divergence

If i change the var.power to a number between 1 and 1.5, the test goes fine, but how do i know what var.power to use?

> glmnycnoc <- glm(NYCNOCactivity ~ StandAge + DistWater + DistOpen + StandSp, data = flagermusdata, family=tweedie(var.power= 1.2, link.power=0), maxit = 100)
> summary(glmnycnoc)

Call:
glm(formula = NYCNOCactivity ~ StandAge + DistWater + DistOpen + 
    StandSp, family = tweedie(var.power = 1.2, link.power = 0), 
    data = flagermusdata, maxit = 100)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-3.2779  -1.4789  -0.9516  -0.2899   8.4236  

Coefficients:
                    Estimate Std. Error t value Pr(>|t|)   
(Intercept)         2.297125   0.742437   3.094  0.00284 **
StandAge           -0.003797   0.005683  -0.668  0.50631   
DistWater          -0.005397   0.002928  -1.843  0.06956 . 
DistOpen           -0.072798   0.063763  -1.142  0.25747   
StandSpBirch       -1.857287   1.670818  -1.112  0.27011   
StandSpForest pine -0.571645   0.851399  -0.671  0.50416   
StandSpOak         -1.682813   1.251839  -1.344  0.18320   
StandSpOther       -1.267486   0.914618  -1.386  0.17020   
StandSpSpruce      -1.852949   1.251323  -1.481  0.14315   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for Tweedie family taken to be 8.034003)

    Null deviance: 421.28  on 78  degrees of freedom
Residual deviance: 309.21  on 70  degrees of freedom
  (12 observations deleted due to missingness)
AIC: NA

Number of Fisher Scoring iterations: 8

Thank you for your time and help.

Gavin Simpson · Accepted Answer · 2021-04-16 20:00:17Z

1

I think you'd be better off leaving the data as counts so you can fit a count-based model and instead correct for the differential recording effort via an offset. Say you have variables count and effort representing the original count data and the recording hours, you could fit the model with an offset using

glm(count ~ x1 + x2 + offset(log(effort)), data = foo, family = poisson())

The offset has the effect of scaling the data without you actually having to do. It turns the model from one of counts to one of counts per unit effort.

If you want to continue in the approach you have started with the Tweedie, perhaps consider checking the model by fitting it using mgcv::gam() and it's tw() family, which can estimate the power parameter while fitting. The syntax follows that of glm():

library('mgcv')
gam(NYCNOCactivity ~ StandAge + DistWater + DistOpen + StandSp,
    data = flagermusdata, family = tw, method = "ML")

You could also use the nicer negative binomial family in gam() which estimates the overdispersion parameter while fitting:

gam(count ~ x1 + x2 + offset(log(effort)), data = foo, family = nb(),
    method = "ML")

edited Apr 16, 2021 at 20:00

answered Apr 16, 2021 at 5:27

Gavin Simpson

45.8k7 gold badges127 silver badges178 bronze badges

1

$\begingroup$ That actually make good sense to use the offset function. Thank you for your reply, much appreciated. $\endgroup$
– Peter Sivertsen
Apr 16, 2021 at 7:53
$\begingroup$ one question, why do you use offset(log(effort)) and not just offset(effort) ? $\endgroup$
– Peter Sivertsen
Apr 16, 2021 at 8:00
$\begingroup$ I forgot to add a count family; rectified now. In the count families, the link used is typically the log link, and th offset variable has to be entered on the scale of the link function. If you changed to use the sqrt link you would have to use offset(sqrt(effort)) $\endgroup$
– Gavin Simpson
Apr 16, 2021 at 14:14
$\begingroup$ Once again, thank you! I tried to do the gam() with nb, but is it "bad" that the outcome gives a negative R-sq value, and is only 19.2% explained deviance to be expected? R-sq.(adj) = -0.0404 Deviance explained = 19.2% -REML = 194.5 Scale est. = 1 n = 79 $\endgroup$
– Peter Sivertsen
Apr 16, 2021 at 14:31
$\begingroup$ I also tried the glm() with family = Poisson(link="log") but that just returns the lowest possible p-value for all variables, so I think that my data does not fit that.? $\endgroup$
– Peter Sivertsen
Apr 16, 2021 at 14:38

| Show 2 more comments

MBorg · Accepted Answer · 2022-05-01 09:51:36Z

Partly linked to Gavin's answer. You can also estimate the var power for a GLM should you not wish to use a GAM, using the statmod package. However, estimating the shape parameter process is faster, simpler and more stable with a GAM, particularly for the full range of values from 1.01 to 1.99.

The lowest shape parameter value that can be estimated with either model is 1.01. For statmod, you need to add a number to the 3rd decimal place to consider 1.01 as a potential value e.g. 1.011.

library(statmod) # Tweedie distibution for glm
library(tweedie) # tweedie.profile() to estimate var power

form <- 'formula = PIPNATactivity ~ StandAge + DistWater + DistOpen + StandSp'

# save selected var power
var.p <- tweedie.profile(formula=as.formula(form), data=flagermusdata, p.vec=c(1.011,1.1,1.2,1.3,1.4,1.5), method='series', do.ci=F) # series method is same as that used by gam()
# model with selected var.p
glm(formula=as.formula(form), family=tweedie(var.power=var.p, link.power=0), data=flagermusdata, maxit=100) # model with selected var.p. link.power=0 specifies a log-link.

Stack Exchange Network

Tweedie GLM, what var.power should I chose?

2 Answers 2

Your Answer

Not the answer you're looking for? Browse other questions tagged
generalized-linear-model
tweedie-distribution
or ask your own question.

Hot Network Questions

Tweedie GLM, what var.power should I chose?

2 Answers 2

Your Answer

Sign up or log in

Post as a guest

Not the answer you're looking for? Browse other questions tagged generalized-linear-modeltweedie-distribution or ask your own question.

Related

Hot Network Questions

Not the answer you're looking for? Browse other questions tagged
generalized-linear-model
tweedie-distribution
or ask your own question.