第五章 討論與結論
A. Lu & Ades model 和貢獻圖
function(data){
N.study = max(data[,'s']) # 儲存必要資訊
N.tr = max(data[,'t']) # 試驗數與治療數
data = data[,c('s','t','n','y','std')] # 保留必要資訊
data = cbind(data, se = data[,'std']/sqrt(data[,'n'])) # 生成 Standard error
No.study = c(); count = 1; # 準備轉成 Contrast-based 的格式
while(count < N.study+1){ # 為了處理 multi-arm 的 variance 而做
Ntrials = length(which(data[,'s'] == count)) – 1
No.study = c(No.study,rep(count,Ntrials)) # 先把每個試驗貢獻幾個相對效果算好
count = count + 1
}
tr.index = array(0,dim=c(length(No.study),N.tr)) # 紀錄每一列是由哪兩個治療相減
tr.e = array(0,dim=c(length(No.study),N.tr)) # 紀錄對應的治療效果
dXY = array(,dim=c(length(No.study),6)) # Contrast-based 的格式
colnames(dXY) = c('study', 'design', 'N.arm' ,'d','Baseline','Comparison')
dXY[,'study'] = No.study; design = 1
V = array(0,dim=c(length(No.study),length(No.study))) # variance 的格式預先設定
for(i in 1:N.study){ # 轉成 Contrast-based 的格式
data.i = which(data[,'s']==i) # 標記試驗在輸入的資料裡的位置
dXY.i = which(dXY[,'study']==i) # 標記試驗在 Contrast-based 的格式的位置
50
L = length(dXY.i) # 試驗 i 提供的相對治療效果數目
dXY[dXY.i,'N.arm'] = rep(L+1,L) # 紀錄 arm 數
temp = data[data.i,'t']; temp.order = order(temp); temp = temp[temp.order] # 抓取資料
temp.y = data[data.i,'y']; temp.y = temp.y[temp.order]
tr.index[dXY.i,temp[1]] = rep(-1,L) # 紀錄基底治療
tr.e[dXY.i,temp[1]] = rep(temp.y[1],L) # 紀錄基底治療的效果
for(tag in 1:L){
tr.index[dXY.i[tag],temp[tag+1]] = 1 # 紀錄成對比較的另一個治療
tr.e[dXY.i[tag],temp[tag+1]] = temp.y[tag+1] # 紀錄成對比較的另一個治療的效果
}
count = 1
for(j in dXY.i){
count = count + 1
dXY[j,'Baseline'] = temp[1] # 在格式中紀錄基底治療
dXY[j,'Comparison'] = temp[count] # 在格式中紀錄成對比較的另一個治療
}
if(i > 1){ # 根據比較的治療分組
for(f in rev(2:i)-1){
dXY.f = which(dXY[,'study']==f)
if(L==length(dXY.f)){
if(sum((tr.index[dXY.i,]-tr.index[dXY.f,])!=0) == 0){
dXY[dXY.i,'design'] = dXY[dXY.f,'design']
}
}
}
51 if(is.na(dXY[dXY[dXY.i[1],'design']]))
design = design + 1; dXY[dXY.i,'design'] = rep(design, L)
}
} else dXY[dXY.i,'design'] = design
temp.se = data[data.i,'se']; temp.se = temp.se[temp.order] # Variance = (Standard error)^2
temp.v = array(0, dim=c(L,L)); temp.v = temp.v + temp.se[1]^2
temp.v = temp.v + diag(L) * (temp.se[2:(1+L)] ^ 2) # 兩個治療的 Standard error 相加
V[dXY.i,dXY.i] = temp.v # 塞入對應的位置
}
dXY[,'d'] = diag(tr.index %*% t(tr.e)) # 計算相對治療效果
W = solve(V) # 權重為 variance 的倒數
net.X = tr.index[,-1] # 拿掉參考治療，將格式轉成 Contrast-based
net.d = solve(t(net.X)%*%W%*%net.X) %*% t(net.X)%*%W%*%dXY[,'d'] # Fixed effect
net.cov = solve(t(net.X)%*%W%*%net.X)
net.sd = sqrt(diag(net.cov))
# Random effect net.Q=t(dXY[,'d']-net.X%*%net.d) %*%W%*%(dXY[,'d']-net.X%*%net.d) # (Pearson residual)^2
net.df = nrow(dXY) - (N.tr-1) # 異質性自由度計算
net.C = rep(0,max(dXY[,'design'])) # 異質性分母計算
for(i in 1:length(net.C)){
index.g = dXY[,'design']==i
temp.g = dXY[index.g,]; index.g.s = unique(temp.g[,'study'])
sum.W=0; sum.W2=0
for(j in index.g.s){
index.g.i = dXY[,'study']==j
52 sum.W = sum.W + W[index.g.i,index.g.i]
sum.W2 = sum.W2 + W[index.g.i,index.g.i] %*% W[index.g.i,index.g.i]
}
net.C[i] = sum(diag(sum.W - sum.W2 %*% solve(sum.W))) # 取 trace
}; net.C = sum(net.C)
net.tau2.pool = (net.Q - net.df)/net.C; net.tau2.pool = as.numeric(net.tau2.pool) # tau^2
if(net.tau2.pool<0){
net.tau2.pool = 0
}
V.r = V + diag(net.tau2.pool/2,length(No.study)); V.r[V.r!=0] = V.r[V.r!=0] + net.tau2.pool/2
W.r = solve(V.r) # Random effect model 的權重
net.d.r = solve(t(net.X)%*%W.r%*%net.X) %*% t(net.X)%*%W.r%*%dXY[,'d']
net.cov.r = solve(t(net.X)%*%W.r%*%net.X)
net.sd.r = sqrt(diag(net.cov.r))
# 貢獻圖
LA.data = format.LA(data) # 由長型轉成寬型的資料格式
rL = sum(LA.data[,'na']*(LA.data[,'na']-1)/2) # 計算全部試驗貢獻多少直接證據
w.dXY = array(0,dim=c(rL,5+N.tr)) # 貢獻圖版本的 Contrast-based 資料格式
colnames(w.dXY) = c('study','d', 'se','Baseline','Comparison',1:N.tr) # 重新命名欄位名稱
t = colnames(LA.data)[seq(1,ncol(LA.data)-1,4)]
n = colnames(LA.data)[seq(2,ncol(LA.data)-1,4)]
y = colnames(LA.data)[seq(3,ncol(LA.data)-1,4)]
std = colnames(LA.data)[seq(4,ncol(LA.data)-1,4)]
count = 0
for(i in 1:nrow(LA.data)){ # 開始紀錄每個試驗貢獻的相對治療效果
53
L = LA.data[i,'na']*(LA.data[i,'na']-1)/2 # 試驗 i 貢獻的相對治療效果數目
temp = matrix(w.dXY[1:L+count,],L)
colnames(temp) = colnames(w.dXY)
temp.c = 1
for(j in 1:(LA.data[i,'na']-1)){
for(k in (j+1):LA.data[i,'na']){
temp[temp.c,'study'] = i # 紀錄試哪一個試驗
temp[temp.c,'Baseline'] = LA.data[i,t[j]] # 紀錄基底治療
temp[temp.c,'Comparison'] = LA.data[i,t[k]] # 紀錄成對比較的另一個治療
temp[temp.c,length(colnames(w.dXY))-N.tr+LA.data[i,t[j]]] = -1 # 以基礎參數紀錄
temp[temp.c,length(colnames(w.dXY))-N.tr+LA.data[i,t[k]]] = 1
temp[temp.c,'d'] = LA.data[i,y[k]] - LA.data[i,y[j]] # 紀錄相對治療效果
temp[temp.c,'se'] = LA.data[i,std[k]]/sqrt(LA.data[i,n[k]]) # 紀錄 Standard error
+ LA.data[i,std[j]]/sqrt(LA.data[i,n[j]])
temp.c = temp.c + 1
}
}
w.dXY[1:L+count,] = temp # 將結果塞到正確的位置
count = count + L
}
w.dXY = w.dXY[order(w.dXY[,'Comparison']),] # 整理資料順序
w.dXY = w.dXY[order(w.dXY[,'Baseline']),]
direct = unique(w.dXY[,c('Baseline','Comparison',as.character(1:N.tr))]) # 直接證據的種類數
direct.d = c(); direct.v = c(); dir.trial = c() # 直接證據
for(i in 1:nrow(direct)){ # 每一種成對比較分開算
54 index = c()
for(j in 1:nrow(w.dXY)){ # 標記是哪一種成對比較
if(identical(w.dXY[j,c('Baseline','Comparison')],direct[i,c('Baseline','Comparison')])){
index = c(index,j)
}
}
dir.trial = c(dir.trial, length(index)) # 提供的相對治療效果數
temp.d = w.dXY[index,'d'] # Fixed effect 相對治療效果
temp.se = w.dXY[index,'se'] # Standard error
temp.w = 1/temp.se^2 # 權重
pooled.d = sum(temp.d*temp.w)/sum(temp.w) # 合併
temp.Q = sum(temp.w*(temp.d-pooled.d)^2) # (Pearson residual)^2
temp.tau2 = (temp.Q - (length(index) - 1)) / (sum(temp.w) - (sum(temp.w^2)/sum(temp.w)))
if(temp.tau2 < 0){ # tau^2
temp.tau2 = 0
}
temp.w.r = 1/(temp.se^2 + temp.tau2) # Random effect 權重
pooled.d.r = sum(temp.d*temp.w.r)/sum(temp.w.r) # 相對治療效果
pooled.v.r = 1/sum(temp.w.r) # variance
direct.d = c(direct.d,pooled.d.r) # 紀錄每個成對比較各自合併後的結果
direct.v = c(direct.v,pooled.v.r) # 紀錄每個成對比較各自合併後的 variance
}; direct.v = diag(direct.v) # variance 轉成
total = array(0,dim=c(N.tr*(N.tr-1)/2,N.tr+3)) # 計算 N.tr 種治療，最多有哪幾種成對比較
colnames(total) = c(1:N.tr,'Baseline','Comparison','Exist')
55 count=1
for(j in 1:(N.tr-1)){
for(k in (j+1):N.tr){
total[count,'Baseline'] = j
total[count,'Comparison'] = k
total[count,j] = -1
total[count,k] = 1
count = count + 1
}
}
design.matrix = direct[,as.character(2:N.tr)]; r.name=rep(NA,nrow(total)) # 轉成 Contrast-basrd
tag = 1
for(i in 1:nrow(total)){ # 紀錄是何種成對比較
if(identical(total[i,as.character(2:N.tr)],design.matrix[tag,])){
r.name[tag] = paste0(re.trtname(total[i,'Baseline']),'vs',re.trtname(total[i,'Comparison']))
if(tag<nrow(design.matrix)){
tag = tag + 1
}
}else{
design.matrix = rbind(design.matrix,total[i,as.character(2:N.tr)])
r.name[nrow(design.matrix)] =
paste0(re.trtname(total[i,'Baseline']),'vs',re.trtname(total[i,'Comparison']))
}
} # 資料準備齊全，最後計算
Xa = design.matrix[1:nrow(direct),]; rownames(Xa) = r.name[1:nrow(Xa)]#; colnames(Xa) =
56 r.name[1:ncol(Xa)]
colnames(direct.v) = r.name[1:ncol(direct.v)]; rownames(direct.v) = r.name[1:ncol(direct.v)]
Wa = solve(direct.v) # Hat 矩陣
H = design.matrix %*% solve(t(Xa)%*%Wa%*%Xa) %*% t(Xa) %*% Wa; H = round(H,3)
rownames(H) = r.name
Contribution = round(100*abs(H)/rowSums(abs(H)),1) # 貢獻比例
Overall.Contribution = round(100*colSums(Contribution)/sum(Contribution),1)
# 輸出 result = list(d.f = net.d, sd.f = net.sd, d.r = net.d.r, sd.r = net.sd.r, tau = sqrt(net.tau2.pool),
d.r.all = total[,2:N.tr] %*% net.d.r,
sd.r.all = sqrt(diag(total[,2:N.tr] %*% net.cov.r %*% t(total[,2:N.tr]))),
direct.comparison = cbind(contrast = r.name[1:ncol(H)], Number = dir.trial) ,
direct.var = round(direct.v,4), design.matrix = design.matrix, H = H,
Contribution = Contribution, Overall.Contribution = Overall.Contribution)
return(result)
}
B. Baseline model
function(data){
N.study = max(data[,'s']) # 抓取必要資訊
N.tr = max(data[,'t'])
data = data[,c('s','t','n','y','std')]
data = cbind(data, se = data[,'std']/sqrt(data[,'n'])) # 計算 Standard error
# 參考治療效果計算
ref = RofA(1,data) # 計算參考治療的加權平均治療效果與 variance
ref.mu = ref$mu.r # 參考治療加權平均治療效果
57
ref.sigma2 = ref$sigma2 # 參考治療 variance
# 捕入參考治療效果 data.add = data
for(i in unique(data[,'s'])){
index = which(data[,'s']==i); temp = data[index,]
if(!is.element(1,temp[,'t'])){
data.add = rbind(data.add, c(i,1,1,ref.mu,se1(i,data),se1(i,data)))
}
}; re = order(data.add[,'t']); data.add = data.add[re,]
re = order(data.add[,'s']); data.add = data.add[re,]
# 與 Lu & Ades model 相同，轉換成 Contrast-based 格式 tr.index = array(0,dim=c(nrow(data.add),N.tr))
dXY = array(,dim=c(nrow(data.add),6))
colnames(dXY) = c('study', 'design', 'N.arm' ,'d','Baseline','Comparison')
dXY[,'study'] = data.add[,'s']; design = 1
V = array(0,dim=c(nrow(data.add),nrow(data.add)))
for(i in 1:N.study){
tag = which(dXY[,'study']==i) # 標記位置
L = length(tag) # 相對治療效果數目
dXY[tag,'N.arm'] = rep(L,L)
temp = data.add[tag,'t']; temp.order = order(temp); temp = temp[temp.order] # 抓取試驗資料
temp.y = data.add[tag,'y']; temp.y = temp.y[temp.order]
tr.index[tag[1],temp[1]] = 1 # 紀錄基底治療(都是參考治療)
tr.index[tag[2:L],temp[1]] = rep(-1,L-1)
tr.index[tag[2:L],temp[2:L]] = diag(rep(1,L-1))
58
dXY[tag[1], 'd'] = temp.y[1] # 紀錄成對比較的另一個治療
dXY[tag[2:L],'d'] = temp.y[2:L] - temp.y[1]
dXY[tag,'Baseline'] = rep(temp[1],L)
dXY[tag[1],'Comparison'] = 0
dXY[tag[2:L],'Comparison'] = temp[2:L]
# 計算是何種分組 if(i > 1){
for(f in rev(2:i)-1){
dXY.f = which(dXY[,'study']==f)
if(L==length(dXY.f)){
if(sum((tr.index[tag,]-tr.index[dXY.f,])!=0) == 0){
dXY[tag,'design'] = dXY[dXY.f,'design']
}
}
}
if(is.na(dXY[dXY[tag[1],'design']])){ # 好像要改成 dXY[dXY.i[1],'design']
design = design + 1; dXY[tag,'design'] = rep(design, L)
}
} else dXY[tag,'design'] = design
# Variance 計算 temp.se = data.add[tag,'se']; temp.se = temp.se[temp.order]
temp.v = array(0, dim=c(L,L)); temp.v = temp.v + temp.se[1]^2
temp.v[2:L,2:L] = temp.v[2:L,2:L] + diag(L-1) * (temp.se[2:L] ^ 2)
if(temp[1]==1){
V[tag,tag] = temp.v + ref.sigma2
59 }else{
V[tag,tag] = temp.v
}
}; del = -which(dXY[,'Comparison']==0 & dXY[,'Baseline']==1) # 拿掉參考治療
dXY = dXY[del,] # 轉成 Contrast-based 格式
tr.index = tr.index[del,-1]
V = V[del,del] # 相對治療效果的 Variance-covariance matrix
W = solve(V) # 權重
net.X = tr.index # Fixed effect
net.d = solve(t(net.X)%*%W%*%net.X) %*% t(net.X)%*%W%*%dXY[,'d']
net.cov = solve(t(net.X)%*%W%*%net.X)
net.sd = sqrt(diag(net.cov))
# Random effect net.Q = t(dXY[,'d']-net.X%*%net.d) %*% W %*% (dXY[,'d']-net.X%*%net.d)
net.df = nrow(dXY) - (N.tr-1)
net.C = rep(0,max(dXY[,'design']))
for(i in 1:length(net.C)){
index.g = which(dXY[,'design']==i)
temp.g = dXY[index.g,]; index.g.s = unique(temp.g[,'study'])
sum.W=0; sum.W2=0
for(j in index.g.s){
index.g.i = which(dXY[,'study']==j)
sum.W = sum.W + W[index.g.i,index.g.i]
sum.W2 = sum.W2 + W[index.g.i,index.g.i] %*% W[index.g.i,index.g.i]
60 }
net.C[i] = sum(diag(sum.W - sum.W2 %*% solve(sum.W)))
}; net.C = sum(net.C)
net.tau2.pool = (net.Q - net.df)/net.C; net.tau2.pool = as.numeric(net.tau2.pool)
if(net.tau2.pool<0){
net.tau2.pool = 0
}
V.r = V + diag(net.tau2.pool/2,nrow(dXY))
V.r[V.r!=0] = V.r[V.r!=0] + net.tau2.pool/2
W.r = solve(V.r)
net.d.r = solve(t(net.X)%*%W.r%*%net.X) %*% t(net.X)%*%W.r%*%dXY[,'d']
net.cov.r = solve(t(net.X)%*%W.r%*%net.X)
net.sd.r = sqrt(diag(net.cov.r))
# 用來將結果轉換成全部的成對比較估計結果 total = array(0,dim=c(N.tr*(N.tr-1)/2,N.tr))
colnames(total) = c(1:N.tr)
count=1
for(j in 1:(N.tr-1)){
for(k in (j+1):N.tr){
total[count,j] = -1
total[count,k] = 1
count = count + 1
}
}
# 輸出
61
result = list(d.f = net.d, sd.f = net.sd, d.r = net.d.r, sd.r = net.sd.r,
d.r.all = total[,2:N.tr] %*% net.d.r,
sd.r.all = sqrt(diag(total[,2:N.tr] %*% net.cov.r %*% t(total[,2:N.tr]))),
tau = sqrt(net.tau2.pool))
return(result)
}
C. Arm-based model
function(data){
N.study = max(data[,'s']) # 抓取資訊
N.tr = max(data[,'t'])
data = data[,c('s','t','n','y','std')]
data = cbind(data, se = data[,'std']/sqrt(data[,'n'])) # 計算 Standard error
V = diag(data[,'se']^2) # Variance
W = solve(V) # 權重
sigma2 = rep(0,N.tr)
for(i in 1:N.tr){
sigma2[i] = RofA(i,data)$sigma2 # 每個治療各自的隨機效應的 variance
}
# Random effect V.r = V
X = array(0, dim = c(nrow(data),N.tr))
for(i in 1:nrow(data)){
X[i,data[i,'t']] = 1
V.r[i,i] = V.r[i,i] + sigma2[data[i,'t']]
62 }
W.r = solve(V.r)
A.r = solve(t(X) %*% W.r %*% X) %*% t(X) %*% W.r %*% data[,'y'] # 絕對治療效果
covA.r = solve(t(X) %*% W.r %*% X) # 絕對治療效果估計值的 variacne-covariance matrix
d.r = A.r[2:N.tr] - A.r[1] # 將絕對治療效果轉成相對治療效果
sd.r = sqrt(diag(covA.r)[2:N.tr] + covA.r[1,1])
# 用來將結果轉換成全部的成對比較估計結果 total = array(0,dim=c(N.tr*(N.tr-1)/2,N.tr))
colnames(total) = c(1:N.tr)
count=1
for(j in 1:(N.tr-1)){
for(k in (j+1):N.tr){
total[count,j] = -1
total[count,k] = 1
count = count + 1
}
}
# 輸出 result = list(d.f = d, sd.f = sd, ref.mu.f = A[1], ref.sigma2.f = covA[1,1],
d.r = d.r, sd.r = sd.r, ref.mu.r = A.r[1], ref.sigma2.r = covA.r[1,1],
d.r.all = total %*% A.r,
sd.r.all = sqrt(diag(total %*% covA.r %*% t(total))))
return(result)
}