1
$\begingroup$

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
$\endgroup$
4
  • 1
    $\begingroup$ Can you post your data? $\endgroup$ Jan 27, 2017 at 21:59
  • $\begingroup$ Added data to the original post, hope that's the right place to put it $\endgroup$
    – Jessica
    Jan 28, 2017 at 0:48
  • $\begingroup$ Before doing analysis, I used relevel to set the reference for Disease as "Control" and "Blood" as reference for the Source $\endgroup$
    – Jessica
    Jan 28, 2017 at 0:52
  • $\begingroup$ From memory, for some reason gamlss uses NBI to refer to the distribution commonly called NB2 (quadratic mean-variance relationship), and NBII for the distribution commonly called NB1 (linear mean-variance relationship). glm.nb fits NB2 only. Try using family = NBI in gamlss. $\endgroup$
    – Mark
    May 11, 2017 at 23:50

0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.