保健統計学特講
2024.5.14
交絡と回帰モデル●
症例数とパワー計算：保健学共通特講 IV のサンプルサイズの設計の話の裏返し●
一般にサンプルサイズが大きいほど検出力（ Statisticalpower ）は大きくなる●
検出力は 1 から第二種の過誤（本当は差があるのにそれを見逃す確率）を引いた値●
多重検定●
３群以上の比較をするとき，２群間の比較を繰り返すと全体としての第一種の過誤（本当は差がないのに誤って差があると判定してしまう確率）が大きくなってしまうので，個々の検定の有意水準あるいは有意確率を調整する（これを検定の多重性の調整という）必要がある。多重比較ともいう●
ボンフェローニの不等式に基づくボンフェローニの方法が最も単純。ホルムの方法や FDR が標準的。
●

交絡 (confounding) とは●●●要因曝露→結果という因果関係を考えているとき，要因とも結果とも関連があって，しかも中間媒介因子ではない因子があるとき，これを交絡因子と呼ぶ。交絡因子の影響によって，注目している因果関係は歪められる。この歪みを交絡と呼ぶ→交絡を制御して解析しないと、要因曝露と結果の関係を正しく評価できない
交絡要因の必要条件（例：肥満⇒高血圧に対する年齢）●
結果（高血圧）に影響している（年齢が高いほど高血圧は増える）●
要因曝露に偏りがある（肥満群と非肥満群で年齢分布が異なっている）●
中間媒介因子ではない（肥満が高齢をもたらすのではない）：中間媒介因子を介する影響を取り除いてしまうと、要因曝露がもたらす結果への効果を過小評価することになるので、中間媒介因子は交絡扱いできない
要因曝露に偏りがあることが交絡因子の必要条件だからといって、検定によって分布に有意差があるという理由だけで、変数を交絡として制御するのは誤り●
RCT では割り付けた変数以外の変数の分布には差がないようにランダム割り付けを行っているので、差がないことは既知。いくら p 値が小さくても、それは差があることを意味しない。第一種の過誤そのもの。
●
観察研究でも、検定によって要因曝露の有無間で分布に有意差がある変数があったとしても、その変数が結果に影響しているかどうかを考慮していないので交絡の必要条件を満たしていない。

DAG の利用
●

●

●

●

●

DAG （ Directed Acyclic Graph ＝有向非巡回グラフ）は、因果推論をする際に、制御すべき交絡因子を見つけるのに役立つ参考文献としてお薦めしたい本●
佐藤俊哉『宇宙怪人しまりす 統計よりも大事なことを学ぶ』朝倉書店、 2024●
林岳彦『はじめての統計的因果推論』岩波書店、 2024観測された変数を四角で囲み、それらの間を矢印の元が原因となる変数、矢印の先が結果となる変数であるように矢印でつなぎ、矢印の向きを辿っていっても元の変数に戻ってこないようなグラフを DAG と呼ぶ
複数の観測された変数について想定される因果関係に基づいた DAGを描くと、「裏口テスト」と呼ばれる手法を使って、制御すべき交絡変数を見出すことができる（注）その研究では観測されていない変数や、原理的に観測不可能な要因も、結果となる変数に影響を与えている可能性も忘れてはならない。

裏口テストの例（佐藤 2024 、 pp.57-63 ）年齢
かぜの重症度

年齢
かぜの重症度ヨクナール使用かぜ症状の回復

ヨクナール使用かぜ症状の回復

解熱作用

解熱作用

想定される因果関係の図（ DAG になっている）

ヨクナールの効果を取り除いてもヨクナール使用とかぜ症状の回復に関連がある（＝裏口道がある）かを調べる→２つの裏口道がある※ 裏口道がすべてなくなるまで、通り道にある変数を除去する。
ただし２つの変数から１つの変数に矢印が入る合流点があると、そこで道は塞がるので裏口道ではなくなることと、除去によって合流が消えて新たな裏口道が生まれる場合がある点に注意

年齢
かぜの重症度ヨクナール使用年齢
かぜの重症度

かぜ症状の回復

ヨクナール使用かぜ症状の回復

解熱作用

解熱作用

年齢を調整する→裏口道が残るので調整不十分（→この場合、この調整は不要）

かぜの重症度を調整する→裏口道が消える＝交絡としてかぜの重症度を調整すればヨクナール使用の効果を調べることができる

交絡の例 (1)

●
●

●

●

スウェーデンの女性とカザフスタンの女性の 1992 年の死亡率
スウェーデンの方がカザフスタンより衛生水準・医療水準・生活水準が高い→年齢別死亡率は全年齢層でスウェーデン＜カザフスタン→当然，平均寿命はスウェーデン＞カザフスタン→が，粗死亡率はスウェーデン (10.5) ＞カザフスタン (6.3)
理由は以下 2 点●
どちらの国でも高齢者の死亡率＞＞若者や成人の死亡率●
高齢者の割合はスウェーデン＞カザフスタン
年齢構造が交絡因子となっている

交絡の例 (2)●

年齢が交絡因子となっているために，喫煙者の方が非喫煙者よりも死亡率が低くなることがある（出典：『ロスマンの疫学』）

交絡の例 (3)●

ダウン症ハイリスクなのは出生順位が後の子か？ それとも母の年齢が高い子か？ （出典『ロスマンの疫学』）

交絡の制御デザインによる制御●
性別や年齢による交絡がありそうな場合は，それらで層別した層別無作為抽出をする，コホート研究なら曝露群と非曝露群を設定するときに年齢や性別をマッチングするか，特定の年齢や性別に限定化して研究すると，それらの影響がデータに表れない●解析方法の工夫による制御●
層別解析：層別に比較する●
標準化：年齢調整死亡率とか●プール化：どの層でも同じ関連や違いがあると仮定して要約統計量を計算する（マンテルヘンツェルの要約オッズ比など）●
多変量解析：交絡因子も独立変数に含めたロジスティック回帰分
析などの一般化線型モデル，傾向スコアを用いた共変量の調整， DID (Difference in Difference) 法など– Rothman は、独立変数に含めることによって、他の独立変数の係数が大きく変わるかどうかが、交絡因子かどうかを判定する目安になるとしている●

ロジスティック回帰分析の例●●

エコポイントチェックデータ
https://minato.sip21c.org/advanced-statistics/ecopx.txt
https://minato.sip21c.org/humeco/ecopoint.html●高月紘（編著）「自分の暮らしがわかるエコロジー・テスト：環境問題は生活のエコ度チェックから」（講談社ブルーバックス）●
CGI 化→表にチェックして，「テスト実行」– エコポイント総合点が 99.8 点満点（丸め誤差で≠ 100 ）– 温暖化問題，廃棄物問題，水質汚染問題，大気汚染問題，有害化学物
質問題のエコポイントがそれぞれ 100 点満点– 若い層を中心とした対象者 356 人の平均が 42.6
環境問題の講演を聞きに来た人たち 182 人の平均が 48.1

質問項目

ロジスティック回帰 (1)Q4 と Q10 の関係には年齢と性別の影響が両方ある（年齢と性別は Q4 と Q10 の関係を交絡している）●
年齢と性別と Q10 を独立変数， Q4 を従属変数としてロジスティック回帰分析●
年齢と性別の影響を調整して，●
Q10 が TRUE の人（いつも / 大体 / 時々公共交通利用＝公共交通利用者）は●
Q10 が FALSE の人（公共交通利用がたまに / 皆無＝マイカー族）に比べて，●どれくらい Q4 が TRUE である（いつも / 大体 / 時々買い物袋持参＝エコバッグ派である）割合が高いかわかる●ロジスティック回帰 (2)●Q4 が TRUE である割合を p とすると，log(p/(1-p))=β0+β1AGE+β2SEX+β3Q10
TRUE を 1 ， FALSE を 0 とすると，他の変数を調整した上での Q10 の効果
は，●
Q10 が TRUE の場合の Q4 が TRUE である割合を p1 とすれば，log(p1/(1-p1))=β0+β1AGE+β2SEX+β3●
Q10 が FALSE の場合の Q4 が TRUE である割合を p0 とすれば，log(p0/(1-p0))=β0+β1AGE+β2SEX●
辺々引けば log(p1/(1-p1)) - log(p0/(1-p0))=β3●
結局この左辺は， log[(p1/(1-p1))/(p0/(1-p0))] となる。これは対数オッズ比である。つまり， AGE と SEX の影響を調整した上で Q10 が TRUE である
人は FALSE の人に比べて， exp(β3) 倍， Q4 が TRUE になりやすい●「性別と年齢の影響を調整した上で，公共交通利用者はマイカー族に比べて，エコバッグ派である可能性が exp(β3) 倍大きい」ということである。 β0～ β3 は，一般化線型モデルの係数として推定される。

ロジスティック回帰 (3)●●R のコードx <- read.delim("https://minato.sip21c.org/advanced-statistics/ecopx.txt") # 読み込み
x$Q04c <- ifelse(x$Q04>=2.5, 1, 0) # スコアなので 1/0 に変換
x$Q10c <- ifelse(x$Q10>=4.8, 1, 0) # スコアなので 1/0 に変換
res <- glm(Q04c ~ Q10c+AGE+SEX, family=binomial(logit), data=x) # ロジスティック回帰
summary(res) # 回帰分析の結果を出力
exp(coef(res)) # 対数オッズ比なので係数の exp() をとってオッズ比を得る
exp(confint(res)) # 信頼区間の exp() をとってオッズ比の 95% 信頼区間を得る
library(fmsb) # fmsb パッケージを読み込む
NagelkerkeR2(res) # モデルの説明力を示す Nagelkerke の R2 の計算
library(stargazer) # きれいに出力するための stargazer パッケージを読み込む
sink("logistic-res.tex") # このファイルを TeXLive でコンパイルすると右の表生成
cat("\\documentclass[a4paper,12pt]{article}\n\\begin{document}\n")stargazer(res, ci=TRUE, apply.coef=exp, type="latex",
notes=sprintf("Nagelkerke's R-squared: %4.2f", NagelkerkeR2(res)$R2))
cat("\\end{document}\n")sink()
出力結果 Estimate Std. Error z value Pr(>|z|)(Intercept) -0.5479 0.3232-1.695 0.09006
Q10c1.3036 0.2886
4.517 6.28e-06 ← 6.28×10-6 の意味
AGE-0.1864 0.1106-1.685 0.09202
SEX0.7565 0.2751
2.750 0.00596
AIC: 391.21(Intercept)
Q10cAGE
SEX
0.5781719 3.6826760 0.8299727 2.1308305
2.5 %
97.5 %(Intercept) 0.3032334 1.081392
Q10c2.1093893 6.560746
AGE0.6667677 1.030981
SEX1.2550376 3.700273$R2 [1] 0.1333421 ← モデルの説明力は 13% 程度

jamovi でやるとロジスティック回帰 (4)●●最適カットオフ 0.639 ，
感度 0.576 ，特異度 0.774
AUC は 0.687（あまり良くない）

1.00.80.6

Sens: 57.6%
Spec: 77.4%
PV+: 44.8%
PV-: 20.9%

0.4●

y = 0.639

0.2●

y <- predict(res,data.frame(Q10c=x$Q10c,AGE=x$AGE, SEX=x$SEX),"response")
library(Epi)
ROC(y, x$Q04c)

Variableest. (s.e.)(Intercept) -2.225 (0.501)
test 4.426 (0.823)
Model: x$Q04c ~ y0.0●

得られたロジスティック回帰式
を使って Q04c を予測するための ROC 分析

Sensitivity●

Area under the curve: 0.687

0.0

0.2

0.4

0.6

1-Specificity

0.8

1.0混合効果モデル●●

アウトカムとなる健康事象が個人レベルで，個人レベルの要因と集団レベルの（サブグループ）要因があり，個人レベルでの回帰モデルの切片と傾きの片方または両方がサブグループによって違う場合，混合効果モデルを当てはめて交絡を調整する
例えば， Demographic and Health Surveys の 2017 年インドネシアのデータ（ IDIR71FL.SAV という SPSS 形式のデータ）には V012 （年齢）， V025 （都市部／農村部）， V106 （教育歴）， V201 （既往出生児数）， V361 （避妊利用）が含まれている。既往出
生児数を他の変数で説明するモデルを考える。
●
●
●
●
●
●

●

●

library(foreign)x <- read.spss("IDIR71FL.SAV")lm(V201 ~ V012+V025+V106+V361, data=x) # 重回帰
glm(V201 ~ V012+V025+V106+V361, data=x, family="poisson") # ポアソン回帰
library(lmerTest)lmer(V201 ~ V012 + V106 + V361 + (1 | V025), data=x, REML=FALSE) # 混合

この結果はどれも大差ない。係数の信頼区間を得るには confint(coef()) を使う。年齢の係数は正，教育歴なしに比べどの教育歴も負，避妊利用は現在利用中に比べたのカテゴリは負，都市部に比べ農村部は正の係数となっている
混合効果モデルにはグループレベルの変数により切片だけが異なるランダム切片モデル（上のコードはそれを仮定している）と，傾きも異なると仮定するモデルがある

結果の比較（重回帰とポアソン回帰）
Call:lm(formula = V201 ~ V012 + V025 + V106 + V361, data = x)
Residuals:
Min
1Q Median
3Q
Max-4.2150 -0.6136 -0.0643 0.3418 10.9900

glm(formula = V201 ~ V012 + V025 + V106 + V361, family = "poisson", data = x)
Deviance Residuals:
Min
1Q Median
3Q
Max-3.4452 -0.7446 -0.3526 0.3175 6.5011
Coefficients:

Estimate Std. Error z value Pr(>|z|)Coefficients:(Intercept)-0.6740534 0.0281197 -23.97 <2e-16 ***Estimate Std. Error t value Pr(>|t|)V012
0.0522313 0.0004569 114.33 <2e-16 ***(Intercept)
0.6291598 0.0482432 13.041 < 2e-16 *** V025Rural
0.1268545 0.0070699 17.94 <2e-16 ***V012
0.0762938 0.0006318 120.759 < 2e-16 *** V106Primary-0.2814560 0.0194853 -14.45 <2e-16 ***V025Rural
0.1844247 0.0105574 17.469 < 2e-16 *** V106Secondary-0.4031578 0.0196109 -20.56 <2e-16 ***V106Primary-0.6379498 0.0394769 -16.160 < 2e-16 *** V106Higher-0.5929119 0.0216679 -27.36 <2e-16 ***V106Secondary-0.9298056 0.0392179 -23.709 < 2e-16 *** V361Used since last birth -0.1458802 0.0102629 -14.21 <2e-16 ***V106Higher-1.2269346 0.0403892 -30.378 < 2e-16 *** V361Used before last birth -0.1405757 0.0100127 -14.04 <2e-16 ***V361Used since last birth -0.3381545 0.0173691 -19.469 < 2e-16 *** V361Never used
-1.2205031 0.0122418 -99.70 <2e-16 ***V361Used before last birth -0.1027063 0.0188352 -5.453 4.98e-08 *** --V361Never used
-1.1855481 0.0133170 -89.025 < 2e-16 *** Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1--Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1(Dispersion parameter for poisson family taken to be 1)
Residual standard error: 1.138 on 49618 degrees of freedom
Multiple R-squared: 0.5461,
Adjusted R-squared: 0.546
F-statistic: 7461 on 8 and 49618 DF, p-value: < 2.2e-16

Null deviance: 90239 on 49626 degrees of freedom
Residual deviance: 37663 on 49618 degrees of freedom
AIC: 129952
Number of Fisher Scoring iterations: 5

線形混合モデルの結果
Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's method[lmerModLmerTest]
Formula: V201 ~ V012 + V106 + V361 + (1 | V025)
AIC
BIC logLik deviance df.resid
153677.8 153765.9 -76828.9 153657.8 49617Scaled residuals:
Min
1Q Median
3QMax-3.7054 -0.5386 -0.0568 0.2999 9.6585
Random effects:
Groups Name
Variance Std.Dev.
V025(Intercept) 0.008444 0.09189
Residual
1.294560 1.13779
Number of obs: 49627, groups: V025, 2
Fixed effects:
Estimate Std. Error
df t value Pr(>|t|)(Intercept)
7.220e-01 8.041e-02 4.630e+00 8.978 0.000423 ***V012
7.629e-02 6.317e-04 4.963e+04 120.761 < 2e-16 ***V106Primary-6.382e-01 3.947e-02 4.963e+04 -16.166 < 2e-16 ***V106Secondary-9.303e-01 3.921e-02 4.963e+04 -23.724 < 2e-16 ***V106Higher-1.228e+00 4.038e-02 4.962e+04 -30.396 < 2e-16 ***V361Used since last birth -3.382e-01 1.737e-02 4.963e+04 -19.471 < 2e-16 ***V361Used before last birth -1.027e-01 1.883e-02 4.963e+04 -5.455 4.92e-08 ***V361Never used
-1.186e+00 1.332e-02 4.963e+04 -89.038 < 2e-16 ***--Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Correlation of Fixed Effects:(Intr) V012 V106Pr V106Sc V106Hg V36slb V36blb
V012-0.340
V106Primary -0.485 0.077
V106Secndry -0.517 0.165 0.943
V106Higher -0.491 0.133 0.912 0.939
V361Usdsnlb -0.059 0.030 0.013 0.013 0.013
V361Usdbflb -0.014 -0.137 0.022 0.023 0.026 0.184
V361Nevrusd -0.229 0.465 0.101 0.094 0.044 0.281 0.176

jamovi の場合
線形回帰
モデル適合度指標
モデルR1

R²0.7390.546

モデル係数 - V201
予測変数

推定値

標準誤差tp

切片 ᵃV0120.62920.07630.0482
6.32E-04

13.04 < .001
120.76 < .001

V025:
Rural – Urban

0.18440.010617.47 < .001

V106:
Primary – No education
Secondary – No education
Higher – No education-0.6379-0.9298-1.2269

0.03950.03920.0404-16.16 < .001
-23.71 < .001
-30.38 < .001

V361:
Used since last birth – Currently using
Used before last birth – Currently using
Never used – Currently using-0.3382-0.1027-1.1855

0.01740.01880.0133-19.47 < .001
-5.45 < .001
-89.03 < .001

jamovi によるポアソン混合効果モデル

DID （時間経過が交絡である場合）
Snow J (1855) が分析したロンドンのコレラのデータ●
1854 年の感染地図と自然実
験（右の 2 社の両方から給水を受けている地域で，水の塩分を調べて，どちらの会社からの給水か特定し， S&V 社からの給水を受けている群で死亡リスクが数倍あることを示した）は有名●
1849 年には両社ともテムズ川下流から採水していたが， 1852 年に L 社は上流に採水地点を移した●

水道会社

1849 年コレ
ラ死者数

1854 年コレラ死者数

S&V

22612458

L

16237

両方から39052547

採水を変えなかった（非介入群） S&V 社から給水を受けている群は死者やや増，採水を変えた（介入群） L 社から給水を受けている群は死者減。ただしその減少が採水点を変えた効果かどうかは，時間経過があるので不明●
差の差を取ると 43 ％減●(37-162)/162-(2458-2261)/
2261=-0.43

症例数とパワー計算 (1)「保健学研究共通特講 IV 」テキスト第 4 章を参照
https://minato.sip21c.org/ebhc-text.pdf●
例えば独立 2 群の平均値の差の t 検定をする場合●
標準偏差，意味のある差， α エラー， β エラーを決めればサンプルサイズを計算できる●⇔サンプルサイズ，標準偏差，実際の差， α エラーがわかればβ エラーがわかる→ 1 － β エラーとして検出力もパワー計算とは検出力を求めること。
●
他の条件が同じなら，サンプルサイズが大きいほど検出力は上がる●

症例数とパワー計算 (2) 計算例
独立 2 群間の両側 t 検定をして，サンプルサイズが各群 100 人，平均値の差が 1.5 ，標準偏差が 4で，有意水準が 0.05 の場合，●
PS によれば（右図）検出力は0.752 。 G*Power は次スライド●R ではpower.t.test(n=100, delta=1.5,sd=4, sig.level=0.05,type="two.sample",alternative="two.sided")と打てば， power=0.7513...
●
EZR では 0.755 となる●
小数点以下３桁目が微妙に違うが気にしなくて良い●G*Power の場合

効果量について
●

●

●

●

統計的仮説検定では，「差がない」とか「関連がない」という帰無仮説の下で，解析対象のデータよりも「差がある」「関連がある」側に外れたデータが偶然得られる確率に相当する p 値が「有意水準」より低いとき，偶然でないと見なして帰無仮説を棄却する●
サンプルサイズを大きくすれば，「偶然でない」差や関連を検出できる可能性（～検出力）は上がるが，差や関連が「ゼロではない」ことの確信度があがるだけ●
p 値からでは，差や関連の大きさは不明
しかし，研究者が本当に知りたいのは，処置やグループの違いが結果にどれくらい大きな差をもたらすのか，あるいは二つの変数の関連はどれくらい大きいのか，という点。その指標が効果量
最近は信頼区間に加えて効果量も示すべきとする Journal が多い●
G*Power では効果量も計算できるし、必要な効果量からサンプルサイズを設計することもできる
参考文献としては，大久保街亜，岡田謙介『伝えるための心理統計：効果量・信頼区間・検定力』勁草書房， ISBN978-4-326-25072-1 がお薦め

2 系統の効果量
効果量とは，帰無仮説が正しくない程度を量的に表す指標●
データに適用する統計モデルにより 2 系統に分かれる●
独立変数と従属変数の区別を含む場合（独立 2 群間の平均値の差の検定も，群分け変数を独立変数，量の変数を従属変数と考える線型モデルで表せるので，この枠組みに含まれる）→独立変数が従属変数に対して影響を及ぼす程度が効果量→ d 族の効果量（ Cohen の d や Hedges の g など）＊ R では effsize パッケージの cohen.d() 関数で計算できる●独立変数と従属変数の区別がない場合→変数間の関連の大きさを表すのが効果量→ r 族の効果量（ピアソンの相関係数，決定係数，イータ 2 乗，オメガ 2 乗など）
＊ R では lsr パッケージの etaSquared() 関数でイータ 2 乗を得る●
一般に効果量は標準化された値なので元データの単位によらない●帰無仮説が完全に正しいとき効果量は 0●

多重検定 (1)
●
●

●

●

●

●

「保健学共通特講 IV 」テキスト第９章を参照
多重比較の方法は多数あるが，有名なのは，ボンフェローニ (Bonferroni) の方法，ホルム (Holm) の方法，シェフェ (Scheffé) の方法，テューキー (Tukey) のHSD ，ダネット (Dunnett) の方法，ウィリアムズ (Williams) の方法である。最近は， FDR(False Discovery Rate) 法も良く使われるようになった。
ボンフェローニの方法とシェフェの方法は検出力が悪いので，特別な場合を除いては使わない方が良い平均値の差の多重比較ならテューキーの HSD （ R では TukeyHSD() 関数を使う），その他の場合はホルムの方法か FDR 法が薦められる
R の pairwise.t.test() や pairwise.prop.test() や fmsb パッケージのpairwise.fisher.test() などでは， p.adjust.method="fdr" のようなオプションとして検定の多重性の調整方法を指定する
等分散を仮定しない方法として Games-Howell 法がある。 R では rosetta パッ
ケージの posthocTGH() 関数で可能である。
●
組み込みデータ chickwts で，餌の違いによる鶏の体重を多重比較するには，rosetta パッケージを呼び出した後，posthocTGH(chickwts$weight, chickwts$feed, method="games-howell")多重検定 (2) EZR での実行例
●

●

●

●

●

R の組み込みデータ（６週間異なる餌で飼育した後のヒヨコ体重 chickwts ：変数 feed が餌，変数 weight が体重）を分析
一元配置分散分析で feed の weight への効果は有意
EZR ではファイル＞パッケージに含まれるデー
タを読み込む＞ datasets の chickwts を選択
し，統計解析＞連続変数の解析＞３群以上の間の平均値の比較と選び，右上画面でオプション選択
Tukey の場合右下の図が表示される→中央の縦棒が「差がゼロ」の線。 95% 同時
信頼区間がそこに掛かっていない組合せが 5%水準で「有意差あり」# R では Tukey は
res <- TukeyHSD(aov(weight~feed,data=chickwts))
plot(res)print(res) # 多重性調整済 p 値表示# 他は（調整法は "" 内に指定）
pairwise.t.test(chickwts$weight,chickwts$feed, p.adjust.method="fdr")# 行列形式で２つずつの餌の組合せについて
検定の多重性調整済 p 値が表示される

95% family-wise confidence level
horsebean-casein
linseed-casein
meatmeal-caseinsoybean-casein
sunflower-casein
linseed-horsebean
meatmeal-horsebean
soybean-horsebean
sunflower-horsebean
meatmeal-linseedsoybean-linseed
sunflower-linseedsoybean-meatmeal
sunflower-meatmeal
sunflower-soybean
-200

-100

0100200