i have fit negative binomial model to my data as follows:

```
> ngbinmodel <- glm.nb( seizure.rate ~ age + treatment, data = epilepsy_reduced)
>     summary(ngbinmodel)

Call:
glm.nb(formula = seizure.rate ~ age + treatment, data = epilepsy_reduced, 
    init.theta = 1.498983674, link = log)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-2.3510  -0.8790  -0.4563   0.4328   1.8916  

Coefficients:
              Estimate Std. Error z value Pr(>|z|)    
(Intercept)  2.0985089  0.5845392   3.590 0.000331 ***
age         -0.0007965  0.0193064  -0.041 0.967092    
treatment   -0.5011593  0.2405658  -2.083 0.037228 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for Negative Binomial(1.499) family taken to be 1)

    Null deviance: 71.217  on 57  degrees of freedom
Residual deviance: 66.875  on 55  degrees of freedom
AIC: 341.12

Number of Fisher Scoring iterations: 1

              Theta:  1.499 
          Std. Err.:  0.362 

 2 x log-likelihood:  -333.119 

```

Now I would like to check if i should include the interaction effect between age and treatment. I have found two methods to do it: 

```
  > intearaction_nbm<-addterm(ngbinmodel, . ~ . * age,test="Chisq")
            >     summary(intearaction_nbm)
                   Df         AIC             LRT            Pr(Chi)      
             Min.   :1   Min.   :339.1   Min.   :0.9383   Min.   :0.3327  
             1st Qu.:1   1st Qu.:339.4   1st Qu.:0.9383   1st Qu.:0.3327  
             Median :1   Median :339.6   Median :0.9383   Median :0.3327  
             Mean   :1   Mean   :339.6   Mean   :0.9383   Mean   :0.3327  
             3rd Qu.:1   3rd Qu.:339.9   3rd Qu.:0.9383   3rd Qu.:0.3327  
             Max.   :1   Max.   :340.2   Max.   :0.9383   Max.   :0.3327  
             NA's   :1                   NA's   :1        NA's   :1

```

and 

```
> ngbinmodel_int <- glm.nb( seizure.rate ~ age*treatment, data = epilepsy_reduced)
>     summary(ngbinmodel_int)
glm.nb(formula = seizure.rate ~ age * treatment, data = epilepsy_reduced, 
    init.theta = 1.531539174, link = log)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-2.3503  -0.8742  -0.3848   0.3403   1.8508  

Coefficients:
              Estimate Std. Error z value Pr(>|z|)  
(Intercept)    1.51361    0.83920   1.804   0.0713 .
age            0.01914    0.02826   0.677   0.4981  
treatment      0.60748    1.12199   0.541   0.5882  
age:treatment -0.03893    0.03850  -1.011   0.3119  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for Negative Binomial(1.5315) family taken to be 1)

    Null deviance: 72.238  on 57  degrees of freedom
Residual deviance: 66.874  on 54  degrees of freedom
AIC: 342.18

Number of Fisher Scoring iterations: 1

              Theta:  1.532 
          Std. Err.:  0.373 

 2 x log-likelihood:  -332.180 

```

I was expecting to obtain the same result from both of the methods. 

**How can i access the regression estimates of `intearaction_nbm`?**

**why are the outcomes different?** According to `intearaction_nbm` i should include the interaction term (the AIC is lower) but according to `ngbinmodel_int` i should not include the interaction term (AIC increases). 

**would discretizing my continuous variable `age` be advised?**