# Tweedie GLM, what var.power should I chose?

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'
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.

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")

• That actually make good sense to use the offset function. Thank you for your reply, much appreciated. Apr 16, 2021 at 7:53
• one question, why do you use offset(log(effort)) and not just offset(effort) ? Apr 16, 2021 at 8:00
• 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)) Apr 16, 2021 at 14:14
• 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 Apr 16, 2021 at 14:31
• 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.? Apr 16, 2021 at 14:38

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