統計ソフトRでAPCモデル分析をするには?

にほんブログ村 科学ブログ 数学へ

年齢・時代・コホート(age-period-cohort、APC)モデルは、 病気の原因を年齢と時代とコホートとに分けて分析する方法だ。

コホートは、同世代とか同級生みたいな意味。

年齢と時代とコホートは、 ほとんど同じようなものなのに、 同時に分析できるのか?

多重共線性の問題が気になるこの3つを同時に扱えるのが、 APCモデルだ。

多重共線性についてはこちら

参照したのは、こちらのページ。

また参考書籍はこちら。 医学データ―デザインから統計モデルまで

使用するデータは年齢区分、時代区分、コホート区分ごとの 肺がん死亡者数と観察人年のデータ。

年齢・時代(AP)モデル、 年齢・コホート(AC)モデル、 年齢・時代・コホートAPC)モデル の順に示していく。

使用する統計モデルは、 ポアソン回帰モデルだ。 対数線形モデルともいう。

ポアソン回帰モデルのポアソンは、 ポアソン分布から来ている。

ポアソン分布とは、 まれにしか起きないイベントの確率分布。

それも0回、1回、2回、3回起こるなど、 離れている確率分布。

離散確率分布という。

ちなみにポアソンは、 ポアソン分布を発表した数学者の名前。

統計ソフトRでAPモデル

年齢・時代区分でデータを概観する

まずデータを読み込んで、 データセットの構造を見てみる。

dat <- read.table("http://toukeier.net/data/lung5-M.txt",header=T)

str(dat)

データは、110件。 Aが年齢区分、 Pが時代区分、 Dが肺がん死亡者数、 Yが観察人年。

> str(dat)
'data.frame':   110 obs. of  4 variables:
 $ A: int  40 40 40 40 40 40 40 40 40 40 ...
 $ P: int  1943 1948 1953 1958 1963 1968 1973 1978 1983 1988 ...
 $ D: int  80 81 73 99 82 97 86 90 116 149 ...
 $ Y: num  694047 754770 769441 749265 757240 ...

データの要約を見てみる。

summary(dat)

年齢は40から85、 時代は1943年から1993年まで、 死亡は7人から2293人、 観察人年は24,656が最小で、1,026,461が最大。

> summary(dat)
       A              P              D                Y          
 Min.   :40.0   Min.   :1943   Min.   :   7.0   Min.   :  24656  
 1st Qu.:50.0   1st Qu.:1953   1st Qu.: 169.2   1st Qu.: 234051  
 Median :62.5   Median :1968   Median : 444.0   Median : 487609  
 Mean   :62.5   Mean   :1968   Mean   : 692.5   Mean   : 461098  
 3rd Qu.:75.0   3rd Qu.:1983   3rd Qu.:1086.2   3rd Qu.: 682706  
 Max.   :85.0   Max.   :1993   Max.   :2293.0   Max.   :1026461  

年齢・時代区分ごとの合計観察人年を集計する。

table(dat$A)
table(dat$P)
tapply(dat$Y, list(dat$A, dat$P), sum)

年齢区分ごと11件のデータ、 時代区分ごと10件のデータ、 であることがわかる。

年齢・時代区分で、観察人年を合計してみる。

> table(dat$A)

40 45 50 55 60 65 70 75 80 85 
11 11 11 11 11 11 11 11 11 11 

> table(dat$P)

1943 1948 1953 1958 1963 1968 1973 1978 1983 1988 1993 
  10   10   10   10   10   10   10   10   10   10   10 

> tapply(dat$Y, list(dat$A, dat$P), sum)
        1943      1948      1953      1958      1963      1968      1973     1978     1983       1988     1993
40 694046.50 754769.50 769440.67 749264.50 757240.00 709558.50 695210.17 756263.0 941402.0 1026461.33 752950.2
45 622256.67 676718.00 738290.50 754357.67 737405.67 747054.83 697976.33 681063.8 741553.2  924431.83 821389.8
50 538964.17 600506.33 653867.50 715819.83 733590.17 717677.33 724880.33 675371.5 659498.3  719660.00 700860.2
55 471016.00 512338.00 571270.67 622413.33 681097.00 699103.17 683242.67 686939.5 640796.8  626481.33 544088.2
60 403172.50 435098.33 474197.50 528106.33 573204.83 627036.33 644142.67 627509.2 630384.5  590669.00 463087.8
65 328690.50 357694.83 386083.00 419562.00 463265.17 501020.00 548399.50 564173.5 548586.2  553420.50 421467.7
70 230090.83 269235.83 294786.67 317388.00 341288.33 373577.00 404348.83 442925.0 458828.3  448996.33 365932.7
75 140110.67 166641.83 195729.83 214930.33 228793.50 245932.00 268415.17 290162.3 319152.7  336461.67 262910.2
80  67778.83  80587.00  98561.33 116116.67 125697.33 136646.17 150131.83 163433.0 175767.2  196465.83 168011.7
85  24656.17  28463.83  34280.50  42136.33  49263.33  56018.17  63742.67  71226.5  77621.5   85376.67  74610.5

年齢・時代区分で、肺がん死亡者数、肺がん死亡率を計算してみる。

tapply(dat$D,list(dat$A,dat$P),sum)

dat$R <- dat$D/dat$Y*10^5

round(tapply(dat$R,list(dat$A,dat$P),sum),1)

結果はこちら。

> tapply(dat$D,list(dat$A,dat$P),sum)
   1943 1948 1953 1958 1963 1968 1973 1978 1983 1988 1993
40   80   81   73   99   82   97   86   90  116  149   91
45  135  163  208  226  252  284  263  251  257  265  251
50  197  292  442  508  560  580  657  608  591  493  446
55  261  404  596  772 1052 1075 1115 1218 1090  995  696
60  213  394  577  955 1342 1682 1654 1826 1885 1497 1113
65  141  273  491  868 1235 1856 2136 2231 2188 2193 1485
70  110  215  300  596  976 1448 1924 2283 2293 2157 1691
75   54  126  167  320  514  860 1213 1559 1824 1640 1221
80   20   57   87  157  220  390  573  753  881  837  716
85    7   10   23   48   72  110  176  213  307  286  262

> round(tapply(dat$R,list(dat$A,dat$P),sum),1)
   1943 1948  1953  1958  1963  1968  1973  1978  1983  1988  1993
40 11.5 10.7   9.5  13.2  10.8  13.7  12.4  11.9  12.3  14.5  12.1
45 21.7 24.1  28.2  30.0  34.2  38.0  37.7  36.9  34.7  28.7  30.6
50 36.6 48.6  67.6  71.0  76.3  80.8  90.6  90.0  89.6  68.5  63.6
55 55.4 78.9 104.3 124.0 154.5 153.8 163.2 177.3 170.1 158.8 127.9
60 52.8 90.6 121.7 180.8 234.1 268.2 256.8 291.0 299.0 253.4 240.3
65 42.9 76.3 127.2 206.9 266.6 370.4 389.5 395.4 398.8 396.3 352.3
70 47.8 79.9 101.8 187.8 286.0 387.6 475.8 515.4 499.8 480.4 462.1
75 38.5 75.6  85.3 148.9 224.7 349.7 451.9 537.3 571.5 487.4 464.4
80 29.5 70.7  88.3 135.2 175.0 285.4 381.7 460.7 501.2 426.0 426.2
85 28.4 35.1  67.1 113.9 146.2 196.4 276.1 299.0 395.5 335.0 351.2

グラフにして眺めてみる。

まず、年齢区分別の時代区分と肺がん死亡率のグラフ。

「通りすがりの人」に教わった賢いスクリプト

split()とかlappy()とかなかなか使えないfunction。

勉強になる。

colは数字では1から8までしか指定できないので、 9番目と10番目は自分で色を選んだ。

plot(dat$P, dat$R, ylim=c(10, 600), type="n", xlab="Period",ylab="Lung Cancer Mortality Rate")
dat2 <- split(dat, dat$A)
colset <- c(1:8,"orange","yellowgreen")
lapply(dat2, function(d) lines(d$P, d$R, type="o", col=colset[d$A/5-7]))
legend(1943, 600, paste("Aged",seq(40,85,5)), col=colset, lty=1)

また、逆に時代区分別の年齢区分と肺がん死亡率のグラフを描く。

1970年代後半から、75歳を超えると肺がん死亡率は減少する という傾向が続いている。

plot(dat$A, dat$R, xlim=c(40, 85), ylim=c(10, 600), type="n", xlab="Age",ylab="Lung Cancer Mortality Rate")
dat2 <- split(dat, dat$P)
colset <- c(1:8,"orange","yellowgreen","navy")
lapply(dat2, function(d) lines(d$A, d$R, type="o", col=colset[(d$P-3)/5-387]))
legend(42, 600, seq(1943,1993,5), col=colset, lty=1)

ポアソン回帰モデルでAPモデルを分析する

AとPを因子型でglm()に投入する。 family=poissonがポアソン回帰モデルの指定。 観察人年の対数をoffset()で指定する。

ap.1 <- glm(D ~ factor(A) + factor(P) + offset(log(Y)), family=poisson, data=dat)
summary(ap.1)

40歳の年齢区分と、 1943年の時代区分をそれぞれ参照カテゴリにして、 相対リスクが対数で計算されている。

> summary(ap.1)

Call:
glm(formula = D ~ factor(A) + factor(P) + offset(log(Y)), family = poisson, 
    data = dat)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-10.400   -3.728   -0.984    3.685   11.203  

Coefficients:
               Estimate Std. Error z value Pr(>|z|)    
(Intercept)   -10.34235    0.04192 -246.71   <2e-16 ***
factor(A)45     0.95258    0.03673   25.93   <2e-16 ***
factor(A)50     1.78237    0.03383   52.69   <2e-16 ***
factor(A)55     2.41412    0.03265   73.94   <2e-16 ***
factor(A)60     2.86259    0.03216   89.01   <2e-16 ***
factor(A)65     3.15159    0.03201   98.47   <2e-16 ***
factor(A)70     3.31784    0.03209  103.40   <2e-16 ***
factor(A)75     3.30980    0.03261  101.50   <2e-16 ***
factor(A)80     3.17640    0.03423   92.81   <2e-16 ***
factor(A)85     2.90983    0.04024   72.32   <2e-16 ***
factor(P)1948   0.39206    0.03629   10.80   <2e-16 ***
factor(P)1953   0.67592    0.03404   19.86   <2e-16 ***
factor(P)1958   1.01434    0.03226   31.44   <2e-16 ***
factor(P)1963   1.26666    0.03130   40.47   <2e-16 ***
factor(P)1968   1.48717    0.03067   48.49   <2e-16 ***
factor(P)1973   1.59239    0.03039   52.40   <2e-16 ***
factor(P)1978   1.67994    0.03020   55.62   <2e-16 ***
factor(P)1983   1.69902    0.03015   56.35   <2e-16 ***
factor(P)1988   1.59958    0.03028   52.83   <2e-16 ***
factor(P)1993   1.52558    0.03078   49.57   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 71776.2  on 109  degrees of freedom
Residual deviance:  2723.5  on  90  degrees of freedom
AIC: 3620.5

Number of Fisher Scoring iterations: 5

Formula内に-1を追加すると切片がないモデルになり、 年齢区分Aの係数の推定値が死亡率の絶対値に変わる。

ap.2 <- glm(D ~ factor(A) -1 + factor(P) + offset(log(Y)), family=poisson, data=dat)
summary(ap.2)

Aの係数が0から3程度だったものが、 -7から-10程度に変わった。

> summary(ap.2)

Call:
glm(formula = D ~ factor(A) - 1 + factor(P) + offset(log(Y)), 
    family = poisson, data = dat)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-10.400   -3.728   -0.984    3.685   11.203  

Coefficients:
               Estimate Std. Error z value Pr(>|z|)    
factor(A)40   -10.34235    0.04192 -246.71   <2e-16 ***
factor(A)45    -9.38977    0.03454 -271.89   <2e-16 ***
factor(A)50    -8.55998    0.03145 -272.17   <2e-16 ***
factor(A)55    -7.92822    0.03020 -262.48   <2e-16 ***
factor(A)60    -7.47976    0.02970 -251.83   <2e-16 ***
factor(A)65    -7.19075    0.02956 -243.26   <2e-16 ***
factor(A)70    -7.02451    0.02970 -236.53   <2e-16 ***
factor(A)75    -7.03255    0.03031 -232.05   <2e-16 ***
factor(A)80    -7.16595    0.03209 -223.33   <2e-16 ***
factor(A)85    -7.43252    0.03847 -193.22   <2e-16 ***
factor(P)1948   0.39206    0.03629   10.80   <2e-16 ***
factor(P)1953   0.67592    0.03404   19.86   <2e-16 ***
factor(P)1958   1.01434    0.03226   31.44   <2e-16 ***
factor(P)1963   1.26666    0.03130   40.47   <2e-16 ***
factor(P)1968   1.48717    0.03067   48.49   <2e-16 ***
factor(P)1973   1.59239    0.03039   52.40   <2e-16 ***
factor(P)1978   1.67994    0.03020   55.62   <2e-16 ***
factor(P)1983   1.69902    0.03015   56.35   <2e-16 ***
factor(P)1988   1.59958    0.03028   52.83   <2e-16 ***
factor(P)1993   1.52558    0.03078   49.57   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 1.0037e+08  on 110  degrees of freedom
Residual deviance: 2.7235e+03  on  90  degrees of freedom
AIC: 3620.5

Number of Fisher Scoring iterations: 5

時代区分Pの参照カテゴリを1968の区分にする。 relevel()を使う。

ap.3 <- glm(D ~ factor(A) -1 + relevel(factor(P),"1968") + offset(log(Y)), family=poisson, data=dat)
summary(ap.3)

1968をゼロとして、 計算しなおされている。

> summary(ap.3)

Call:
glm(formula = D ~ factor(A) - 1 + relevel(factor(P), "1968") + 
    offset(log(Y)), family = poisson, data = dat)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-10.400   -3.728   -0.984    3.685   11.203  

Coefficients:
                               Estimate Std. Error  z value Pr(>|z|)    
factor(A)40                    -8.85517    0.03267 -271.034  < 2e-16 ***
factor(A)45                    -7.90259    0.02232 -354.007  < 2e-16 ***
factor(A)50                    -7.07280    0.01708 -414.106  < 2e-16 ***
factor(A)55                    -6.44105    0.01455 -442.648  < 2e-16 ***
factor(A)60                    -5.99259    0.01342 -446.388  < 2e-16 ***
factor(A)65                    -5.70358    0.01313 -434.461  < 2e-16 ***
factor(A)70                    -5.53734    0.01338 -413.986  < 2e-16 ***
factor(A)75                    -5.54537    0.01462 -379.299  < 2e-16 ***
factor(A)80                    -5.67877    0.01795 -316.396  < 2e-16 ***
factor(A)85                    -5.94534    0.02776 -214.208  < 2e-16 ***
relevel(factor(P), "1968")1943 -1.48717    0.03067  -48.493  < 2e-16 ***
relevel(factor(P), "1968")1948 -1.09512    0.02481  -44.134  < 2e-16 ***
relevel(factor(P), "1968")1953 -0.81125    0.02137  -37.958  < 2e-16 ***
relevel(factor(P), "1968")1958 -0.47284    0.01842  -25.674  < 2e-16 ***
relevel(factor(P), "1968")1963 -0.22051    0.01667  -13.227  < 2e-16 ***
relevel(factor(P), "1968")1973  0.10522    0.01488    7.071 1.54e-12 ***
relevel(factor(P), "1968")1978  0.19276    0.01449   13.300  < 2e-16 ***
relevel(factor(P), "1968")1983  0.21184    0.01439   14.724  < 2e-16 ***
relevel(factor(P), "1968")1988  0.11241    0.01465    7.670 1.71e-14 ***
relevel(factor(P), "1968")1993  0.03840    0.01566    2.453   0.0142 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 1.0037e+08  on 110  degrees of freedom
Residual deviance: 2.7235e+03  on  90  degrees of freedom
AIC: 3620.5

Number of Fisher Scoring iterations: 5

APモデルの年齢効果をグラフにする

年齢区分の推定値(係数, coef, estimateすべて同じもの) と標準誤差を取り出して、 推定値は黒の太実線で、 95%信頼区間の上限・下限は赤の実線で描画する。

ap.cf <- summary(ap.3)$coef
age.cf <- ap.cf[1:10,1:2]
am <- seq(40,85,5)+2.5
plot(am,exp(age.cf[,1]),type="l",log="y", lwd=3, ylab="Lung cancer mortality rate", xlab="Age")
lines(am,exp(age.cf[,1]-1.96*age.cf[,2]),col="red")
lines(am,exp(age.cf[,1]+1.96*age.cf[,2]),col="red")

APモデルの時代効果をグラフに表す

APモデルの結果から、 時代効果の推定値と標準誤差を取り出し、 参照カテゴリ"1968"に推定値0、標準誤差0を挿入し、 データセットを準備する。

上下対称になるように、0.2(5分の1)から5.0(5倍)にY軸を指定する。

黒い太実線が推定値、 赤い細実線が95%信頼区間の上限と下限。

X軸上、1970.5(1968のカテゴリ1968から1973の中点が1970.5)で、 Y軸1.0と交差しているのがわかる。 exp(Y)=1はY=0。

1968の時代区分を参照カテゴリ=0。 つまり、exp(0)=1としたことが確認できる。

per.cf <- ap.cf[11:20,1:2]
per.cf <- rbind(per.cf[1:5,],c(0,0),per.cf[6:10,])
pm <- seq(1943,1993,5)+2.5
plot(pm, exp(per.cf[,1]),type="l", log="y", lwd=3, ylim=c(0.2,5.0), ylab="Relative risk of lung cancer death for period", xlab="Period")
lines(pm, exp(per.cf[,1]-1.96*per.cf[,2]),col="red")
lines(pm, exp(per.cf[,1]+1.96*per.cf[,2]),col="red")

統計ソフトRでACモデル

年齢・コホート区分でデータを概観する

コホート区分変数を作成する。

年齢・コホート区分ごとの観察人年データを集計する。

dat$C <- dat$P-dat$A
table(dat$C)
tapply(dat$Y, list(dat$A, dat$C), sum)

コホート区分は、20区分できる。

観察人年が計算できない年齢・コホート区分はNAと返される。

データがない区分なので仕方ない。

> table(dat$C)

1858 1863 1868 1873 1878 1883 1888 1893 1898 1903 1908 1913 1918 1923 1928 1933 1938 1943 1948 1953 
   1    2    3    4    5    6    7    8    9   10   10    9    8    7    6    5    4    3    2    1 

> tapply(dat$Y, list(dat$A, dat$C), sum)
       1858     1863     1868      1873      1878      1883      1888     1893     1898      1903     1908     1913     1918     1923     1928     1933     1938
40       NA       NA       NA        NA        NA        NA        NA       NA       NA 694046.50 754769.5 769440.7 749264.5 757240.0 709558.5 695210.2 756263.0
45       NA       NA       NA        NA        NA        NA        NA       NA 622256.7 676718.00 738290.5 754357.7 737405.7 747054.8 697976.3 681063.8 741553.2
50       NA       NA       NA        NA        NA        NA        NA 538964.2 600506.3 653867.50 715819.8 733590.2 717677.3 724880.3 675371.5 659498.3 719660.0
55       NA       NA       NA        NA        NA        NA 471016.00 512338.0 571270.7 622413.33 681097.0 699103.2 683242.7 686939.5 640796.8 626481.3 544088.2
60       NA       NA       NA        NA        NA 403172.50 435098.33 474197.5 528106.3 573204.83 627036.3 644142.7 627509.2 630384.5 590669.0 463087.8       NA
65       NA       NA       NA        NA 328690.50 357694.83 386083.00 419562.0 463265.2 501020.00 548399.5 564173.5 548586.2 553420.5 421467.7       NA       NA
70       NA       NA       NA 230090.83 269235.83 294786.67 317388.00 341288.3 373577.0 404348.83 442925.0 458828.3 448996.3 365932.7       NA       NA       NA
75       NA       NA 140110.7 166641.83 195729.83 214930.33 228793.50 245932.0 268415.2 290162.33 319152.7 336461.7 262910.2       NA       NA       NA       NA
80       NA 67778.83  80587.0  98561.33 116116.67 125697.33 136646.17 150131.8 163433.0 175767.17 196465.8 168011.7       NA       NA       NA       NA       NA
85 24656.17 28463.83  34280.5  42136.33  49263.33  56018.17  63742.67  71226.5  77621.5  85376.67  74610.5       NA       NA       NA       NA       NA       NA
       1943      1948     1953
40 941402.0 1026461.3 752950.2
45 924431.8  821389.8       NA
50 700860.2        NA       NA
55       NA        NA       NA
60       NA        NA       NA
65       NA        NA       NA
70       NA        NA       NA
75       NA        NA       NA
80       NA        NA       NA
85       NA        NA       NA

年齢区分別、コホート別の肺がん死亡者数と肺がん死亡率を計算してみる。

tapply(dat$D,list(dat$A,dat$C),sum)

dat$R <- dat$D/dat$Y*10^5

round(tapply(dat$R,list(dat$A,dat$C),sum),1)

結果はこちら。

> tapply(dat$D,list(dat$A,dat$C),sum)
   1858 1863 1868 1873 1878 1883 1888 1893 1898 1903 1908 1913 1918 1923 1928 1933 1938 1943 1948 1953
40   NA   NA   NA   NA   NA   NA   NA   NA   NA   80   81   73   99   82   97   86   90  116  149   91
45   NA   NA   NA   NA   NA   NA   NA   NA  135  163  208  226  252  284  263  251  257  265  251   NA
50   NA   NA   NA   NA   NA   NA   NA  197  292  442  508  560  580  657  608  591  493  446   NA   NA
55   NA   NA   NA   NA   NA   NA  261  404  596  772 1052 1075 1115 1218 1090  995  696   NA   NA   NA
60   NA   NA   NA   NA   NA  213  394  577  955 1342 1682 1654 1826 1885 1497 1113   NA   NA   NA   NA
65   NA   NA   NA   NA  141  273  491  868 1235 1856 2136 2231 2188 2193 1485   NA   NA   NA   NA   NA
70   NA   NA   NA  110  215  300  596  976 1448 1924 2283 2293 2157 1691   NA   NA   NA   NA   NA   NA
75   NA   NA   54  126  167  320  514  860 1213 1559 1824 1640 1221   NA   NA   NA   NA   NA   NA   NA
80   NA   20   57   87  157  220  390  573  753  881  837  716   NA   NA   NA   NA   NA   NA   NA   NA
85    7   10   23   48   72  110  176  213  307  286  262   NA   NA   NA   NA   NA   NA   NA   NA   NA

> round(tapply(dat$R,list(dat$A,dat$C),sum),1)
   1858 1863 1868  1873  1878  1883  1888  1893  1898  1903  1908  1913  1918  1923  1928  1933  1938 1943 1948 1953
40   NA   NA   NA    NA    NA    NA    NA    NA    NA  11.5  10.7   9.5  13.2  10.8  13.7  12.4  11.9 12.3 14.5 12.1
45   NA   NA   NA    NA    NA    NA    NA    NA  21.7  24.1  28.2  30.0  34.2  38.0  37.7  36.9  34.7 28.7 30.6   NA
50   NA   NA   NA    NA    NA    NA    NA  36.6  48.6  67.6  71.0  76.3  80.8  90.6  90.0  89.6  68.5 63.6   NA   NA
55   NA   NA   NA    NA    NA    NA  55.4  78.9 104.3 124.0 154.5 153.8 163.2 177.3 170.1 158.8 127.9   NA   NA   NA
60   NA   NA   NA    NA    NA  52.8  90.6 121.7 180.8 234.1 268.2 256.8 291.0 299.0 253.4 240.3    NA   NA   NA   NA
65   NA   NA   NA    NA  42.9  76.3 127.2 206.9 266.6 370.4 389.5 395.4 398.8 396.3 352.3    NA    NA   NA   NA   NA
70   NA   NA   NA  47.8  79.9 101.8 187.8 286.0 387.6 475.8 515.4 499.8 480.4 462.1    NA    NA    NA   NA   NA   NA
75   NA   NA 38.5  75.6  85.3 148.9 224.7 349.7 451.9 537.3 571.5 487.4 464.4    NA    NA    NA    NA   NA   NA   NA
80   NA 29.5 70.7  88.3 135.2 175.0 285.4 381.7 460.7 501.2 426.0 426.2    NA    NA    NA    NA    NA   NA   NA   NA
85 28.4 35.1 67.1 113.9 146.2 196.4 276.1 299.0 395.5 335.0 351.2    NA    NA    NA    NA    NA    NA   NA   NA   NA

肺がん死亡率のグラフを描いてみる。

まず、年齢別のコホート区分と肺がん死亡率のグラフ。

コホートによって、まったく様相が異なることがわかる。

1868年(日本なら明治元年)生まれの75歳は肺がん死亡率がとても低いのに比べ、 1908年(日本なら日露戦争直後)生まれの75歳の肺がん死亡率はもっとも高いのがわかる。

plot(dat$C, dat$R, ylim=c(10, 600), type="n", xlab="Cohort",ylab="Lung Cancer Mortality Rate")
dat2 <- split(dat, dat$A)
colset <- c(1:8,"orange","yellowgreen")
lapply(dat2, function(d) lines(d$C, d$R, type="o", col=colset[d$A/5-7]))
legend(1858, 600, paste("Aged",seq(40,85,5)), col=colset, lty=1)

次にコホート区分別の年齢区分と肺がん死亡率のグラフ。

データがある程度そろっている1883年から1933年のコホートに限定した。

1903年生まれから、肺がん死亡率の上昇は鈍化し、 1923年生まれがピークで、1928年生まれから減少に転じた。

dat1 <- subset(dat, dat$C>=1883 & dat$C<=1993)
plot(dat1$A, dat1$R, xlim=c(40, 85), ylim=c(10, 600), type="n", xlab="Age",ylab="Lung Cancer Mortality Rate")
dat2 <- split(dat1, dat1$C)
colset <- c(1:8,"orange","yellowgreen","navy")
lapply(dat2, function(d) lines(d$A, d$R, type="o", col=colset[(d$C-3)/5-375]))
legend(42, 600, seq(1883,1933,5), col=colset, lty=1)

ポアソン回帰モデルでACモデルを分析する

APモデルと同じように、 AとCを因子型でモデルに投入する。 Formulaに-1を追加して、 年齢区分の推定値を絶対値に変換しておく。

ac.2 <- glm(D ~ -1 + factor(A) + factor(C) + offset(log(Y)), family=poisson, data=dat)
summary(ac.2)

結果はこちら。

> summary(ac.2)

Call:
glm(formula = D ~ -1 + factor(A) + factor(C) + offset(log(Y)), 
    family = poisson, data = dat)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-7.2822  -2.0274   0.3573   2.0545   5.2834  

Coefficients:
               Estimate Std. Error z value Pr(>|z|)    
factor(A)40   -11.83501    0.38038 -31.114  < 2e-16 ***
factor(A)45   -10.86658    0.37948 -28.636  < 2e-16 ***
factor(A)50   -10.00034    0.37915 -26.376  < 2e-16 ***
factor(A)55    -9.32333    0.37903 -24.598  < 2e-16 ***
factor(A)60    -8.80577    0.37898 -23.236  < 2e-16 ***
factor(A)65    -8.42760    0.37896 -22.239  < 2e-16 ***
factor(A)70    -8.16176    0.37896 -21.537  < 2e-16 ***
factor(A)75    -8.04871    0.37899 -21.237  < 2e-16 ***
factor(A)80    -8.05099    0.37913 -21.236  < 2e-16 ***
factor(A)85    -8.16687    0.37796 -21.608  < 2e-16 ***
factor(C)1863   0.01046    0.42031   0.025 0.980152    
factor(C)1868   0.51345    0.38845   1.322 0.186240    
factor(C)1873   0.82684    0.38231   2.163 0.030560 *  
factor(C)1878   1.05336    0.38054   2.768 0.005639 ** 
factor(C)1883   1.41904    0.37972   3.737 0.000186 ***
factor(C)1888   1.91197    0.37927   5.041 4.63e-07 ***
factor(C)1893   2.28073    0.37909   6.016 1.78e-09 ***
factor(C)1898   2.55794    0.37900   6.749 1.49e-11 ***
factor(C)1903   2.76315    0.37895   7.292 3.06e-13 ***
factor(C)1908   2.83415    0.37894   7.479 7.48e-14 ***
factor(C)1913   2.81410    0.37901   7.425 1.13e-13 ***
factor(C)1918   2.86228    0.37902   7.552 4.30e-14 ***
factor(C)1923   2.91551    0.37906   7.691 1.45e-14 ***
factor(C)1928   2.86546    0.37917   7.557 4.12e-14 ***
factor(C)1933   2.86314    0.37936   7.547 4.44e-14 ***
factor(C)1938   2.72290    0.37983   7.169 7.57e-13 ***
factor(C)1943   2.68759    0.38066   7.060 1.66e-12 ***
factor(C)1948   2.85099    0.38263   7.451 9.27e-14 ***
factor(C)1953   2.81411    0.39456   7.132 9.87e-13 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 1.0037e+08  on 110  degrees of freedom
Residual deviance: 8.2963e+02  on  81  degrees of freedom
AIC: 1744.7

Number of Fisher Scoring iterations: 4

コホート区分Cの参照カテゴリを1908の区分にする。 端のカテゴリはデータが少なくモデルが不安定なので、 データがしっかりしている中央のカテゴリに変更する。

ac.3 <- glm(D ~ -1 + factor(A) + relevel(factor(C), "1908") + offset(log(Y)), family=poisson, data=dat)
summary(ac.3)

結果はこちら。

> summary(ac.3)

Call:
glm(formula = D ~ -1 + factor(A) + relevel(factor(C), "1908") + 
    offset(log(Y)), family = poisson, data = dat)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-7.2822  -2.0274   0.3573   2.0545   5.2834  

Coefficients:
                               Estimate Std. Error  z value Pr(>|z|)    
factor(A)40                    -9.00086    0.03491 -257.802  < 2e-16 ***
factor(A)45                    -8.03243    0.02319 -346.355  < 2e-16 ***
factor(A)50                    -7.16619    0.01710 -419.037  < 2e-16 ***
factor(A)55                    -6.48918    0.01418 -457.703  < 2e-16 ***
factor(A)60                    -5.97162    0.01277 -467.563  < 2e-16 ***
factor(A)65                    -5.59346    0.01221 -458.001  < 2e-16 ***
factor(A)70                    -5.32761    0.01229 -433.384  < 2e-16 ***
factor(A)75                    -5.21456    0.01348 -386.874  < 2e-16 ***
factor(A)80                    -5.21684    0.01690 -308.730  < 2e-16 ***
factor(A)85                    -5.33272    0.02721 -196.012  < 2e-16 ***
relevel(factor(C), "1908")1858 -2.83415    0.37894   -7.479 7.48e-14 ***
relevel(factor(C), "1908")1863 -2.82369    0.18322  -15.412  < 2e-16 ***
relevel(factor(C), "1908")1868 -2.32070    0.08718  -26.619  < 2e-16 ***
relevel(factor(C), "1908")1873 -2.00731    0.05296  -37.902  < 2e-16 ***
relevel(factor(C), "1908")1878 -1.78079    0.03780  -47.110  < 2e-16 ***
relevel(factor(C), "1908")1883 -1.41511    0.02813  -50.307  < 2e-16 ***
relevel(factor(C), "1908")1888 -0.92218    0.02115  -43.595  < 2e-16 ***
relevel(factor(C), "1908")1893 -0.55342    0.01751  -31.606  < 2e-16 ***
relevel(factor(C), "1908")1898 -0.27621    0.01537  -17.968  < 2e-16 ***
relevel(factor(C), "1908")1903 -0.07100    0.01412   -5.027 4.99e-07 ***
relevel(factor(C), "1908")1913 -0.02005    0.01372   -1.461 0.144019    
relevel(factor(C), "1908")1918  0.02813    0.01419    1.983 0.047369 *  
relevel(factor(C), "1908")1923  0.08136    0.01499    5.427 5.72e-08 ***
relevel(factor(C), "1908")1928  0.03131    0.01751    1.788 0.073699 .  
relevel(factor(C), "1908")1933  0.02899    0.02128    1.363 0.173032    
relevel(factor(C), "1908")1938 -0.11124    0.02839   -3.919 8.91e-05 ***
relevel(factor(C), "1908")1943 -0.14656    0.03794   -3.863 0.000112 ***
relevel(factor(C), "1908")1948  0.01684    0.05424    0.311 0.756136    
relevel(factor(C), "1908")1953 -0.02004    0.11049   -0.181 0.856109    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 1.0037e+08  on 110  degrees of freedom
Residual deviance: 8.2963e+02  on  81  degrees of freedom
AIC: 1744.7

Number of Fisher Scoring iterations: 4

ACモデルの年齢効果をグラフに表す

分析結果から、 推定値と標準誤差を取り出し、 年齢区分の横軸を作り、 黒実線で推定値、 赤細線で95%信頼区間を描画する。

ac.cf <- summary(ac.3)$coef
age.ac.cf <- ac.cf[1:10,1:2]
am <- seq(40,85,5)+2.5
plot(am,exp(age.ac.cf[,1]),type="l",log="y",lwd=3, ylab="Lung cancer mortality rate", xlab="Age")
lines(am,exp(age.ac.cf[,1]-1.96*age.ac.cf[,2]),col="red")
lines(am,exp(age.ac.cf[,1]+1.96*age.ac.cf[,2]),col="red")

ACモデルのコホート効果をグラフに表す

解析結果から、 コホートの推定値と標準誤差を取り出し、 参照カテゴリ"1908"に推定値0、標準誤差0を挿入する。

グラフが上下対称になるように 0.05(20分の1)から20.0(20倍)にY軸を指定。

cht.ac.cf <- ac.cf[11:29,1:2]
cht.ac.cf <- rbind(cht.ac.cf[1:10,],c(0,0),cht.ac.cf[11:19,])
cm <- seq(1858,1953,5)
plot(cm,exp(cht.ac.cf[,1]),type="l",log="y",lwd=3,ylim=c(0.05,20.0), ylab="Relative risk of lung cancer death for cohort", xlab="Cohort")
lines(cm, exp(cht.ac.cf[,1]-1.96*cht.ac.cf[,2]),col="red")
lines(cm, exp(cht.ac.cf[,1]+1.96*cht.ac.cf[,2]),col="red")

両端ともデータが少ないことを物語っていて、 95%信頼区間が広くなっている。

統計ソフトRでAPCモデル

ポアソン回帰モデルでAPCモデルを分析する

A, P, Cすべてを因子型でFormulaに投入する。 Formulaに-1を追加して、 切片なしモデルで年齢区分の推定値を死亡率の絶対値にする。 また、時代区分、コホート区分の参照カテゴリをそれぞれ、 1968と1908の区分にする。

apc.3 <- glm (D~ -1 + factor(A) + relevel(factor(P), "1968") + relevel(factor(C), "1908") + offset(log(Y)), family=poisson, data=dat)
summary(apc.3)

結果は、統計学的有意でないカテゴリが散見される。 時代区分、コホート区分ともに、 参照カテゴリより現代になると統計学的有意でなく、 参照カテゴリと比較して変化しないと言える。

> summary(apc.3)

Call:
glm(formula = D ~ -1 + factor(A) + relevel(factor(P), "1968") + 
    relevel(factor(C), "1908") + offset(log(Y)), family = poisson, 
    data = dat)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-4.0478  -0.8303   0.0357   0.8045   3.1655  

Coefficients: (1 not defined because of singularities)
                               Estimate Std. Error  z value Pr(>|z|)    
factor(A)40                    -8.88084    0.05166 -171.899  < 2e-16 ***
factor(A)45                    -7.92087    0.04360 -181.663  < 2e-16 ***
factor(A)50                    -7.07551    0.03083 -229.476  < 2e-16 ***
factor(A)55                    -6.42381    0.02071 -310.198  < 2e-16 ***
factor(A)60                    -5.93048    0.01600 -370.693  < 2e-16 ***
factor(A)65                    -5.57284    0.02036 -273.723  < 2e-16 ***
factor(A)70                    -5.32555    0.02989 -178.145  < 2e-16 ***
factor(A)75                    -5.23156    0.04134 -126.559  < 2e-16 ***
factor(A)80                    -5.25056    0.05391  -97.398  < 2e-16 ***
factor(A)85                    -5.38676    0.06900  -78.070  < 2e-16 ***
relevel(factor(P), "1968")1943 -0.48206    0.07035   -6.852 7.29e-12 ***
relevel(factor(P), "1968")1948 -0.35244    0.05609   -6.283 3.31e-10 ***
relevel(factor(P), "1968")1953 -0.30889    0.04328   -7.137 9.52e-13 ***
relevel(factor(P), "1968")1958 -0.17921    0.03109   -5.765 8.16e-09 ***
relevel(factor(P), "1968")1963 -0.09615    0.02083   -4.616 3.91e-06 ***
relevel(factor(P), "1968")1973  0.01905    0.01940    0.982   0.3261    
relevel(factor(P), "1968")1978  0.05178    0.02877    1.800   0.0718 .  
relevel(factor(P), "1968")1983  0.04057    0.03987    1.017   0.3090    
relevel(factor(P), "1968")1988 -0.07113    0.05157   -1.379   0.1678    
relevel(factor(P), "1968")1993 -0.14006    0.06237   -2.246   0.0247 *  
relevel(factor(C), "1908")1858 -2.29805    0.40016   -5.743 9.31e-09 ***
relevel(factor(C), "1908")1863 -2.33887    0.21621  -10.817  < 2e-16 ***
relevel(factor(C), "1908")1868 -1.87612    0.13402  -13.999  < 2e-16 ***
relevel(factor(C), "1908")1873 -1.61119    0.10348  -15.569  < 2e-16 ***
relevel(factor(C), "1908")1878 -1.44007    0.08500  -16.941  < 2e-16 ***
relevel(factor(C), "1908")1883 -1.13903    0.06939  -16.415  < 2e-16 ***
relevel(factor(C), "1908")1888 -0.71654    0.05502  -13.022  < 2e-16 ***
relevel(factor(C), "1908")1893 -0.41829    0.04203   -9.952  < 2e-16 ***
relevel(factor(C), "1908")1898 -0.20382    0.02993   -6.810 9.73e-12 ***
relevel(factor(C), "1908")1903 -0.04489    0.01935   -2.320   0.0203 *  
relevel(factor(C), "1908")1913 -0.02874    0.01791   -1.605   0.1086    
relevel(factor(C), "1908")1918  0.02011    0.02773    0.725   0.4683    
relevel(factor(C), "1908")1923  0.07705    0.03907    1.972   0.0486 *  
relevel(factor(C), "1908")1928  0.02919    0.05134    0.569   0.5696    
relevel(factor(C), "1908")1933  0.02670    0.06399    0.417   0.6765    
relevel(factor(C), "1908")1938 -0.11619    0.07759   -1.498   0.1343    
relevel(factor(C), "1908")1943 -0.15336    0.09175   -1.672   0.0946 .  
relevel(factor(C), "1908")1948  0.01978    0.10815    0.183   0.8549    
relevel(factor(C), "1908")1953       NA         NA       NA       NA    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 1.0037e+08  on 110  degrees of freedom
Residual deviance: 2.0855e+02  on  72  degrees of freedom
AIC: 1141.6

Number of Fisher Scoring iterations: 4

年齢・時代・コホート効果をそれぞれグラフに表す

まず、推定値と標準誤差を取り出す。

apc.cf <- summary(apc.3)$coef

年齢効果をグラフに表す。 高年齢区分で信頼区間が広くなっている。

age.apc.cf <- apc.cf[1:10,1:2]
am <- seq(40,85,5)+2.5
plot(am, exp(age.apc.cf[,1]), type="l", log="y", lwd=3, ylab="Lung cancer mortality rate", xlab="Age")
lines(am,exp(age.apc.cf[,1]-1.96*age.apc.cf[,2]),col="red")
lines(am,exp(age.apc.cf[,1]+1.96*age.apc.cf[,2]),col="red")

時代効果をグラフに表す。 参照カテゴリ"1968"に推定値0、標準誤差0を挿入。 上下対称になるように、0.2(5分の1)から5.0(5倍)にY軸を指定。

per.apc.cf <- apc.cf[11:20,1:2]
per.apc.cf <- rbind(per.apc.cf[1:5,],c(0,0),per.apc.cf[6:10,])
pm <- seq(1943,1993,5)+2.5
plot(pm, exp(per.apc.cf[,1]),type="l", log="y", lwd=3, ylim=c(0.2,5.0),ylab="Relative risk of lung cancer death for period",xlab="Period")
lines(pm, exp(per.apc.cf[,1]-1.96*per.apc.cf[,2]),col="red")
lines(pm, exp(per.apc.cf[,1]+1.96*per.apc.cf[,2]),col="red")

コホート効果をグラフに表す。 参照カテゴリ"1908"に推定値0、標準誤差0を挿入。 上下対称になるように、0.05(20分の1)から20.0(20倍)にY軸を指定。

cht.apc.cf <- apc.cf[21:38,1:2]
cht.apc.cf <- rbind(cht.apc.cf[1:10,],c(0,0),cht.apc.cf[11:18,])
cm <- seq(1858,1948,5)
plot(cm,exp(cht.apc.cf[,1]),type="l",log="y",lwd=3,ylim=c(0.05,20.0),ylab="Relative risk of lung cancer death for cohort",xlab="Cohort")
lines(cm, exp(cht.apc.cf[,1]-1.96*cht.apc.cf[,2]),col="red")
lines(cm, exp(cht.apc.cf[,1]+1.96*cht.apc.cf[,2]),col="red")

以上が、統計ソフトRでのAPCモデル解析結果グラフ図示。

統計ソフトRのライブラリでAPCモデルはできる?

ライブラリでAPCモデルを簡単にしかもきれいなグラフ描画が できないものだろうか?

実はEpiというライブラリで実現可能だ。

初回はEpiをインストールして、

install.packages("Epi")

Epiを呼び出し、 apc.fit()でポアソン回帰モデルを適用する。 apc.plot()で結果のグラフ表示をしてくれる。

library(Epi)
apc1 <- apc.fit (A=A, P=P, D=D, Y=Y, data=dat)
apc.plot(apc1)

print()で、計算結果が表示できる。

> print(apc1)
$`Type`
[1] "ML of APC-model Poisson with log(Y) offset : ( ACP ):\n"

$Model

Call:  glm(formula = D ~ MA + MPr + MCr - 1, family = poisson, offset = log(Y))

Coefficients:
      MA       MA1       MA2       MA3       MA4      MPr1      MPr2      MPr3       MCr      MCr1      MCr2  
-7.36188   1.60570   1.64203   3.15021   1.26784   0.28354   0.19235   0.36059   0.01963   0.25778   0.20859  
    MCr3  
 0.76182  

Degrees of Freedom: 110 Total (i.e. Null);  98 Residual
Null Deviance:      100400000 
Residual Deviance: 419.3        AIC: 1300

$Age
      Age         Rate         2.5%        97.5%
 [1,]  40 0.0001320909 0.0001269518 0.0001374379
 [2,]  45 0.0002896177 0.0002820562 0.0002973820
 [3,]  50 0.0006350055 0.0006248035 0.0006453740
 [4,]  55 0.0013494813 0.0013299702 0.0013692788
 [5,]  60 0.0023778182 0.0023486272 0.0024073720
 [6,]  65 0.0032649683 0.0032147255 0.0033159963
 [7,]  70 0.0044643607 0.0043994078 0.0045302727
 [8,]  75 0.0047922282 0.0047260472 0.0048593359
 [9,]  80 0.0048039824 0.0047052431 0.0049047938
[10,]  85 0.0048157655 0.0046553485 0.0049817103

$Per
       Per      P-RR      2.5%     97.5%
 [1,] 1943 0.7594381 0.7364244 0.7831710
 [2,] 1948 0.8188221 0.8020877 0.8359057
 [3,] 1953 0.8828497 0.8730873 0.8927212
 [4,] 1958 0.9518838 0.9459791 0.9578253
 [5,] 1963 1.0233283 1.0114140 1.0353829
 [6,] 1968 1.0810581 1.0675224 1.0947655
 [7,] 1973 1.1062089 1.0934462 1.1191205
 [8,] 1978 1.0942631 1.0799575 1.1087583
 [9,] 1983 1.0434487 1.0310683 1.0559778
[10,] 1988 0.9493130 0.9408789 0.9578228
[11,] 1993 0.8400780 0.8250205 0.8554104

$Coh
       Coh       C-RR       2.5%      97.5%
 [1,] 1858 0.06757542 0.06271255 0.07281537
 [2,] 1863 0.09278350 0.08696768 0.09898825
 [3,] 1868 0.12739512 0.12059620 0.13457735
 [4,] 1873 0.17491813 0.16721045 0.18298110
 [5,] 1878 0.24016895 0.23179774 0.24884248
 [6,] 1883 0.32976069 0.32120232 0.33854710
 [7,] 1888 0.45277341 0.44463086 0.46106507
 [8,] 1893 0.62167433 0.61356987 0.62988585
 [9,] 1898 0.83878971 0.82728495 0.85045447
[10,] 1903 1.01904189 1.00603099 1.03222106
[11,] 1908 1.05319722 1.04120907 1.06532339
[12,] 1913 1.06037791 1.04630862 1.07463637
[13,] 1918 1.16117323 1.14846175 1.17402540
[14,] 1923 1.21151927 1.19583133 1.22741302
[15,] 1928 1.20272146 1.18831276 1.21730486
[16,] 1933 1.18413176 1.16303925 1.20560679
[17,] 1938 1.16582939 1.13444565 1.19808134
[18,] 1943 1.14780991 1.10550646 1.19173215
[19,] 1948 1.13006894 1.07693312 1.18582649
[20,] 1953 1.11260219 1.04893119 1.18013806

$Drift
                exp(Est.)     2.5%    97.5%
APC (D-weights)  1.019828 1.019231 1.020426
A-d              1.023573 1.023057 1.024088

$Ref
 Per  Coh 
  NA 1913 

$Anova

Analysis of deviance for Age-Period-Cohort model

                  Resid. Df Resid. Dev Df Deviance  Pr(>Chi)    
Age                     105    15242.0                          
Age-drift               104     6564.0  1   8678.0 < 2.2e-16 ***
Age-Cohort              101     1016.4  3   5547.6 < 2.2e-16 ***
Age-Period-Cohort        98      419.3  3    597.1 < 2.2e-16 ***
Age-Period              101     2910.5 -3  -2491.3 < 2.2e-16 ***
Age-drift               104     6564.0 -3  -3653.5 < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

$Knots
$Knots$`Age`
[1] 50 60 65 70 75

$Knots$Per
[1] 1958 1968 1978 1983 1993

$Knots$Coh
[1] 1893 1903 1913 1918 1928


attr(,"class")
[1] "apc"

 

まとめ

統計ソフトRでAPCモデルを分析してみた。

glm()を使って解析し、 plot()を使って図示することはできるが、 ライブラリEpiを使えば、 あっというまに解析出来て、 きれいなグラフも描ける。

Epiはとても便利なライブラリだ。