myscratchpad

# 変数の名前を変える
library(reshape)
dat <- iris
dat <- rename(dat, c(Species="newname"))
names(dat)

# 変数の再コード化
## サンプルデータ
x <- trunc(rnorm(30, m=100, sd=10))
## 90未満、90以上110未満、110以上でlow, middle, highにする
# cut関数よりも楽で以上、未満、以下が定義しやすい
nx <- x
nx[x >= 110] <- "high"
nx[x >= 90 & x < 110] <- "middle"
nx[x < 90] <- "low"
x
nx <- factor(nx) # 順序尺度にするときはordered(nx, levels=c("high", "middle", "low"))

# cut関数でコード化
x <- rnorm(1000, mean=100, sd=10)
## 中央値分割
cut(x, breaks=c(-Inf, median(x), Inf), label=c("LOW", "HIGH"), right=T, ordered_result=T)
## ordered_result=Tとすると順序尺度になる。分割は以上-未満で行われる。right=Tとしないとへんな結果になる

## 中央値分割
nx <- x
nx[x >= median(x)] <- "high"
nx[x < median(x)] <- "low"
nx <- factor(nx)
summary(nx)

# 数値変数を因子に変換する
## サンプルデータ
dat <- data.frame(matrix(sample(1:5, 100, replace=T), 10))
## 1から5列目だけ因子変数にする
dat2 <- dat
dat2[1:5] <- data.frame(lapply(dat[1:5], factor))
sapply(dat2, class) # 確認

# ダミー変数の生成
x <- gl(5, 4, labels=c("high", "mid", "low", "out", "ukn"))
lvs <- levels(x)
nlv <- nlevels(x)
x.dmy <- data.frame(matrix(0, nrow=20, ncol=5))
for (i in 1:nlv) {
x.dmy[,i] <- as.integer(x==lvs[i])
}
data.frame(x.dmy, x)

# 変数の水準のシャッフル
x <- gl(6, 5, labels=c("G", "T", "D", "C", "B", "S"))
x
x[order(rnorm(length(x)))]

# ベクトルをコンマ区切りの文字列にする
x <- gl(6, 5, labels=c("G", "T", "D", "C", "B", "S"))
paste(x, collapse=",")

SPSSのフリーのクローン、PSPPというのがあるらしい。
SPSSのsavというデータファイルも開くことができる
http://www.gnu.org/software/pspp/

RやOpenOfficeでもそうだけど、この世には限られた自分の時間を
こういうことに使う人がいて、有料のソフトやコンテンツはさらに高機能化が
求められるだろう。

by関数とAnovaのerror引数を使ってプールした誤差項による検定を行う
こちらのページを参考にした。感謝。卒論でやったという。最近の卒論はレベルが高いのう
データはテクニカルブックp95
2要因とも対応がない場合

# データ入力
dat <- data.frame(
s=factor(paste("s", 1:40, sep="")),
a=factor(c(rep("a1", 20), rep("a2", 20))),
b=factor(rep(c(rep("b1",5), rep("b2",5), rep("b3",5), rep("b4",5)),2)),
result = c(3,3,1,3,5, 4,3,4,5,7, 6,6,6,4,8, 5,7,8,7,9, 3,5,2,4,6, 2,6,3,6,4, 3,2,3,6,5, 2,3,3,4,6))

# 合計を出す
aggregate(dat[4], list(dat[,3], dat[,2]), sum)

# lmとcarのAnovaで分析
options(contrasts = c("contr.sum", "contr.sum"))
lmres <- lm(result~a*b, data=dat)
library(car)
Anovares <- Anova(lmres, type=3)

# a要因についてプールした誤差項で検定
by(dat, dat$a, function(x) Anova(lm(result~b, data=x), type=2, error=lmres))
# b要因についてプールした誤差項で検定
by(dat, dat$b, function(x) Anova(lm(result~a, data=x), type=2, error=lmres))
## セルサイズが等しいため、テクニカルブックp101どおりの結果になる
## by関数がこんなに便利とは、Anovaで誤差項が指定できるとは知らなかった

## 問題
上記の分析で、type=3とするとプールした誤差項を使わず、そのまんま水準別で検定したのと同じ結果になる というか、linear.hypothesisでerror.SSやerror.dfの情報を使っていないような…
car:::linear.hypothesis.lm
car:::linear.hypothesis.default

…分散分析ばっかりふえるのう…

Andrew Yonelinasの再認記憶2過程信号検出理論 (dual-process signal-detection model) の推定をRでやってみる
本人のホームページ: http://psychology.ucdavis.edu/labs/Yonelinas/index.html
ExcelソルバーでのサンプルをRでやっただけ

# データはExcelファイル (DPSDSSE) に従う
hit <- c(0.87, 0.75, 0.61, 0.47, 0.34) # recognition課題のヒット率
fa <- c(0.71, 0.5, 0.31, 0.14, 0.05) # フォールス・アラーム
est.v.start <- c(rep(0.5, 7)) # 推定するのは7つ。ro, d', cが5つ。初期値を全部0.5にする

# 推定式から出るヒット率、フォールス・アラーム率の予測値と、それぞれの実測値との二乗誤差を最小化するパラメータを出す
# 最適化のための関数定義
fc <- function(x){
ro <- x[1]
rn <- 0
dp <- x[2]
cv <- x[3:7]
oldvar <- 1
prehit <- ro+((1-ro)*(pnorm(cv, -dp, oldvar)))
prefa <- (1-rn)*(pnorm(cv, 0,1))
difhit <- (hit-prehit)^2
diffa <- (fa-prefa)^2
sum(difhit+diffa)
}

# optim関数で最適化
(opres <- optim(par=est.v.start, fn=fc, method="BFGS"))
est.v <- opres$par
names(est.v) <- c("ro", "dprime", "c1","c2","c3","c4","c5")
est.v # 最終的な推定値

# 改めて予測を出す
ro <- est.v[1]
rn <- 0
dp <- est.v[2]
cv <- est.v[3:7]
oldvar <- 1
prehit <- ro+((1-ro)*(pnorm(cv, -dp, oldvar)))
prefa <- (1-rn)*(pnorm(cv, 0,1))
pre.o <- cbind(prehit=prehit, hit, prefa=prefa, fa) # 予測、実測をまとめる
rownames(pre.o) <- NULL
pre.o
fc(est.v) # 最終的な二乗誤差


# プロット
plot(fa, hit, xlim=c(0,1), ylim=c(0,1), xaxs="i", yaxs="i")
# 曲線を描く。
# 適当なx軸データ (fa) をつくり、推定値に基づいてhitを予測する。その上でラインを引く。
plfa <- c(0.001,0.01,0.02,0.035,0.05,0.1,0.15,0.2,0.25,0.3,0.35,0.4,0.45,0.5,0.55,0.6,0.65,0.7,0.75,0.8,0.85,0.9,0.95,0.99999)
plcv <- qnorm(plfa, mean=0, sd=1)
plhit <- plfa+ro+(((1-ro)*(pnorm(plcv, -dp, oldvar)))-((1-rn)*(pnorm(plcv, 0, 1))))
sum(plhit>1) # 1を超えるかどうか調べる
lines(spline(plfa, plhit))
abline(0,1) # 対角線
# もっとうまい描きかたがあるだろうな…

基本の初歩の入門
L(p) = p³(1-p)^(5-3)
のL(p)を最大化するpを求める
# よくある例: 5回コインを投げて3回表、2回裏が出る確率 (L(p)) を最大化する、コインの表の出る確率 (p) は？

# Excelのソルバー
A2に初期値として0.5を入れる
B2に式を書く
=(A2^3)*((1-A2)^2)
ソルバーを起動し、
目的セル: $B$2
目標値: 最大化
変化させるセル: $A$2
制約条件: $A$2<=1, $A$2>=0
として実行をクリック

A2 (p) は0.5、B2 (L(p)) は0.03456が記入される


# 同じことをRでやるにはoptim関数を使う
fc <- function(x){
(x^3)*((1-x)^2) # 関数定義。上の=(A2^3)*((1-A2)^2)　と同じ
}
optim(par=0.5, fn=fc, control=list(fnscale=-1))
## parは初期値、fnは関数、controlは制御構文。fnscale=-1は最大化しろということ
## 返り値でconvergenceが0だと収束したということ

# 追記。Excelのソルバーは準ニュートン法か共役勾配 (傾斜) 法を用いている。オプションで指定できる。
# Rではmethodで指定する。大概制約条件があるので、"L-BFGS-B" であろう (デフォルトは"Nelder-Mead") 。Rjpwikiの記事を参考。感謝
optim(par=0.5, fn=fc, method="L-BFGS-B", upper=1, lower=0, control=list(fnscale=-1))


# optimize関数もある
optimize(f=fc, interval=c(0,1), maximum=T)
## fは関数, intervalは範囲、maximumは最大化