檢討與建議

國

立 政 治 大

㈻㊫學

•

N a tio na

l C h engchi U ni ve rs it y

第五章 結論與建議

第一節 結論

由前述各章節的說明與結果，我們總結如下：

一、在男女自殺死亡率的比較方面，就年齡方面來看，女性的自殺死亡率在 10 歲到 24 歲時顯著比男性高，特別在 15 到 19 歲時差異達到最大，20 歲之後 差異開始變小，到了 25 至 34 歲，兩性已無顯著差異，但是 35 歲後男性開 始顯著大於女性，並且隨著年齡增長兩性的差異越大，直到 60 歲之後差異 才開始變小，到 70 歲時兩性已無顯著差異。從年代方面來看，男女的自殺 死亡率在 1959 年到 1973 年間並無顯著的差異，1974 年開始到 2008 年間呈 現一個 J 型曲線狀，在 1974 到 1988 年女性的自殺死亡率顯著大於男性並 於 1979 年到 1983 年來到最低點，也就是差異最大，之後差異開始變小，

於 1989 年時兩性已無顯著差異，直到 1994 年開始男性的自殺死亡率反而 開始顯著大於女性，而且隨著年代越近差異越大，在 2008 年差異達到最 大。

二、在參數估計方法方面，非條件概似函數法與條件概似函數法兩者估計出來 的參數相去不遠 ，但前者的過度離散度的確比後者的大很多。

第二節 檢討與建議

當資料有分層時，在做年齡-年代-世代分析時，以往文獻多半直接對各分

層資料藉由繪圖的方式來進行比較。這些方式所得出的結果可能會過於主觀，

有鑑於此，本文採用 Held and Riebler (2010)的多元年齡-年代-世代模型來對男 女自殺死亡率作比較。我們同時採用條件概似函數法與非條件概似函數法來配 適多元年齡-年代-世代模型，並比較兩種方法之差異。

然而，在配適多元年齡-年代-世代模型時須設定一致效應，因此得事先瞭 解資料在哪個效應上男女的趨勢是較一致的，如果各層資料在三種效應的趨勢

都不一致，或是取錯一致效應則有可能導致錯的估計，這也是使用 Held and

Riebler (2010)所建議之方法的最大限制。

男女在三種效應上的差異趨勢，其背後所隱藏的原因，也是值得深入研究

• 國

立 政 治 大

㈻㊫學

•

N a tio na

l C h engchi U ni ve rs it y

與探討的 。像是為甚麼男女自殺死亡的相對風險在年代方面，從 1974 年到 2008 年呈現一個 J 型曲線狀，在年齡方面為什麼在 24 歲以下女性顯著大於男 性，35 歲之後男性則開始顯著大於女性，這些現象都需要社會學家，精神科醫 師等專業人士，藉由實際的臨床經驗及知識來進行解釋。

• 國

立 政 治 大

㈻㊫學

•

N a tio na

l C h engchi U ni ve rs it y

參考文獻

林佳瑩，蔡毓智 (2005)：台灣地區自殺趨勢研究：1976-2001。家醫研究，

3(1) ，28-38

Ananth,C.V., Misra,D.P., Demissie,K., and Smulian,J.C. (2001),“Rates of Preterm Delivery among Black Women and White Women in the United States over Two Decades: An Age-Period-Cohort Analysis”, American of Journal

Epidemiology 154(7): 657–665.

Burnham,K.P., and Anderson,D.R. (2002), Model Selection and Multivariate

Inference: A Practical Information-Theoretic Approach, 2nd edn., New York:

Springer

Fienberg,S.E., and Mason,W.M.（1978）,“Identification and Estimation of Age-Period-Cohort Models in the Analysis of Discrete Archival Data”, Sociological

Methodology 8: 1–67.

Frost,W.H. (1939),“The Age Selection of Mortality from Tuberculosis in Successive Decades”, Amreican Joural of Hygiene 30: 92–96.

Fu,W.J.（2000）,“Ridge Estimator in Singular Design with Application to Age-Period-Cohort Analysis of Disease Rates”, Communications in Statistics–Theory

and Method 29: 263–278.

Glenn,N.D. (2005),“Age, Period, and Cohort Effects”, Encyclopedia of Social

Measurement 1: 27–32.

Granizo,J.J., Guallar,E., and Artalejo,F.R. (1996),“Age-Period-Cohort Analysis of Suicide Mortality Rates in Spain, 1959-1991”,International Journal of

Epidemiology 25(4): 814–820.

• 國

立 政 治 大

㈻㊫學

•

N a tio na

l C h engchi U ni ve rs it y

Gross,V.A., Bopp,M., Gostynski,M., Lauber,C., Gutzwiller,F., and Rössler,W.

(2006),”Age-Period-Cohort Analysis of Swiss Suicide Data, 1881–2000”, Eur

Arch Psychiatry Clin Neurosci 256 : 207–214

Gunnell,D., Middleton,N., Whitley,E., Dorling,D., and Frankel,S. (2003), ”Influence of Cohort Effects on Patterns of Suicide in England and

Wales,1950-1999”, British Journal of Psychiatry 182:164-170

Held,L., and Riebler,A. (2010),“A Conditional Approach for Inference in Multivariate Age-Period-Cohort Models”, Statistical Methods in Medical

Research 0(0): 1–19.

Kupper,L.L., Janis,J.M., Salama,I.A., Yoshizawa,C.N., and Greenberg,B.G.

（1983）,“Age-Period-Cohort Analysis: An Illustration of the Problems in Assessing Interaction in One Observation Per Cell Data”, Communications in

Statistics–Theory and Method 12: 2779–2807.

Riebler,A., and Held,L. (2010),“The Analysis of Heterogeneous Time Trends in Multivariate Age-Period-Cohort Models”, Biostatistics 11: 57-69

Wu,R., and Cheng,Y. (2008),“Trends in Suicide Mortality in Taiwan, 1959-2006”, Taiwan J Public Health 27(2): 110–120.

Yang,Y. (2008),“Trends in U.S. Adult Chronic Disease Mortality, 1960-1999:

Age, Period, and Cohort variations”, Demography 45(2): 387–416.

Yang,Y., Fu,W.J., and Land,K.C.(2004),“A Methodological Comparison of Age-Period-Cohort Models: The Intrinsic Estimator and Conventional Generalized Linear Models”, Sociological Methodology 34: 75–110.

Yang,Y., Fu, W.J., Land,K.C., and Wohl,S.S. (2008),“The Intrinsic Estimator for Age-Period-Cohort Analysis: What It Is and How to Use It”, American Journal

of Sociology 113(6): 1697–1736.

•

1959-1963 1974-1978 1989-1993 2004-2008

Apc

Period

average log relative risk

-3

1875-1883 1880-1888 1885-1893 1890-1898 1895-1903 1900-1908 1905-1913 1910-1918 1915-1923 1920-1927 1925-1993 1930-1938 1935-1943 1940-1948 1945-1953 1950-1958 1955-1963 1960-1968 1965-1973 1970-1978 1975-1983 1980-1988 1985-1993 1990-1998

Apc

average log relative risk

-2

10-14 25-29 40-44 55-59 70-74 80+

aPc

Age

average log relative risk

-2

1875-1883 1880-1888 1885-1893 1890-1898 1895-1903 1900-1908 1905-1913 1910-1918 1915-1923 1920-1927 1925-1993 1930-1938 1935-1943 1940-1948 1945-1953 1950-1958 1955-1963 1960-1968 1965-1973 1970-1978 1975-1983 1980-1988 1985-1993 1990-1998

aPc

average log relative risk

-0.4 -0.20.00.20.40.60.8

10-14 25-29 40-44 55-59 70-74 80+

apC

Age

average log relative risk

-0.4 -0.20.00.20.40.60.8

1959-1963 1974-1978 1989-1993 2004-2008

apC

Period

average log relative risk

表6-1 平均相對風險-非條件概似函數法

•

1959-1963 1974-1978 1989-1993 2004-2008 Apc

Period

average log relative risk

-2 -1 0 1

1875-1883 1880-1888 1885-1893 1890-1898 1895-1903 1900-1908 1905-1913 1910-1918 1915-1923 1920-1927 1925-1993 1930-1938 1935-1943 1940-1948 1945-1953 1950-1958 1955-1963 1960-1968 1965-1973 1970-1978 1975-1983 1980-1988 1985-1993 1990-1998

Apc

average log relative risk

-1.0

10-14 25-29 40-44 55-59 70-74 80+

aPc

Age

average log relative risk

-1.0

1875-1883 1880-1888 1885-1893 1890-1898 1895-1903 1900-1908 1905-1913 1910-1918 1915-1923 1920-1927 1925-1993 1930-1938 1935-1943 1940-1948 1945-1953 1950-1958 1955-1963 1960-1968 1965-1973 1970-1978 1975-1983 1980-1988 1985-1993 1990-1998

aPc

average log relative risk

-0.2

10-14 25-29 40-44 55-59 70-74 80+

apC

Age

average log relative risk

-0.2

1959-1963 1974-1978 1989-1993 2004-2008 apC

Period

average log relative risk

表6-2 平均相對風險-條件概似函數法

•

########## unconditional approach##############

#########mapc-aPc#####################

J <- ncol(data11) # number of age groups I <- nrow(data11) # number of periods K <- J+I-1 # number of cohorts

G <- 1 # grid factor (here agegroups and periods are provided in 5-year

intervals)

P <- rep(rep(1:J, I),2) # period index A <- rep(rep(1:I, each=J),2) # age index C <- P-A+I # cohort index

R <- rep(1:2, each=I*J) # strata index

########x matrix##########

x=matrix(0,ncol=(I-1)*2+J-1+(K-1)*2,nrow=2*I*J) for(i in 1:(I*J)){if(A[i]!=I)x[i,A[i]]=1

else{if(A[i]==I)x[i,c(1:(I-1))]=-1}}

for(i in ((I*J)+1):(2*I*J)){if(A[i]!=I)x[i,I-1+A[i]]=1

else{if(A[i]==I)x[i,c(I:((I-1)*2))]=-1}}

for(i in 1:(I*J*2)){if(P[i]!=J)x[i,((I-1)*2+P[i])]=1 else{if(P[i]==J)x[i,c(((I-1)*2+1):((I-1)*2+J-1))]=-1}}

for(i in 1:(I*J)){if(C[i]!=K)x[i,((I-1)*2+J-1+C[i])]=1

else{if(C[i]==K)x[i,c(((I-1)*2+J):((I-1)*2+J-1+K-1))]=-1}}

for(i in ((I*J)+1):(2*I*J)){if(C[i]!=K)x[i,((I-1)*2+J-1+K-1+C[i])]=1

else{if(C[i]==K)x[i,c(((I-1)*2+J-1+K):((I-1)*2+J-1+(K-1)*2))]=-1}}

m1=c(rep(1,(I*J)),rep(0,(I*J))) m2=c(rep(0,(I*J)),rep(1,(I*J))) x=cbind(m1,m2,x)

N=rbind(data12,data22) #data of population Y=rbind(data11,data21) #data of suicide death

#######將data轉為vector#############

y=c()

for(i in 1:(I*2)) {y=c(y,Y[i,])}

y=as.numeric(y) n=c()

for(i in 1:(I*2)) {n=c(n,N[i,])}

n=as.numeric(n)

######設定限制式##########

############e1=e2###############

for(i in 1:(2*I*J)){

if(x[i,(32+2)]+x[i,(33+2)]==1)x[i,(32+2)]=1 else{if(x[i,(32+2)]+x[i,(33+2)]==-2)x[i,(32+2)]=-2 else{x[i,(32+2)]=0}}}

x=x[,-(33+2)]

######fit model########

lm1=glm.fit(x,y,family=poisson(),offset=log(n),intercept=F)

•

coe=as.matrix(lm1$coefficients)#model coefficients summary=summary.glm(lm1)

cov=summary$cov.unscaled #unscaled covariance

intercept=coe[1:2]

coe=coe[-c(1:2)]

coe=c(coe[1:32],coe[32],coe[(1+32):(2*(I-1)+J-1+(K-1)*2-1)]) coea1=coe[1:(I-1)]

coea2=coe[I:(2*(I-1))]

coep=coe[(2*(I-1)+1):(2*(I-1)+J-1)]

coec1=coe[(2*(I-1)+J):(2*(I-1)+J-1+K-1)]

coec2=coe[(2*(I-1)+J-1+K):(2*(I-1)+J-1+2*(K-1))]

coe1=c(intercept[1],coea1,0-sum(coea1),coep,0-sum(coep),coec1,0-sum(coec1)) coe2=c(intercept[2],coea2,0-sum(coea2),coep,0-sum(coep),coec2,0-sum(coec2)) #######參數差#######

result=coe1-coe2

resulta=c(result[2:16])#年齡估計參數 resultc=c(result[27:50])#世代的估計參數

cov=summary.glm(lm1,correlation = T)$cov.unscaled

cov=(3529/df.residual(lm1))*cov#乘上過度離散度(scaled covariance)

#########年齡部份的標準差##########

sta=c()

for(i in

1:(I-1)){sta=c(sta,c(0,0,rep(0,i-1),1,rep(0,I-2),-1,rep(0,(I-1)-i+2*(K-1)+(J-2)))%*%cov%*%t(t(c(0,0,rep(0,i-1),1,rep(0,I-2),-1,rep(0,(I-1)-i+2*(K-1)+(J-2))))))}

lasta=c(0,0,rep(-1,I-1),rep(1,I-1),rep(0,(J-2)+2*(K-1)))%*%cov%*%t(t(c(0,0,rep(-1,I-1),rep(1,I-1),rep(0,(J-2)+2*(K-1)))))

sta=c(sqrt(sta),sqrt(lasta))

#########世代部份的標準差##########

stc=c()

for(i in 1:(K-1)){stc=c(stc,c(0,0,rep(0,2*(I-1)+(J-2)+i-1),1,rep(0,K-2),-1,rep(0,K-1-i))%*%cov%*%t(t(c(0,0,rep(0,2*(I-1)+(J-2)+i-1),1,rep(0,K-2),-1,rep(0,K-1-i)))))}

lastc=c(0,0,rep(0,2*(I-1)+(J-2)),rep(-1,K-1),rep(1,K-1))%*%cov%*%t(t(c(0,0,rep(0,2*(I-1)+(J-2)),rep(-1,K-1),rep(1,K-1)))) stc=c(sqrt(stc),sqrt(lastc))

quant <- -qnorm((1-0.95)/2)#

#####plot######

######age##########

plot(resulta,type="n",axes=FALSE,ann=FALSE,ylim=c(-3,3)) lines(resulta,col=3,lwd=2)

axis(2, col.axis="blue", las=1) axis(1, at=1:15, lab=c("10-14",

"15-19",

•

"80+"), col.axis="blue")

title(main="aPc ", font.main=4, col.main="red") title(xlab="Age", col.lab="red")

title(ylab="difference",col.lab="red")

axis(2, col.axis="blue", las=1) axis(1, at=1:24, lab=c("1875-1883",

"1880-1888",

"1990-1998"), col.axis="blue", las=2)

title(main="aPc ", font.main=4, col.main="red")

•

for(i in 1:(I-1)){sta=c(sta,c(1,-1,rep(0,i-1),1,rep(0,I-2),-1,rep(0,(I-1)-i+2*(K-1)+(J-2)))%*%cov%*%t(t(c(1,-1,rep(0,i-1),1,rep(0,I-2),-1,rep(0,(I-1)-i+2*(K-1)+(J-2))))))}

lasta=c(1,-1,rep(-1,I-1),rep(1,I-1),rep(0,(J-2)+2*(K-1)))%*%cov%*%t(t(c(1,-1,rep(-1,I-1),rep(1,I-1),rep(0,(J-2)+2*(K-1)))))

sta=c(sqrt(sta),sqrt(lasta)) stc=c()

for(i in 1:(K-1)){stc=c(stc,c(1,-1,rep(0,2*(I-1)+(J-2)+i-1),1,rep(0,K-2),-1,rep(0,K-1-i))%*%cov%*%t(t(c(1,-1,rep(0,2*(I-1)+(J-2)+i-1),1,rep(0,K-2),-1,rep(0,K-1-i)))))}

lastc=c(1,-1,rep(0,2*(I-1)+(J-2)),rep(-1,K-1),rep(1,K-1))%*%cov%*%t(t(c(1,-1,rep(0,2*(I-1)+(J-2)),rep(-1,K-1),rep(1,K-1))))

stc=c(sqrt(stc),sqrt(lastc))

result=coe1-coe2

axis(2, col.axis="blue", las=1) axis(1, at=1:15, lab=c("10-14",

"15-19",

"80+"), col.axis="blue")

title(main="aPc ", font.main=4, col.main="red") title(xlab="Age", col.lab="red")

title(ylab="average log relative risk",col.lab="red") lines(resulta+quant*sta,lty=2,col=3,lwd=2)

•

axis(2, col.axis="blue", las=1) axis(1, at=1:24, lab=c("1875-1883",

"1880-1888",

"1990-1998"), col.axis="blue", las=2)

title(main="aPc ", font.main=4, col.main="red")

title(ylab="average log relative risk",col.lab="red") lines(resultc+quant*stc,lty=2,col=3,lwd=2)

lines(resultc-quant*stc,lty=2,col=3,lwd=2) abline(0,0)

#########mapc-Apc#####################

x=matrix(0,ncol=I-1+(J-1)*2+(K-1)*2,nrow=2*I*J) for(i in 1:(I*J*2)){if(A[i]!=I)x[i,A[i]]=1

else{if(A[i]==I)x[i,c(1:(I-1))]=-1}}

for(i in 1:(I*J)){if(P[i]!=J)x[i,((I-1)+P[i])]=1 else{if(P[i]==J)x[i,c(((I-1)+1):((I-1)+J-1))]=-1}}

for(i in ((I*J)+1):(2*I*J)){if(P[i]!=J)x[i,I-1+J-1+P[i]]=1 else{if(P[i]==J)x[i,c((I-1+J):((I-1)+2*(J-1)))]=-1}}

for(i in 1:(I*J)){if(C[i]!=K)x[i,((I-1)+(J-1)*2+C[i])]=1

else{if(C[i]==K)x[i,c((I+(J-1)*2):((I-1)+(J-1)*2+K-1))]=-1}}

for(i in ((I*J)+1):(2*I*J)){if(C[i]!=K)x[i,((I-1)+(J-1)*2+K-1+C[i])]=1 else{if(C[i]==K)x[i,c(((I-1)+(J-1)*2+K):((I-1)+(J-1)*2+(K-1)*2))]=-1}}

m1=c(rep(1,(I*J)),rep(0,(I*J)))

•

m2=c(rep(0,(I*J)),rep(1,(I*J))) x=cbind(m1,m2,x)

N=rbind(data12,data22) Y=rbind(data11,data21) y=c()

for(i in 1:(I*2)) {y=c(y,Y[i,])}

y=as.numeric(y) n=c()

for(i in 1:(I*2)) {n=c(n,N[i,])}

n=as.numeric(n)

############e1=e2###############

for(i in 1:(2*I*J)){

if(x[i,(1+2)]+x[i,(2+2)]==1)x[i,(1+2)]=1 else{if(x[i,(1+2)]+x[i,(2+2)]==-2)x[i,(1+2)]=-2 else{x[i,(1+2)]=0}}}

x=x[,-(2+2)]

lm1=glm.fit(x,y,family=poisson(),offset=log(n),intercept=F) coe=as.matrix(lm1$coefficients)

summary=summary.glm(lm1) cov1=summary$cov.unscaled intercept=coe[1:2]

coe=coe[-c(1:2)]

coe=c(coe[1:1],coe[1],coe[(1+1):((I-1)+(J-1)*2+(K-1)*2-1)]) coea=coe[1:(I-1)]

coep1=coe[((I-1)+1):((I-1)+J-1)]

coep2=coe[((I-1)+J):((I-1)+2*(J-1))]

coec1=coe[((I-1)+(J-1)*2+1):((I-1)+2*(J-1)+(K-1))]

coec2=coe[((I-1)+2*(J-1)+K):((I-1)+2*(J-1)+2*(K-1))]

coe1=c(intercept[1],coea,0-sum(coea),coep1,0-sum(coep1),coec1,0-sum(coec1)) coe2=c(intercept[2],coea,0-sum(coea),coep2,0-sum(coep2),coec2,0-sum(coec2)) result=coe1-coe2

resultp=c(result[17:26]) resultc=c(result[27:50])

cov=summary.glm(lm1,correlation = T)$cov.unscaled cov=(4099/df.residual(lm1))*cov

stp=c()

for(i in 1:(J-1)){stp=c(stp,c(0,0,rep(0,I-2+i-1),1,rep(0,J-2),-1,rep(0,2*(K-1)+J-1-i))%*%cov%*%t(t(c(0,0,rep(0,I-2+i-1),1,rep(0,J-2),-1,rep(0,2*(K-1)+J-1-i)))))}

lastp=c(0,0,rep(0,I-2),rep(-1,J-1),rep(1,J-1),rep(0,2*(K-1)))%*%cov%*%t(t(c(0,0,rep(0,I-2),rep(-1,J-1),rep(1,J-1),rep(0,2*(K-1))))) stp=c(sqrt(stp),sqrt(lastp))

stc=c()

for(i in 1:(K-1)){stc=c(stc,c(0,0,rep(0,(I-2)+(2*(J-1))+i-1),1,rep(0,K-2),-1,rep(0,K-1-i))%*%cov%*%t(t(c(0,0,rep(0,(I-2)+(2*(J-1))+i-1),1,rep(0,K-2),-1,rep(0,K-1-i)))))}

•

stc=c(sqrt(stc),sqrt(lastc))

quant <- -qnorm((1-0.95)/2) #####plot######

par(mfrow=c(3,2))

plot(resultp,type="n",axes=FALSE,ann=FALSE,ylim=c(-3,3)) lines(resultp,col=4,lwd=2)

axis(2, col.axis="blue", las=1) axis(1, at=1:10, lab=c("1959-1963",

"1964-1968",

"2004-2008"), col.axis="blue")

title(main="Apc ", font.main=4, col.main="red") title(xlab="Period", col.lab="red")

title(ylab="difference",col.lab="red")

axis(2, col.axis="blue", las=1) axis(1, at=1:24, lab=c("1875-1883",

"1880-1888",

•

"1990-1998"), col.axis="blue", las=2)

title(main="Apc ", font.main=4, col.main="red") title(ylab="difference",col.lab="red")

cov=summary.glm(lm1,correlation = T)$cov.unscaled cov=(4099/df.residual(lm1))*cov

stp=c()

for(i in 1:(J-1)){stp=c(stp,c(1,-1,rep(0,I-2+i-1),1,rep(0,J-2),-1,rep(0,2*(K-1)+J-1-i))%*%cov%*%t(t(c(1,-1,rep(0,I-2+i-1),1,rep(0,J-2),-1,rep(0,2*(K-1)+J-1-i)))))}

lastp=c(1,-1,rep(0,I-2),rep(-1,J-1),rep(1,J-1),rep(0,2*(K-1)))%*%cov%*%t(t(c(1,-1,rep(0,I-2),rep(-1,J-1),rep(1,J-1),rep(0,2*(K-1)))))

stp=c(sqrt(stp),sqrt(lastp)) stc=c()

for(i in 1:(K-1)){stc=c(stc,c(1,-1,rep(0,(I-2)+(2*(J-1))+i-1),1,rep(0,K-2),-1,rep(0,K- 1-i))%*%cov%*%t(t(c(1,-1,rep(0,(I-2)+(2*(J-1))+i-1),1,rep(0,K-2),-1,rep(0,K-1-i)))))}

lastc=c(1,-1,rep(0,(I-2)+2*(J-1)),rep(-1,K-1),rep(1,K-1))%*%cov%*%t(t(c(1,-1,rep(0,(I-2)+2*(J-1)),rep(-1,K-1),rep(1,K-1))))

stc=c(sqrt(stc),sqrt(lastc)) inter1=rep(1,J)

par(mfrow=c(3,2))

plot(resultp,type="n",axes=FALSE,ann=FALSE,ylim=c(min(resultc-quant*stc),max(resultc+quant*stc)))

lines(resultp,col=4,lwd=2) axis(2, col.axis="blue", las=1) axis(1, at=1:10, lab=c("1959-1963",

"1964-1968",

"1969-1973",

•

"2004-2008"), col.axis="blue")

title(main="Apc ", font.main=4, col.main="red") title(xlab="Period", col.lab="red")

title(ylab="average log relative risk",col.lab="red") lines(resultp+quant*stp,lty=2,col=4,lwd=2) axis(2, col.axis="blue", las=1) axis(1, at=1:24, lab=c("1875-1883",

"1880-1888",

"1990-1998"), col.axis="blue", las=2)

title(main="Apc ", font.main=4, col.main="red")

title(ylab="average log relative risk",col.lab="red") lines(resultc+quant*stc,lty=2,col=4,lwd=2)

lines(resultc-quant*stc,lty=2,col=4,lwd=2) abline(0,0)

#########mapc-apC#####################

•

for(i in 1:(I*J)){if(A[i]!=I)x[i,A[i]]=1

else{if(A[i]==I)x[i,c(1:(I-1))]=-1}}

for(i in ((I*J)+1):(2*I*J)){if(A[i]!=I)x[i,I-1+A[i]]=1 else{if(A[i]==I)x[i,c(I:((I-1)*2))]=-1}}

for(i in 1:(I*J)){if(P[i]!=J)x[i,(2*(I-1)+P[i])]=1

else{if(P[i]==J)x[i,c((2*(I-1)+1):(2*(I-1)+J-1))]=-1}}

for(i in ((I*J)+1):(2*I*J)){if(P[i]!=J)x[i,2*(I-1)+J-1+P[i]]=1 else{if(P[i]==J)x[i,c((2*(I-1)+J):(2*(I-1)+2*(J-1)))]=-1}}

for(i in 1:(I*J*2)){if(C[i]!=K)x[i,((I-1)*2+(J-1)*2+C[i])]=1

else{if(C[i]==K)x[i,c(((I-1)*2+(J-1)*2+1):((I-1)*2+(J-1)*2+K-1))]=-1}}

m1=c(rep(1,(I*J)),rep(0,(I*J))) m2=c(rep(0,(I*J)),rep(1,(I*J))) x=cbind(m1,m2,x)

N=rbind(data12,data22) Y=rbind(data11,data21) y=c()

for(i in 1:(I*2)) {y=c(y,Y[i,])}

y=as.numeric(y) n=c()

for(i in 1:(I*2)) {n=c(n,N[i,])}

n=as.numeric(n)

############e1=e2###############

for(i in 1:(2*I*J)){

if(x[i,(58+2)]+x[i,(59+2)]==1)x[i,(58+2)]=1 else{if(x[i,(58+2)]+x[i,(59+2)]==-2)x[i,(58+2)]=-2 else{x[i,(58+2)]=0}}}

x=x[,-(59+2)]

lm1=glm.fit(x,y,family=poisson(),offset=log(n),intercept=F) coe=as.matrix(lm1$coefficients)

summary=summary.glm(lm1) cov=summary$cov.unscaled intercept=coe[1:2]

coe=coe[-c(1:2)]

coe=c(coe[1:58],coe[58],coe[(58+1):((I-1)*2+(J-1)*2+(K-1)-1)]) coea1=coe[1:(I-1)]

coea2=coe[((I-1)+1):((I-1)*2)]

coep1=coe[(2*(I-1)+1):(2*(I-1)+J-1)]

coep2=coe[((I-1)*2+J):(2*(I-1)+2*(J-1))]

coec=coe[(2*(I-1)+2*(J-1)+1):(2*(I-1)+2*(J-1)+K-1)]

coe1=c(intercept[1],coea1,0-sum(coea1),coep1,0-sum(coep1),coec,0-sum(coec)) coe2=c(intercept[2],coea2,0-sum(coea2),coep2,0-sum(coep2),coec,0-sum(coec))

•

result=coe1-coe2 resulta=c(result[2:16]) resultp=c(result[17:26])

cov=summary.glm(lm1,correlation = T)$cov.unscaled cov=(3614/df.residual(lm1))*cov

sta=c()

for(i in 1:(I-1)){sta=c(sta,c(0,0,rep(0,i-1),1,rep(0,I-2),-1,rep(0,(I-1)-i+2*(J-1)+K-2))%*%cov%*%t(t(c(0,0,rep(0,i-1),1,rep(0,I-2),-1,rep(0,(I-1)-i+2*(J-1)+K-2)))))}

lasta=c(0,0,rep(-1,I-1),rep(1,I-1),rep(0,2*(J-1)+(K-2)))%*%cov%*%t(t(c(0,0,rep(-1,I-1),rep(1,I-1),rep(0,2*(J-1)+(K-2)))))

sta=c(sqrt(sta),sqrt(lasta)) stp=c()

for(i in 1:(J-1)){stp=c(stp,c(0,0,rep(0,2*(I-1)+i-1),1,rep(0,J-2),-1,rep(0,K-2+J-1-i))%*%cov%*%t(t(c(0,0,rep(0,2*(I-1)+i-1),1,rep(0,J-2),-1,rep(0,K-2+J-1-i)))))}

lastp=c(0,0,rep(0,2*(I-1)),rep(-1,J-1),rep(1,J-1),rep(0,K-2))%*%cov%*%t(t(c(0,0,rep(0,2*(I-1)),rep(-1,J-1),rep(1,J-1),rep(0,K-2)))) stp=c(sqrt(stp),sqrt(lastp))

####plot####

plot(resulta,type="n",axes=FALSE,ann=FALSE,ylim=c(-3,3)) lines(resulta,col=5,lwd=2)

axis(2, col.axis="blue", las=1) axis(1, at=1:15, lab=c("10-14",

"15-19",

"80+"), col.axis="blue")

title(main="apC ", font.main=4, col.main="red") title(xlab="Age", col.lab="red")

title(ylab="difference",col.lab="red")

axis(2, col.axis="blue", las=1) axis(1, at=1:10, lab=c("1959-1963",

•

"2004-2008"), col.axis="blue")

title(main="apC ", font.main=4, col.main="red") title(xlab="Period", col.lab="red")

title(ylab="difference",col.lab="red")

cov=summary.glm(lm1,correlation = T)$cov.unscaled cov=(3614/df.residual(lm1))*cov

sta=c()

for(i in 1:(I-1)){sta=c(sta,c(1,-1,rep(0,i-1),1,rep(0,I-2),-1,rep(0,(I-1)-i+2*(J-1)+K-2))%*%cov%*%t(t(c(1,-1,rep(0,i-1),1,rep(0,I-2),-1,rep(0,(I-1)-i+2*(J-1)+K-2)))))}

lasta=c(1,-1,rep(-1,I-1),rep(1,I-1),rep(0,2*(J-1)+(K-2)))%*%cov%*%t(t(c(1,-1,rep(-1,I-1),rep(1,I-1),rep(0,2*(J-1)+(K-2)))))

sta=c(sqrt(sta),sqrt(lasta)) stp=c()

for(i in 1:(J-1)){stp=c(stp,c(1,-1,rep(0,2*(I-1)+i-1),1,rep(0,J-2),-1,rep(0,K-2+J-1-i))%*%cov%*%t(t(c(1,-1,rep(0,2*(I-1)+i-1),1,rep(0,J-2),-1,rep(0,K-2+J-1-i)))))}

lastp=c(1,-1,rep(0,2*(I-1)),rep(-1,J-1),rep(1,J-1),rep(0,K-2))%*%cov%*%t(t(c(1,-1,rep(0,2*(I-1)),rep(-1,J-1),rep(1,J-1),rep(0,K-2))))

stp=c(sqrt(stp),sqrt(lastp)) inter1=rep(1,I) axis(2, col.axis="blue", las=1) axis(1, at=1:15, lab=c("10-14",

"15-19",

•

"80+"), col.axis="blue")

title(main="apC ", font.main=4, col.main="red") title(xlab="Age", col.lab="red")

title(ylab="average log relative risk",col.lab="red") lines(resulta+quant*sta,lty=2,col=5,lwd=2) axis(2, col.axis="blue", las=1) axis(1, at=1:10, lab=c("1959-1963",

"1964-1968",

"2004-2008"), col.axis="blue")

title(main="apC ", font.main=4, col.main="red") title(xlab="Period", col.lab="red")

title(ylab="average log relative risk",col.lab="red") lines(resultp+quant*stp,lty=2,col=5,lwd=2)

lines(resultp-quant*stp,lty=2,col=5,lwd=2) abline(0,0)

####################conditional approach##############

library(VGAM)

##############data#####################

•

calculate_dispersion <- function(object){

# divide the sum of the squared Pearson residuals # by the residual degrees of freedom.

disp <- sum(resid(object,type="pearson")^2)/(df.residual(object))

return(disp) }

####################aPc########################

x=matrix(0,ncol=(I-1)+(K-1),nrow=I*J) for(i in 1:(I*J)){if( A[i]!=I)x[i,A[i]]=1 else{if(A[i]==I)x[i,c(1:(I-1))]=-1}}

for(i in 1:(I*J)){if(C[i]!=K)x[i,((I-1)+C[i])]=1 else{if(C[i]==K)x[i,c(I:(I-1+K-1))]=-1}}

###################fit model######################

offterm=log(m/f)

lm2=vgam(cbind(md,fd)~x, family=multinomial, offset=offterm, qr.arg=T) coe=coef(lm2)

##############estimate_difference###########

resulta=c(coe[-1][1:I-1],0-sum(coe[-1][1:I-1]))

resultc=c(coe[-1][I:(K+I-1-1)],0-sum(coe[-1][I:(K+I-1-1)]))

###############difference######################

#########stand.error#############

cov=vcov(lm2)

covmat.scaled <- calculate_dispersion(lm2)*cov sta=c()

for(i in

1:(I-1)){sta=c(sta,c(0,rep(0,i-1),1,rep(0,(I-1)-i+(K-1)))%*%covmat.scaled%*%t(t(c(0,rep(0,i-1),1,rep(0,(I-1)-i+(K-1))))))}

lasta=c(0,rep(-1,I-1),rep(0,K-1))%*%covmat.scaled%*%t(t(c(0,rep(-1,I-1),rep(0,K-1))))

sta=c(sqrt(sta),sqrt(lasta)) stc=c()

for(i in

1:(K-1)){stc=c(stc,c(0,rep(0,((I-1)+i-1)),1,rep(0,(K-1-i)))%*%covmat.scaled%*%t(t(c(0,rep(0,((I-1)+i-1)),1,rep(0,(K-1-i))))))}

lastc=c(0,rep(0,I-1),rep(-1,K-1))%*%covmat.scaled%*%t(t(c(0,rep(0,I-1),rep(-1,K-1))))

•

stc=c(sqrt(stc),sqrt(lastc)) se=sqrt(diag(covmat.scaled))

########difference###########

#############PLOT###################

plot(resulta,type="n",axes=FALSE,ann=FALSE,ylim=c(-2,1.5)) lines(resulta,col=3,lwd=2)

axis(2, col.axis="blue", las=1) axis(1, at=1:15, lab=c("10-14",

"15-19",

"80+"), col.axis="blue")

title(main="aPc ", font.main=4, col.main="red") title(xlab="Age", col.lab="red")

title(ylab="difference",col.lab="red")

axis(2, col.axis="blue", las=1) axis(1, at=1:24, lab=c("1875-1883",

"1880-1888",

•

"1990-1998"), col.axis="blue", las=2)

title(main="aPc ", font.main=4, col.main="red")

resulta=coe[1]+c(coe[-1][1:I-1],0-sum(coe[-1][1:I-1]))

resultc=coe[1]+c(coe[-1][I:(K+I-1-1)],0-sum(coe[-1][I:(K+I-1-1)]))

#########stand.error_平均相對風險############

cov=vcov(lm2)

covmat.scaled <- calculate_dispersion(lm2)*cov sta=c()

for(i in

1:(I-1)){sta=c(sta,c(1,rep(0,i-1),1,rep(0,(I-1)-i+(K-1)))%*%covmat.scaled%*%t(t(c(1,rep(0,i-1),1,rep(0,(I-1)-i+(K-1))))))}

lasta=c(1,rep(-1,I-1),rep(0,K-1))%*%covmat.scaled%*%t(t(c(1,rep(-1,I-1),rep(0,K-1))))

sta=c(sqrt(sta),sqrt(lasta)) stc=c()

for(i in

1:(K-1)){stc=c(stc,c(1,rep(0,((I-1)+i-1)),1,rep(0,(K-1-i)))%*%covmat.scaled%*%t(t(c(1,rep(0,((I-1)+i-1)),1,rep(0,(K-1-i))))))}

lastc=c(1,rep(0,I-1),rep(-1,K-1))%*%covmat.scaled%*%t(t(c(1,rep(0,I-1),rep(-1,K-1))))

stc=c(sqrt(stc),sqrt(lastc))

#########qaic#############

dispersion=3442/208

qaic_1 <- -2*logLik(lm2)/dispersion + 2*(length(resid(lm2)) - df.residual(lm2)+1)

#############PLOT###################

par(mfrow=c(3,2))

plot(resulta,type="n",axes=FALSE,ann=FALSE,ylim=c(min(resultc-quant*stc),max(resultc+quant*stc)))

lines(resulta,col=3,lwd=2) axis(2, col.axis="blue", las=1) axis(1, at=1:15, lab=c("10-14",

"15-19",

"20-24",

"25-29",

"30-34",

•

"80+"), col.axis="blue")

title(main="aPc ", font.main=4, col.main="red") title(xlab="Age", col.lab="red")

title(ylab="average log relative risk",col.lab="red") lines(resulta+quant*sta,lty=2,col=3,lwd=2) axis(2, col.axis="blue", las=1) axis(1, at=1:24, lab=c("1875-1883",

"1880-1888",

"1990-1998"), col.axis="blue", las=2)

title(main="aPc ", font.main=4, col.main="red")

title(ylab="average log relative risk",col.lab="red") lines(resultc+quant*stc,lty=2,col=4,lwd=2)

lines(resultc-quant*stc,lty=2,col=4,lwd=2)

•

####################conditional Apc##############

x=matrix(0,ncol=(J-1)+(K-1),nrow=I*J) for(i in 1:(I*J)){if( P[i]!=J)x[i,P[i]]=1 else{if(P[i]==J)x[i,c(1:(J-1))]=-1}}

for(i in 1:(I*J)){if(C[i]!=K)x[i,((J-1)+C[i])]=1 else{if(C[i]==K)x[i,c(J:(J-1+K-1))]=-1}}

###################fit model######################

offterm=log(m/f)

lm2=vgam(cbind(md,fd)~x, family=multinomial, offset=offterm, qr.arg=T) coe=coef(lm2)

#########difference#######

########estimate#######

resultp=c(coe[-1][1:J-1],0-sum(coe[-1][1:J-1]))

resultc=c(coe[-1][J:(K+J-1-1)],0-sum(coe[-1][J:(K+J-1-1)]))

#########stand.error#############

cov=vcov(lm2)

covmat.scaled <- calculate_dispersion(lm2)*cov stp=c()

for(i in

1:(J-1)){stp=c(stp,c(0,rep(0,i-1),1,rep(0,K-1+J-1-i))%*%covmat.scaled%*%t(t(c(0,rep(0,i-1),1,rep(0,K-1+J-1-i)))))}

lastp=c(0,rep(-1,J-1),rep(0,K-1))%*%covmat.scaled%*%t(t(c(0,rep(-1,J-1),rep(0,K-1))))

stp=c(sqrt(stp),sqrt(lastp)) stc=c()

for(i in

1:(K-1)){stc=c(stc,c(0,rep(0,J-1+i-1),1,rep(0,K-1-i))%*%covmat.scaled%*%t(t(c(0,rep(0,J-1+i-1),1,rep(0,K-1-i)))))}

lastc=c(0,rep(0,J-1),rep(-1,K-1))%*%covmat.scaled%*%t(t(c(0,rep(0,J-1),rep(-1,K-1))))

stc=c(sqrt(stc),sqrt(lastc)) se=sqrt(diag(covmat.scaled))

###########plot#############

par(mfrow=c(3,2))

plot(resultp,type="n",axes=FALSE,ann=FALSE,ylim=c(-2,1.5)) lines(resultp,col=4,lwd=2)

axis(2, col.axis="blue", las=1) axis(1, at=1:10, lab=c("1959-1963",

"1964-1968",

"2004-2008"), col.axis="blue")

title(main="Apc ", font.main=4, col.main="red") title(xlab="Period", col.lab="red")

•

axis(2, col.axis="blue", las=1) axis(1, at=1:24, lab=c("1875-1883",

"1880-1888",

"1990-1998"), col.axis="blue", las=2)

title(main="Apc ", font.main=4, col.main="red") title(ylab="difference",col.lab="red")

resultp=coe[1]+c(coe[-1][1:J-1],0-sum(coe[-1][1:J-1]))

resultc=coe[1]+c(coe[-1][J:(K+J-1-1)],0-sum(coe[-1][J:(K+J-1-1)]))

#########stand.error#############

cov=vcov(lm2)

covmat.scaled <- calculate_dispersion(lm2)*cov stp=c()

for(i in

1:(J-1)){stp=c(stp,c(1,rep(0,i-1),1,rep(0,K-1+J-1-i))%*%covmat.scaled%*%t(t(c(1,rep(0,i-1),1,rep(0,K-1+J-1-i)))))}

•

stp=c(sqrt(stp),sqrt(lastp)) stc=c()

for(i in

1:(K-1)){stc=c(stc,c(1,rep(0,J-1+i-1),1,rep(0,K-1-i))%*%covmat.scaled%*%t(t(c(1,rep(0,J-1+i-1),1,rep(0,K-1-i)))))}

lastc=c(1,rep(0,J-1),rep(-1,K-1))%*%covmat.scaled%*%t(t(c(1,rep(0,J-1),rep(-1,K-1))))

stc=c(sqrt(stc),sqrt(lastc))

#########qaic#############

dispersion=3442/208

qaic_2<- -2*logLik(lm2)/dispersion + 2*(length(resid(lm2)) - df.residual(lm2)+1)

###########plot#############

par(mfrow=c(3,2))

plot(resultp,type="n",axes=FALSE,ann=FALSE,ylim=c(min(resultc-quant*stc),max(resultc+quant*stc)))

lines(resultp,col=4,lwd=2) axis(2, col.axis="blue", las=1) axis(1, at=1:10, lab=c("1959-1963",

"1964-1968",

"2004-2008"), col.axis="blue")

title(main="Apc ", font.main=4, col.main="red") title(xlab="Period", col.lab="red")

title(ylab="average log relative risk",col.lab="red") lines(resultp+quant*stp,lty=2,col=4,lwd=2) axis(2, col.axis="blue", las=1) axis(1, at=1:24, lab=c("1875-1883",

"1880-1888",

•

"1990-1998"), col.axis="blue", las=2)

title(main="Apc ", font.main=4, col.main="red") title(ylab="average log relative risk",col.lab="red") lines(resultc+quant*stc,lty=2,col=4,lwd=2)

lines(resultc-quant*stc,lty=2,col=4,lwd=2) abline(0,0)

####################apC#######################

x=matrix(0,ncol=(I-1)+(J-1),nrow=I*J) for(i in 1:(I*J)){if( A[i]!=I)x[i,A[i]]=1 else{if(A[i]==I)x[i,c(1:(I-1))]=-1}}

for(i in 1:(I*J)){if(P[i]!=J)x[i,((I-1)+P[i])]=1 else{if(P[i]==J)x[i,c(I:(I-1+J-1))]=-1}}

###################fit model######################

offterm=log(m/f)

lm2=vgam(cbind(md,fd)~x, family=multinomial, offset=offterm, qr.arg=T) coe=coef(lm2)

###############difference##########

###########estimate###############

resulta=c(coe[-1][1:I-1],0-sum(coe[-1][1:I-1]))

resultp=c(coe[-1][I:(I+J-1-1)],0-sum(coe[-1][I:(I+J-1-1)])) #########stand.error#############

cov=vcov(lm2)

covmat.scaled <- calculate_dispersion(lm2)*cov sta=c()

for(i in

1:(I-1)){sta=c(sta,c(0,rep(0,i-1),1,rep(0,(I-1)-i+(J-1)))%*%covmat.scaled%*%t(t(c(0,rep(0,i-1),1,rep(0,(I-1)-i+(J-1))))))}

lasta=c(0,rep(-1,I-1),rep(0,J-1))%*%covmat.scaled%*%t(t(c(0,rep(-1,I-1),rep(0,J-1))))

sta=c(sqrt(sta),sqrt(lasta)) stp=c()

for(i in

1:(J-1)){stp=c(stp,c(0,rep(0,(I-1)+i-1),1,rep(0,J-1-•

stp=c(sqrt(stp),sqrt(lastp)) se=sqrt(diag(covmat.scaled))

#######plot#########

plot(resulta,type="n",axes=FALSE,ann=FALSE,ylim=c(-2,1.5)) lines(resulta,col=5,lwd=2)

axis(2, col.axis="blue", las=1) axis(1, at=1:15, lab=c("10-14",

"15-19",

"80+"), col.axis="blue")

title(main="apC ", font.main=4, col.main="red") title(xlab="Age", col.lab="red")

title(ylab="difference",col.lab="red")

axis(2, col.axis="blue", las=1) axis(1, at=1:10, lab=c("1959-1963",

"1964-1968",

"2004-2008"), col.axis="blue")

title(main="apC ", font.main=4, col.main="red") title(xlab="Period", col.lab="red")

title(ylab="difference",col.lab="red")

lines(resultp+quant*stp,lty=2,col=5,lwd=2)

lines(resultp-quant*stp,lty=2,col=5,lwd=2)

•

resulta=coe[1]+c(coe[-1][1:I-1],0-sum(coe[-1][1:I-1]))

resultp=coe[1]+c(coe[-1][I:(I+J-1-1)],0-sum(coe[-1][I:(I+J-1-1)])) #########stand.error#############

cov=vcov(lm2)

covmat.scaled <- calculate_dispersion(lm2)*cov sta=c()

for(i in

1:(I-1)){sta=c(sta,c(1,rep(0,i-1),1,rep(0,(I-1)-i+(J-1)))%*%covmat.scaled%*%t(t(c(1,rep(0,i-1),1,rep(0,(I-1)-i+(J-1))))))}

lasta=c(1,rep(-1,I-1),rep(0,J-1))%*%covmat.scaled%*%t(t(c(1,rep(-1,I-1),rep(0,J-1))))

sta=c(sqrt(sta),sqrt(lasta)) stp=c()

for(i in

1:(J-1)){stp=c(stp,c(1,rep(0,(I-1)+i-1),1,rep(0,J-1-i))%*%covmat.scaled%*%t(t(c(1,rep(0,(I-1)+i-1),1,rep(0,J-1-i)))))}

lastp=c(1,rep(0,I-1),rep(-1,J-1))%*%covmat.scaled%*%t(t(c(1,rep(0,I-1),rep(-1,J-1))))

stp=c(sqrt(stp),sqrt(lastp)) #########qaic#############

dispersion=3442/208

qaic_3<- -2*logLik(lm2)/dispersion + 2*(length(resid(lm2))- df.residual(lm2)+1)

#######plot#########

plot(resulta,type="n",axes=FALSE,ann=FALSE,ylim=c(min(resulta-quant*sta),max(resulta+quant*sta)))

lines(resulta,col=5,lwd=2) axis(2, col.axis="blue", las=1) axis(1, at=1:15, lab=c("10-14",

"15-19",

"80+"), col.axis="blue")

title(main="apC ", font.main=4, col.main="red") title(xlab="Age", col.lab="red")

title(ylab="average log relative risk",col.lab="red")

lines(resulta+quant*sta,lty=2,col=5,lwd=2)

• 國

立 政 治 大

㈻㊫學

•

N a tio na

l C h engchi U ni ve rs it y

lines(resulta-quant*sta,lty=2,col=5,lwd=2) abline(0,0)

plot(resultp,type="n",axes=FALSE,ann=FALSE,ylim=c(min(resulta-quant*sta),max(resulta+quant*sta)))

lines(resultp,col=5,lwd=2) axis(2, col.axis="blue", las=1) axis(1, at=1:10, lab=c("1959-1963",

"1964-1968",

"1969-1973",

"1974-1978",

"1979-1983",

"1984-1988",

"1989-1993",

"1994-1998",

"1999-2003",

"2004-2008"), col.axis="blue")

title(main="apC ", font.main=4, col.main="red") title(xlab="Period", col.lab="red")

title(ylab="average log relative risk",col.lab="red") lines(resultp+quant*stp,lty=2,col=5,lwd=2)

lines(resultp-quant*stp,lty=2,col=5,lwd=2)

abline(0,0)

在文檔中 台灣地區男女自殺死亡率之比較研究 - 政大學術集成 (頁 51-80)