No.22397 Re: 直線性の適合度検定 【青木繁伸】 2017/06/17(Sat) 15:10
もう40年も前にやっていたことで，細かいことは忘れてしまいましたね。
直線性の適合度検定も何通りかあったはずで，「この方法では，まず，物質の投与量に対する腫瘍発生率の変化の直線性について適合度検定」というのがどれか定かではないですが，たとえば，コクラン・アーミテージの検定は
http://aoki2.si.gunma-u.ac.jp/lecture/Hiritu/Armitage.html
で説明してあります。そのページにしたがってプログラムしてもよし，
http://aoki2.si.gunma-u.ac.jp/R/Cochran-Armitage.html
そこからリンクしている（R で計算してみる）ように R にすでにある prop.trend.test を使ってもよいでしょう。
http://aoki2.si.gunma-u.ac.jp/lecture/Hiritu/Armitage-r.html
http://www.tamagaki.com/math/Statistics606.html の記事とは違います（一様性の分解という意味では関係あり）。
No.22398 Re: 直線性の適合度検定 【やまが】 2017/06/17(Sat) 21:38
青木先生
ご回答いただきありがとうございます。
教えていただいたリンク先を確認し試してみます。
結果が得られ次第再度記入させていただきます。
No.22402 Re: 直線性の適合度検定+最尤法 【やまが】 2017/06/30(Fri) 23:00
青木先生
直線性の適合度検定では大変お世話になりました。
さて，現在，上述のTD50を求めるため，
腫瘍発生確率： p=1-exp(-(b0+b1x))のパラメータ（b0, b1）を推定するため最尤法を行おうとしていますが，うまくいっておりません。元々，
http://www.yukms.com/biostat/haga/download/archive/likelihood/Likelihood.pdf
を参照し，Excelでは計算できていたのですが，
Rでも計算できるように取り組んでおります。
http://www.agr.nagoya-u.ac.jp/~seitai/document/R2009/t090614glmIntro.pdf
のロジスティック回帰を参考に，p=1-exp(-(a + b・x))とすればglmで同様に計算できるのではないかと考え，以下を試してみました。
Ln(1/(1-p))=a + b・x> x <-c(0,1,3,9)しかしながら，正しくない値が得られてしまいました。
> r <-c(0,1,3,38)
> n <-c(49,50,50,50)
> p <- r/n
> plot(x,p)
>
> tumor<-cbind(n,n-r)
> glm.tumor<-glm(tumor~1+x,family=binomial)
>
> glm.tumor$coefficients["x"]
x
0.1539377
> glm.tumor$coefficients["(Intercept)"]
(Intercept)
-0.1495911
>
> x <- seq(min(x),max(x),by=0.01)
> eta.pred<- glm.tumor$coefficients["x"]*x+glm.tumor$coefficients["(Intercept)"]
> p.pred <- 1-exp(-eta.pred)
>
> lines(x, p.pred)
そこで，対数尤度関数を最大化するパラメータを求めるために，
http://abrahamcow.hatenablog.com/entry/2016/04/09/085412
http://www.agr.nagoya-u.ac.jp/~seitai/document/R2009/t090614glmIntro.pdf
を参考にoptimを用いて以下を試みました。> x <-c(0,1,3,9)今度は正しい値が得られたのですが，
> r <-c(0,1,3,38)
> n <-c(49,50,50,50)
>
> LL<-function(par){
+ beta0<-par[1]
+ beta1<-par[2]
+ eta <- beta0+beta1*x
+ p<- 1-exp(-eta)
+ NLL<- sum(dbinom(r,n,p,log=TRUE))}
>
> opt1<-optim(c(1,1),LL,control=list(fnscale=-1))
There were 50 or more warnings (use warnings() to see the first 50)
> #opt1<-optim(c(1,1),LL,control=list(fnscale=-1),hessian=TRUE)
> opt1$par
[1] 5.405072e-09 9.328301e-02
> opt2<-optim(opt1$par,LL,control=list(fnscale=-1),method="BFGS",hessian=TRUE)
Error in optim(opt1$par, LL, control = list(fnscale = -1), method = "BFGS", :
non-finite finite-difference value [1]
In addition: Warning message:
In dbinom(r, n, p, log = TRUE) : NaNs produced
> my.std<-sqrt(-diag(solve(opt1$hessian)))
Error in array(x, c(length(x), 1L), if (!is.null(names(x))) list(names(x), :
'data' must be of a vector type, was 'NULL'
> up<-qnorm(0.995,opt1$par,my.std)
Error in qnorm(0.995, opt1$par, my.std) : object 'my.std' not found
> low<-qnorm(0.005,opt1$par,my.std)
Error in qnorm(0.005, opt1$par, my.std) : object 'my.std' not found
> round(matrix(c(low,up),2,2),3)
Error in matrix(c(low, up), 2, 2) : object 'low' not found
完全に収束していないのかhessianが得られませんでした。
（ちなみに以下のようにしてもだめでした）opt1<-optim(c(1,1),LL,control=list(fnscale=-1),hessian=TRUE)大変申し訳ございません。
何らかの解決法がございましたらご教示いただけないでしょうか？
よろしくお願い申し上げます。
No.22403 Re: 直線性の適合度検定 【青木繁伸】 2017/07/01(Sat) 10:10
いくつか間違いがありますねtumor <- cbind(n, n-r)ではなく，tumor <- cbind(r，n)です。しかし，glm(tumor~1+x, family=binomial)で得られるのは，「ロジスティック回帰」の結果なので，eta.pred<- glm.tumor$coefficients["x"]*x+glm.tumor$coefficients["(Intercept)"]で予測値は計算できません。
p.pred <- 1-exp(-eta.pred) # ロジスティック曲線の式ではない
ロジスティック回帰の結果を示すには，y = 1/(1+exp(-predict(glm.tumor)))とすれば，「ロジスティック回帰」はうまくいっていることがわかるでしょう。
points(x, y, col=2, pch=19)
x2 <- seq(min(x),max(x),by=0.01)
eta.pred<- glm.tumor$coefficients["x"]*x2+glm.tumor$coefficients["(Intercept)"]
p.pred <- 1/(1+exp(-eta.pred))
lines(x2, p.pred)
No.22404 Re: 直線性の適合度検定 【青木繁伸】 2017/07/01(Sat) 15:50
optim を使う場合，対数尤度は
http://www.ma.utexas.edu/users/mks/RA/1hit.html
によれば，LL2 = function(par) {のようで，sum(dbinom(r, n, p, log = TRUE)) とはちがうのだけど，得られる答えは同じになる。hessian = TRUE を指定すると
p = 1-exp(-par[1]-par[2]* x)
valid = (p != 0) & (p != 1)
p = p[valid]
n = n[valid]
r = r[valid]
sum(r*log(p) + (n-r)*log(1-p))
}optim(c(1, 1), LL, control = list(fnscale = -1), hessian = TRUE)解が得られない。
LL2 を使ったans = optim(c(1, 1), LL2 , control= list(fnscale=-1), hessian=TRUE)では，エラーメッセージ付きで hessian は求まるが，1×1 である。
ans
なお，1-exp(-par[1]-par[2]* x) は原点を通らないので，1-exp(-par[1]* x) とすると，LL3 = function(q) {として，optimize を使うのがよさそうではある。
p = 1-exp(-q*x)
valid = (p != 0) & (p != 1)
p = p[valid]
n = n[valid]
r = r[valid]
sum(r*log(p) + (n-r)*log(1-p))
}ans2 = optimize(LL3, lower=0, upper=1, maximum=TRUE)
ans2
No.22405 Re: 直線性の適合度検定 【やまが】 2017/07/02(Sun) 21:44
青木先生
丁寧なご回答をいただきありがとうございます。
大変助かりました。
確かに，Log(0)はよくないと気づかされました。
次は，このパラメータの信頼区間を求める必要があるのですが，
もう少し自分で考えてみようと思います。
結局再度質問させていただくことになるかもしれませんが，
今後ともよろしくお願いいたします。
No.22410 Re: 直線性の適合度検定 【やまが】 2017/07/07(Fri) 22:23
いつもお世話になっております。
先日教えていただいた方法をもとに，
以下の資料のp.3の方法を参考に信頼区間を求めました。
http://www.agr.nagoya-u.ac.jp/~seitai/document/R2009/t090614glmIntro.pdfx <-c(0,1,3,9)一応うまくいっているようなのですが，
r <-c(0,1,3,38)
n <-c(49,50,50,50)
LL3 <- function(q) {
p <- 1-exp(-q*x)
valid = (p != 0) & (p != 1)
p = p[valid]
n = n[valid]
r = r[valid]
sum(r*log(p) + (n-r)*log(1-p))
}
ans2 <- optimize(LL3, lower=0, upper=1, maximum=TRUE)
ans2
q<-ans2$maximum
n.sim <- 10000
beta <- data.frame(beta1 = rep(NA, n.sim))
eta <- q*x
p <- 1 - exp(-eta)
for(i in 1:n.sim){
r <- rbinom(length(p), n, p)
ans22 <- optimize(LL3, lower=0, upper=1, maximum=TRUE)
ans22$maximum
beta[i, ] <- ans22$maximum
}
beta.ci <- apply(beta, 2, quantile, prob = c(0.005, 0.995))
beta.ci
log(2)/ans2$maximum
log(2)/beta.ci[1,1]
log(2)/beta.ci[2,1]
この方法だとどうしても計算のたびに異なる結果になります。
そこで，以下の資料のp.36,37を参考に求められないか検討しております。
http://www.yukms.com/biostat/haga/download/archive/likelihood/Likelihood.pdf
しかしながら，p.36に，
”＃最尤推定値π^の対数尤度と帰無仮説πの対数尤度との差の2 倍は，近似的に，自由度1のカイ2 乗分布に従う”
という記述があるのですが，なぜそうなるのか基本的な考え方がわかっておりません。
もう少し詳しく解説されている書式やwebページ，または考え方のヒントがございましたらご教示いただけないでしょうか？
また，この方法以外に信頼区間を求める良い方法がございましたら併せてご教示いただけないでしょうか？
お手数をおかけして申し訳ございません。
よろしくお願い申し上げます。
● 「統計学関連なんでもあり」の過去ログ--- 048 の目次へジャンプ
● 「統計学関連なんでもあり」の目次へジャンプ
● 直前のページへ戻る