> mod=glm(y~offset(log(years))+as.factor(gender)+age,family=Gamma(link="log"),
+         data=pm,control = glm.control(maxit = 50))
> shape=gamma.shape(mod)
> summary(mod,dispersion = 1/shape$alpha)

Call:
glm(formula = y ~ offset(log(years)) + as.factor(gender) + age, 
    family = Gamma(link = "log"), data = pm, control = glm.control(maxit = 50))

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-3.8207  -1.2145  -0.5334   0.1910  15.1410  

Coefficients:
                    Estimate Std. Error z value Pr(>|z|)    
(Intercept)        4.6730931  0.0134128   348.4   <2e-16 ***
as.factor(gender)M 0.7806667  0.0024625   317.0   <2e-16 ***
age                0.0642592  0.0001908   336.8   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for Gamma family taken to be 1.238619)

    Null deviance: 1519880  on 852449  degrees of freedom
Residual deviance: 1251784  on 852447  degrees of freedom
AIC: 20497381

Number of Fisher Scoring iterations: 8

> summary(mod)

Call:
glm(formula = y ~ offset(log(years)) + as.factor(gender) + age, 
    family = Gamma(link = "log"), data = pm, control = glm.control(maxit = 50))

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-3.8207  -1.2145  -0.5334   0.1910  15.1410  

Coefficients:
                    Estimate Std. Error t value Pr(>|t|)    
(Intercept)        4.6730931  0.0200320   233.3   <2e-16 ***
as.factor(gender)M 0.7806667  0.0036777   212.3   <2e-16 ***
age                0.0642592  0.0002849   225.5   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for Gamma family taken to be 2.762759)

    Null deviance: 1519880  on 852449  degrees of freedom
Residual deviance: 1251784  on 852447  degrees of freedom
AIC: 20497381

Number of Fisher Scoring iterations: 8

> pscl::pR2(mod)
          llh       llhNull            G2      McFadden          r2ML          r2CU 
-1.024869e+07 -1.034377e+07  1.901660e+05  9.192299e-03  1.999506e-01  1.999506e-01 
> drop1(mod)
Single term deletions

Model:
y ~ offset(log(years)) + as.factor(gender) + age
                  Df Deviance      AIC
<none>                1251784 20497381
as.factor(gender)  1  1367517 20539269
age                1  1383277 20544973
> exp(confint(mod))
Waiting for profiling to be done...
                        2.5 %     97.5 %
(Intercept)        103.423488 110.819024
as.factor(gender)M   2.167202   2.198753
age                  1.065840   1.066889

The confidence interval from confint is definitely not correct. I will need to recalculate based on the new standard erorrs.

                    age, family=Gamma(link="log"),
            data=pm, control = glm.control(maxit = 50))
     shape = gamma.shape(mod)
     summary(mod, dispersion = 1/shape$alpha)
    
    Call:
    glm(formula = y ~ offset(log(years)) + as.factor(gender) + 
               age, 
        family = Gamma(link = "log"), data = pm, 
         control = glm.control(maxit = 50))
    
    Deviance Residuals: 
        Min       1Q   Median       3Q      Max  
    -3.8207  -1.2145  -0.5334   0.1910  15.1410  
    
    Coefficients:
                        Estimate Std. Error z value Pr(>|z|)    
    (Intercept)        4.6730931  0.0134128   348.4   <2e-16 ***
    as.factor(gender)M 0.7806667  0.0024625   317.0   <2e-16 ***
    age                0.0642592  0.0001908   336.8   <2e-16 ***
    ---
    Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
    
    (Dispersion parameter for Gamma family taken to be 1.238619)
    
        Null deviance: 1519880  on 852449  degrees of freedom
    Residual deviance: 1251784  on 852447  degrees of freedom
    AIC: 20497381
    
    Number of Fisher Scoring iterations: 8
    
     summary(mod)
    
    Call:
    glm(formula = y ~ offset(log(years)) + as.factor(gender) + 
          age, 
        family = Gamma(link = "log"), data = pm, 
        control = glm.control(maxit = 50))
    
    Deviance Residuals: 
        Min       1Q   Median       3Q      Max  
    -3.8207  -1.2145  -0.5334   0.1910  15.1410  
    
    Coefficients:
                        Estimate Std. Error t value Pr(>|t|)    
    (Intercept)        4.6730931  0.0200320   233.3   <2e-16 ***
    as.factor(gender)M 0.7806667  0.0036777   212.3   <2e-16 ***
    age                0.0642592  0.0002849   225.5   <2e-16 ***
    ---
    Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
    
    (Dispersion parameter for Gamma family taken to be 2.762759)
    
        Null deviance: 1519880  on 852449  degrees of freedom
    Residual deviance: 1251784  on 852447  degrees of freedom
    AIC: 20497381
    
    Number of Fisher Scoring iterations: 8
    
    > pscl::pR2(mod)
              llh       llhNull            G2      McFadden          r2ML          r2CU 
    -1.024869e+07 -1.034377e+07  1.901660e+05  9.192299e-03  1.999506e-01  1.999506e-01 
    > drop1(mod)
    Single term deletions
    
    Model:
    y ~ offset(log(years)) + as.factor(gender) + age
                      Df Deviance      AIC
    <none>                1251784 20497381
    as.factor(gender)  1  1367517 20539269
    age                1  1383277 20544973
    > exp(confint(mod))
    Waiting for profiling to be done...
                            2.5 %     97.5 %
    (Intercept)        103.423488 110.819024
    as.factor(gender)M   2.167202   2.198753
    age                  1.065840   1.066889

The confidence interval from confint is definitely not correct. I will need to recalculate based on the new standard errors.

Edit3: I realized there is no way to adjusted for MLE dispersion when I used drop1 to get likelihood p.

Edit2: I am looking for a more detailed answer with references and perhaps formula.