サクサクとメインの表を作っていく作業はダイナミックで楽しいところ
まずはLDLの変化量をありのまま出してみる 対応のあるt検定
> with(Dataset, (t.test(LDL0, LDL3, alternative='two.sided', conf.level=.95, paired=TRUE)))
Paired t-test
data: LDL0 and LDL3
t = -2.748, df = 12, p-value = 0.01767
alternative hypothesis: true mean difference is not equal to 0
95 percent confidence interval:
-29.789220 -3.441549
sample estimates:
mean difference
-16.61538
Paired t-test
P=0.01767 よし、とりあえず差はありそうだ。
いや、新しい変数をつくってみよう 3か月でのLDL上昇
Dataset$LDL_diff <- with(Dataset, LDL3- LDL0)
ヒストグラムをつくってみよう。
with(Dataset, Hist(LDL_diff, scale="frequency", breaks="Sturges", col="darkgray"))
併用ステロイドの減量によってならLDLは下がるはずだが上がっているな。
交絡と考えているPSL、CRPの変化も変数として作成。
Dataset$PSL_diff <- with(Dataset, PSL3- PSL0)
Dataset$CRP_diff <- with(Dataset, CRP3- CRP0)
PSLの変化量を箱ひげ図にしてみよう。
Boxplot( ~ PSL_diff, data=Dataset, id=list(method="y"))
回帰分析してみようとしたらIL6阻害薬が認識されていない。
IL6阻害薬をカテゴリカル変数に指定。
editDataset(Dataset)
Dataset <- within(Dataset, {
IL6 <- as.factor(IL6)
})
うまくいかない。ダミー変数をつくろう。TCZ、SARをデータセットに直接編集して追加。
そもそもモデルの作り方が間違っているな。
LDLの変化量=IL6阻害薬の種類+CRP変化量+PSL変化量
これじゃ、IL6阻害薬の種類の比較になってしまうので、ダミー変数を使おうがそこにはほとんど差がないという結果しかでないよな。
さて、ではどのようにしよう。岡山大学の津田先生から多変量解析の出発点は層別分析だ。丁寧に解析してみなさい。そう言われたことを思い出した。
PSL使用例8例と未使用例5例に分けてまずLDL変化量を解析しなおしてみよう。
> numSummary(Dataset[,"LDL_diff", drop=FALSE], groups=Dataset$PSLumu, statistics=c("mean", "sd", "IQR", "quantiles"), quantiles=c(0,
+ .25,.5,.75,1))
mean sd IQR 0% 25% 50% 75% 100% LDL_diff:n
ari 15.25 22.35269 36.5 -17 -4.75 25 31.75 39 8
nasi 18.80 23.27445 21.0 -3 2.00 16 23.00 56 5
とりあえずやはりPSLがない症例のほうがLDL上昇は高そうにみえる。これはPSL併用症例はPSL減量によりLDL低下効果がありそうにみえる。とりあえず定性的に評価してみる。
> t.test(LDL_diff~PSLumu, alternative='two.sided', conf.level=.95, var.equal=FALSE, data=Dataset)
Welch Two Sample t-test
data: LDL_diff by PSLumu
t = -0.27164, df = 8.3546, p-value = 0.7925
alternative hypothesis: true difference in means between group ari and group nasi is not equal to 0
95 percent confidence interval:
-33.46568 26.36568
sample estimates:
mean in group ari mean in group nasi
15.25 18.80
統計的には有意差ないが、とりあえずこれを箱ひげ図で比較
> Boxplot(LDL_diff ~ PSLumu, data=Dataset, id=list(method="y"))
他にCRPでの層別もしたいから、開始前CRPの分布も確認しておこう。
> with(Dataset, Hist(CRP0, scale="frequency", breaks="Sturges", col="darkgray"))
意外と開始前CRPは低めに偏っていることがわかる。要約。
> numSummary(Dataset[,"CRP0", drop=FALSE], statistics=c("mean", "sd", "IQR", "quantiles"), quantiles=c(0,.25,.5,.75,1))
mean sd IQR 0% 25% 50% 75% 100% n
1.736154 3.42934 2.46 0.01 0.04 0.1 2.5 12.49 13
CRPを2層に分離、中央値の0.1でわけます。
> Dataset$CRPsou <- with(Dataset, ifelse(CRP0 >= 0.1, "0.1以上", "0.1未満")
+ )
さてLDLの解析です。
> t.test(LDL_diff~CRPsou, alternative='two.sided', conf.level=.95, var.equal=FALSE, data=Dataset)
Welch Two Sample t-test
data: LDL_diff by CRPsou
t = -0.49162, df = 10.526, p-value = 0.6331
alternative hypothesis: true difference in means between group 0.1以上 and group 0.1未満 is not equal to 0
95 percent confidence interval:
-32.87932 20.92693
sample estimates:
mean in group 0.1以上 mean in group 0.1未満
13.85714 19.83333
箱ひげ図で比較
> Boxplot(LDL_diff ~ CRPsou, data=Dataset, id=list(method="y"))
こうみるとCRPの層でわけて解析すると平均は差があるが、分散が大きいだけで中央値は同じくらいなのだ。
CRPはPSLとの関係もあるので結果の解釈はとても複雑怪奇だ。
もう少し数が増えればさらに層を増やしてみることができるのだが。
PSL未使用でCRPが0.1未満の症例はどれだけいるだろう。
一応3症例いる。一応この3症例でLDL変化量をみてみよう。理論上はこれで炎症やPSLによる影響を排することができる。
とりあえず手動でデータセットをいじってみる。まず今までのデータセットを保存。
> save("Dataset", file="C:/Users/Owner/OneDrive/デスクトップ/Lipid paradox/Dataset20250625.RData")
もとデータのEXCELで3症例を限定して再度データセット取り込み、LDL 変化量の変数作成
それぞれ-3、16、2。ばらつきがすごい。
そして、要約。
> numSummary(nasinasi[,"LDLDif", drop=FALSE], statistics=c("mean", "sd", "IQR", "quantiles"), quantiles=c(0,.25,.5,.75,1))
mean sd IQR 0% 25% 50% 75% 100% n
5 9.848858 9.5 -3 -0.5 2 9 16 3
3症例の要約に意味があるかまったくわからない。いやまったくない。
いちおう検定。
> with(nasinasi, (t.test(LDL0, LDL3, alternative='two.sided', conf.level=.95, paired=TRUE)))
Paired t-test
data: LDL0 and LDL3
t = -0.87932, df = 2, p-value = 0.472
alternative hypothesis: true mean difference is not equal to 0
95 percent confidence interval:
-29.46592 19.46592
sample estimates:
mean difference
-5
やはり有意差なし。これはさすがに数の問題だろう。
さて、そろそろ一旦終わるか。LDLに関しては様々なことが影響しているが、少なくともIL6阻害薬使用によって、PSLが減ろうがCRPが陰性だろうが、上昇する方向であることに間違いはないだろう。そろそろ材料は揃ってきた。だんだん抄録を書いていかないと締め切りに間に合わなくなるぞ。
