myscratchpad

フリー統計ソフト EZR on Rコマンダー (β版)
Rコマンダーのカスタマイズ機能を利用して、医療統計でしばしば用いる以下のような統計解析機能を組み込んだ統計ソフト EZR

http://www.jichi.ac.jp/saitama-sct/SaitamaHP.files/statmed.html#_pageHead 

今はもうほとんどこっちを更新してます。

以下はこのブログへのアクセスに使われた検索ワード
Excelのエラーバーの検索が多いｗ

1     R
2     r
3     エラーバー
4     R (Yahoo検索)
5     ロジスティック
6     spss
7     エクセル
8     棒グラフ
9    2007
10     excel
11     グラフ
12     交互作用
13     効果量
14     多重比較
15     信頼区間
16     lavaan
17     lrm
18     ロジスティック回帰
19     因子分析
20     最尤推定

広告対策 

Quick-Rを真似しながらまとめてみる。ブログはブログで思いついたことをメモしていく
http://eau.uijin.com/

# この記事と同じくoptim関数で
# 井関先生のサイトを参考にした。感謝。

# Bill Prinzmetalのサイトにある論文とExcelファイルを参考
# # Dodson, C. S., Prinzmetal, W., & Shimamura, A. (1998). Using Excel estimate parameters from observed data: An exa mple from source memory data. Behavioral Research Methods, Instruments & Computers, 30, 517-526.


# 一番単純な例でやってみる

# 初期値
values <- rep(0.5, 7)
values

# 観測値
obs <- c(612,151,77,123,643,74,19,18,383)
# 観測回答率
act.prb <- c(
# Sorce=Eric
obs[1]/sum(obs[1:3]), # response=eric
obs[2]/sum(obs[1:3]), # response=julie
obs[3]/sum(obs[1:3]), # response=new
# Sorce=Julie
obs[4]/sum(obs[4:6]), # response=eric
obs[5]/sum(obs[4:6]), # response=julie
obs[6]/sum(obs[4:6]), # response=new
# Sorce=New
obs[7]/sum(obs[7:9]), # response=eric
obs[8]/sum(obs[7:9]), # response=julie
obs[9]/sum(obs[7:9])  # response=new
)

# データの確認
data.frame(obs, act.prb)

# 関数定義
prd.prb <- vector("numeric", 9)
fnc <- function(x, ifm=F){
d1 <- x[1]; d2 <- x[1] # 等値制約のため同じ値を入れる
ld1 <- x[3]; ld2 <- x[3] # 等値制約のため同じ値を入れる
a <- x[5]; g <- x[5] # 等値制約のため同じ値を入れる
b <- x[7]

prd.prb[1] <- (d1*ld1) + (d1*(1-ld1)*a) + ((1-d1)*b*g)
prd.prb[2] <- (d1*(1-ld1)*(1-a)) + ((1-d1)*b*(1-g))
prd.prb[3] <- (1-d1) * (1-b)
prd.prb[4] <- (d2*(1-ld2)*a) + ((1-d2)*b*g)
prd.prb[5] <- (d2*ld2) +(d2*(1-ld2)*(1-a)) + ((1-d2)*b*(1-g))
prd.prb[6] <- (1-d2) * (1-b)
prd.prb[7] <- b * g
prd.prb[8] <- b * (1-g)
prd.prb[9] <- 1 - b

x1 <- 2*obs
x2 <- log(act.prb/prd.prb)
g2v <- x1*x2
estv <- c(d1, d2, ld1, ld2, a, g, b)
reslist <- list(prd.prb, estv, g2v, sum(g2v))
names(reslist) <- c("prd.prb", "estimates", "G2", "sumG2")
if(ifm==T) return(reslist) else return(sum(g2v))
}

(opres <- optim(par=values, fn=fnc, method="BFGS"))
est.v <- opres$par

# 改めて予測値を出す
x <- est.v
d1 <- x[1]; d2 <- x[1] # 等値制約のため同じ値を入れる
ld1 <- x[3]; ld2 <- x[3] # 等値制約のため同じ値を入れる
a <- x[5]; g <- x[5] # 等値制約のため同じ値を入れる
b <- x[7]
est.v <- c(d1,d2, ld1, ld2, a, g, b)

fnc(est.v, ifm=T)
# 適合度の評価
dfv <- 2 # 自由度は2
g2 <- fnc(est.v, ifm=T)$sumG2
pchisq(g2, dfv)

# 前の記事から進歩しておらず、相変わらずconstrOptim関数の使い方がわからない。たぶんこれで制約を色々入れないとダメだろうな…
# Multinomial Processing Treeモデル専用のパッケージmptというのがあるようだ。よくこんなマニアックなものを…。これも使い方がよくわからんのう
http://homepages.uni-tuebingen.de/florian.wickelmaier/wickelmaier_software.html#mpt


# 参考文献
Batchelder, W. H., & Riefer, D. M. (1990). Multinomial processing models of source monitoring. Psychological Review, 97, 548-564.
Batchelder, W. H., & Riefer, D. M. (1999). Theoretical and empirical review of multinomial process tree modeling. Psychonomic Bulletin & Review, 6, 57-86.
Dodson, C. S., Prinzmetal, W., & Shimamura, A. P. (1998). Using Excel to estimate parameters from observed data: An example from source memory data. Behavior Research Methods, Instruments, & Computers, 30, 517-526.
Jacoby, L., Bishara, A., Hessels, S., & Toth, J. (2005). Aging, subjective experience, and cognitive control: Dramatic false remembering by older Adults. Journal of Experimental Psychology: General, 134, 131-148.
Wickelmaier, F. (2010). mpt: Multinomial Processing Tree (MPT) Models. R package version 0.2-0. http://CRAN.R-project.org/package=mpt ワークショップのpdf
Psychoco 2011: International Workshop on Psychometric Computing