I have a dataset consisting of number of mutations per person. Covariates I think may be related are age and source of the DNA sample. I want to determine whether the disease the person has significantly affects the the number of mutations they have.
I determined using these two guides (http://www.ats.ucla.edu/stat/r/dae/poissonreg.htm, Is there a test to determine whether GLM overdispersion is significant?) that the Poisson model is not a good fit to the data because of overdispersion, so I moved onto negative binomial regression.
First, I ran glm.nb(Mutation ~ Disease + Source + Age, data=newdata, maxit = 500)
which yielded:
Estimate SE z-value Pr(>|z|)
(Intercept) 0.15 0.56 0.26 7.94E-01
Disease1 -0.08 0.40 -0.20 8.43E-01
Disease2 1.75 0.50 3.51 4.44E-04
Disease3 2.04 0.53 3.85 1.20E-04
Disease4 0.54 0.44 1.23 2.18E-01
Disease5 -0.42 0.55 -0.76 4.45E-01
Sourcefib 0.64 0.42 1.54 1.25E-01
SouceLCL 0.35 0.25 1.40 1.60E-01
Sourceunknown 0.37 0.67 0.56 5.79E-01
Age 0.03 0.01 3.50 4.69E-04
Unfortunately, when I tried to generate confidence intervals using confint(, maxit=500), I got the error Warning message: glm.fit: algorithm did not converge
(Edit: it turns out I was able to get around this by fixing a typo/data coding error, but my original question remains).
I searched but could not find a solution to this convergence error, so I found two other alternative R packages that perform negative binomial regression, the ml.nb2
function in the COUNT
package and gamlss
in the gamlss
package.
ml.nb2(Mutation ~ Disease + Source + Age, data=newdata)
Estimate SE Z LCL UCL p-value
(Intercept) 0.75 0.62 1.20 -0.48 1.97 2.32E-01
Disease1 -0.16 0.42 -0.37 -0.99 0.67 7.12E-01
Disease2 1.42 0.48 2.93 0.47 2.37 3.38E-03
Disease3 1.37 0.62 2.19 0.15 2.59 2.83E-02
Disease4 0.64 0.48 1.34 -0.29 1.57 1.79E-01
Disease5 -0.50 0.59 -0.85 -1.65 0.65 3.94E-01
Sourcefib 0.05 0.41 0.13 -0.75 0.85 9.00E-01
SouceLCL 0.09 0.27 0.33 -0.45 0.63 7.42E-01
Sourceunknown -1.47 0.62 -2.37 -2.68 -0.25 1.79E-02
Age 0.02 0.01 2.40 0.00 0.04 1.66E-02
alpha 0.94 0.16 5.97 0.63 1.25 2.44E-09
gamlss(Mutation ~ Disease + Source + Age, data=newdata, family=NBII)
Estimate SE t-value Pr(>|t|)
(Intercept) 1.15 0.51 2.24 2.80E-02
Disease1 -0.12 0.33 -0.36 7.24E-01
Disease2 0.94 0.41 2.30 2.41E-02
Disease3 0.41 0.42 0.98 3.32E-01
Disease4 -0.04 0.37 -0.10 9.24E-01
Disease5 -0.56 0.50 -1.12 2.66E-01
Sourcefib 0.62 0.34 1.85 6.90E-02
SouceLCL 0.48 0.23 2.12 3.72E-02
Sourceunknown 0.60 0.51 1.17 2.45E-01
Age 0.02 0.01 2.26 2.67E-02
Why do these different implementations give me such different estimates for the coefficients, especially Disease2 and Disease3, and also differ (particularly gamlss
) in assessing which covariates are significant?
Data:
Disease Source Mutation Age
Disease4 Blood 47 47
Disease4 Blood 6 40
Disease4 Blood 10 39
Disease4 Blood 5 31
Disease4 Blood 1 36
Disease4 Blood 2 59
Disease4 Blood 0 33
Disease4 Blood 1 28
Disease5 Blood 1 58
Disease5 Blood 19 52
Disease5 Blood 0 54
Disease5 Blood 1 52
Disease5 Blood 1 50
Disease5 Blood 1 52
Control Blood 10 60
Disease3 Blood 14 12
Disease3 Blood 1 30
Disease3 Blood 65 10
Disease3 Blood 3 62
Disease2 Blood 2 37
Disease2 Blood 4 25
Disease2 Blood 6 37
Disease2 Blood 3 36
Disease2 Blood 13 16
Disease1 LCL 9 52
Disease1 LCL 7 55
Disease1 LCL 16 35
Disease1 LCL 3 51
Disease1 LCL 8 26
Disease1 LCL 5 49
Disease1 LCL 8 45
Disease1 LCL 1 38
Disease1 LCL 4 52
Disease1 LCL 2 39
Disease1 LCL 1 21
Disease4 LCL 4 41
Disease4 LCL 1 25
Disease4 LCL 4 46
Disease4 LCL 45 48
Disease4 LCL 3 35
Disease4 LCL 6 41
Disease4 LCL 3 NA
Disease4 LCL 1 34
Disease5 LCL 3 NA
Disease5 LCL 3 29
Control LCL 2 22
Control LCL 6 48
Control LCL 2 19
Control LCL 7 54
Control LCL 3 NA
Control LCL 2 23
Control LCL 20 51
Disease3 LCL 9 39
Disease3 LCL 3 10
Disease3 LCL 9 28
Disease2 LCL 18 23
Disease2 LCL 194 47
Disease2 LCL 5 28
Disease1 Blood 5 56
Disease1 Blood 1 30
Disease1 Blood 2 24
Disease1 Blood 3 25
Disease1 Blood 4 61
Disease1 Blood 2 41
Disease1 Blood 3 36
Disease1 Blood 3 40
Disease1 Blood 4 39
Disease1 Blood 12 42
Disease1 Blood 1 15
Disease1 Blood 1 38
Disease1 Blood 0 15
Disease1 Blood 2 50
Disease1 Blood 4 62
Disease1 Blood 1 41
Disease1 Blood 5 37
Disease1 Blood 17 48
Disease1 Blood 2 26
Control Blood 8 56
Control Blood 2 52
Control Blood 1 NA
Control Blood 0 NA
Control Blood 2 NA
Control Blood 4 NA
Control Blood 2 NA
Control Blood 12 NA
Disease1 fib 3 40
Disease1 fib 1 18
Disease1 fib 4 58
Disease1 fib 10 36
Disease1 fib 7 24
Disease1 fib 7 27
Disease1 fib 15 58
Disease2 fib 23 17
Disease1 unknown 4 19
Disease1 unknown 3 17
Disease1 unknown 4 50
gamlss
usesNBI
to refer to the distribution commonly called NB2 (quadratic mean-variance relationship), andNBII
for the distribution commonly called NB1 (linear mean-variance relationship).glm.nb
fits NB2 only. Try usingfamily = NBI
ingamlss
. $\endgroup$