親戚が集まったときの会で、ふいにそんなことを言われた。その人としては減量をがんばったということを言いたかったのだろうが、自分自身、糖質制限ダイエットに対してエビデンスを持っていないことと、1日で2kg 減るっていうのは、偶然「2日間連続で体重を測ったら2kg 減っていた」というポジティブな観測を強烈に覚えているだけの観測バイアスなので「水分だけでも1〜2kg 程度は体重の日内変動はある。1ヶ月間毎日体重測ってプロットしてから出直せ」と返してしまったのだが、本当に糖質制限ダイエットをするだけで1日で2kg 減ったらどうしよう、と思ったので、データで考えてみる。
ある程度の期間、それなりに管理された測定データが欲しいのだが、実は某氏の体重データがある。某氏はデータ取得と解析の知識があるので、データのとり方は「毎朝、起床時に自宅の体重計で測定する」というものらしい。出張などで測れないときは欠測である。170cm 半ばの某氏は、さすがにラーメンを食べ過ぎたので減量しようと思ったようである。ただし、食生活を変えたり、運動習慣を強化した、というようなダイエット行動がいまいちわからなかったのでよくわからない。
ここで、体重を測っていなかったとしても、決して体重が存在しなかった、というわけではない。観測していなかっただけである。そういう場合に、状態空間モデルを使って、「真の体重」というのもは常にあるが、それを観測するときに誤差が生じるし、場合によっては欠測にもなる、というのをstan でやる。
真の体重 は、ふたつ前までの状態の差分で決まるとする。このとき、差分に対して適当なパラメータが決まり、その分散は固定された
という実数とすると
つまり
である。
実際に観測した体重 (大文字) は、測定時に着ている服とか、体内の水分バランスとか、体重計のキャリブレーション具合とかで、なんらかの観測誤差を含むと考えられるので、適当な分散
を用いて
でサンプリングされる、とする。毎回の観測の条件を揃えると、 は小さくて済む。
欠測は、rstan ではNA を受け付けないので、欠測の日を先にindex で作成しておいて、rstan のmodel ブロック内でfor loop で回避する。
,
をcauchy 分布からサンプリングしているが、なるべくばらついて大げさな値が出て欲しかったのでやってみた。適当である。
でうまく収束したようだ。計算時間は1000 iter で1 chain 3分程度。
1日での体重減少の事後分布を見てみると、中央値で0.03 kg 程度減少する。95% 信用区間では、実は有意に体重が減少するという結果ではないようだ。しかし、200日程度のデータがあるので、分布の正の領域と負の領域が半々でない以上、負の方向(痩せる方向)へ傾いていくのは妥当だろう。
quantile(c(apply(ex$w, 1, diff)), c(a/2, 0.5, 1-a/2))
2.5% 50% 97.5% -0.08474003 -0.03705018 0.01160505
何も考えず、データ観測期間内での体重減少から日割りすると、-0.04 kg 程度だった。
diff(range(dat$weight, na.rm=TRUE))/nrow(dat)
[1] 0.04056604
さて、ここで気になるのが、「1日で2kg 痩せる」というのがどれくらいの出来事なのだろうか、ということである。体重の減少は、「真の体重減少」に加えて、そのときの偶然の誤差を加味した「体重計の指し示す値」で決まる。ここでいう偶然の誤差は、先に言った通り、体重計のキャリブレーションや、そのときの服装、直前に食事をしたり排泄をしたかどうか、が含まれる。
この「誤差」の作用は、上のモデルの によって表現される。ここから、2kg 減ったというのが95% 信頼区間(ここでは頻度主義的に) でどうか、というのを調べるには
qnorm(1-a/2, sd=s)
[1] 0.8670893
となり、0.87kg 以上の観測誤差があるとき、2kg の1日減量やばい、ということになる。よって、糖質制限ダイエットの真の1日体重減少量は1.2 kg 程度となる。
(体重の減少、が興味の対象なので、片側では? という話もあるが、このときは0.72 kg である)
単純に、糖質制限ダイエットの真の減量効果が1kg/day としよう。1ヶ月で30kg 痩せる。1ヶ月で30kg 痩せたら自分としては悪性腫瘍を疑って医療機関の受診を勧める。
というのはまあ笑い話なのだが、糖質制限はするがたんぱく質、脂質はどんだけ摂取してもよいらしい。糖質がないと脳内でエネルギーにはならないけど? と聞くと、午前中はツライらしく、脂質のエネルギーって9kcal/g だけど? と聞くと、脂質は肌の恒常性の維持に必要だからむしろたくさん摂取しろという。謎()
library(rstan) rstan_options(auto_write=TRUE) options(mc.cores=parallel::detectCores()) dat <- read.csv("weight.csv") code <- " data{ int Day; real<lower=0, upper=200> W[Day]; int Na; int na[Na]; } parameters{ real<lower=0, upper=200> w[Day]; # weight real<lower=0> b; # sigma for hidden model real<lower=0> s; # sigma for observation } model{ for(i in 3:Day){ w[i] ~ normal(2*w[i-1] - w[i-2], b); } for(i in 1:Na){ W[ na[i] ] ~ normal(w[ na[i] ], s); } b ~ cauchy(0, 0.5); s ~ cauchy(0, 0.5); } " m <- stan_model(model_code=code) W = replace(dat$weight, is.na(dat$weight), 1) standata <- list(Day=nrow(dat), W=W, Na=sum(!is.na(dat$weight)), na=which(!is.na(dat$weight))) fit <- sampling(m, standata) ex <- extract(fit) a <- 0.05 w <- apply(ex$w, 2, median) b <- median(ex$b) s <- median(ex$s) ci <- apply(ex$w, 2, quantile, c(a/2, 1-a/2)) par(mar=c(5, 4.5, 2, 2), cex.lab=1.5, cex.axis=1.5) xi <- seq(nrow(dat)) yl <- range(dat$weight, ci, na.rm=TRUE) plot(w, type="l", ylim=yl, xlab="Day", ylab="Weight [kg]") polygon(c(xi, rev(xi)), c(ci[1,], rev(ci[2,])), col=grey(0.8), border=NA) lines(w) points(dat$weight, pch=16, col=4, cex=0.8) points(which(is.na(dat$weight)), w[which(is.na(dat$weight))], pch=16, col=2) legend("topright", legend=c("実測値", "欠測時"), col=c(4, 2), pch=16, cex=3)
