×
[PR]上記の広告は3ヶ月以上新規記事投稿のないブログに表示されています。新しい記事を書く事で広告が消えます。
cprp <- function(x, p=1) {
xn <- addmargins(x)
cp <- t(t(xn)/xn[nrow(xn),]*p)[-nrow(xn), -ncol(xn)]
rp <- (xn/xn[, ncol(xn)]*p)[-nrow(xn), -ncol(xn)]
pt <- prop.table(x)
return(list(adm=xn, cp=cp, rp=rp, pt=pt))
}
x <- matrix(c(10,10,20, 15, 5,20, 25,10,5, 5,5,10), nr=3,
dimnames=list(paste("r", 1:3, sep=""), paste("c", 1:4, sep="")))
x
cprp(x)
PR
right 引数
logical; if TRUE, the histogram cells are right-closed (left open) intervals.
デフォルトはright=TRUEで、より上～以下
right=FALSEにすると以上～未満
区間の下限(左側)は区間に含まれるが、上限(右側)は含まれない
デフォルトはこの逆
x <- c(1,2, 4,5)
par(mfrow=c(1,2))
hist(x)
hist(x, right=F)
# ちょっと復習
## 心理統計学の基礎、p. 9 変化量データ。
y1 <- c(4, 3, -3, 4, 1, 4, 0, 6, 2, -6, 2, 6, 6, -2, 1, 1, 6, -2, 7, -1) # 男子
y2 <- c(9, 0, 6, 5, 5, 5, 2, 11, 5, 3, 4, 7, 4, 10, -2, 6, 1, 2, 4, 4) # 女子
n1 <- 20
n2 <- 20
m1 <- mean(y1)
m2 <- mean(y2)
(tres <- t.test(y1, y2, var.equal=T))
# 標本分散の場合
s1 <- var(y1)*((n1-1)/(n1))
s2 <- var(y2)*((n2-1)/(n2))
sps <- sqrt(((s1*n1+s2*n2)/(n1+n2-2)))
(m1-m2)/sps
# 不偏分散の場合
v1 <- var(y1)
v2 <- var(y2)
spe <- sqrt(((v1*(n1-1)+v2*(n2-1))/(n1+n2-2)))
(m1-m2)/spe
# つまるところ、プールされた標準偏差の分子は2群それぞれの平方和を足したもので、それを(n1+n2-2) = 自由度 で割る
# Peasonのrは対応のないt検定だけ
## t.testの結果オブジェクト、もしくはt値とサンプルサイズを与えて効果量を算出する
tef <- function(tres=NULL, tv, n1, n2) {
if (is.null(tres)) {
tv <- tv
df <- n1+n2-2
}
else {
tv <- tres$statistic
df <- tres$parameter
}
#dv1 <- (2*tv)/sqrt(df)
dv2 <- (2*tv)/sqrt(df+2)
r1 <- sqrt(tv^2/(tv^2+df))
#print("# G*Powerと一致するのはdv2")
#c(dv1=dv1, dv2=dv2, r=r1)
res <- c(dv2, r1)
names(res) <- c("Cohen's d", "Pearson's r")
res
}
tef(tres)
tef(tv=2.20, n1=50, n2=60)
# 芝・南風原 (1990). 行動科学における統計解析法 p251 より
n1v <- c(50,30,29,25,20,20,15,15)
n2v <- c(60,35,30,30,26,24,18,15)
tvv <- c(2.20,0.91,0.76,1.15,-0.06,0.93,1.64,0.38)
dvs <- list()
for (i in 1:8) {
n1i <- n1v[i]
n2i <- n2v[i]
tvi <- tvv[i]
x <- tef(tv=tvi, n1=n1i, n2=n2i)
dvs[[i]] <- x
}
dvs
# 効果量dを他の色んな効果量に変換するExcelシート
http://www.psychsystems.net/Manuals/StatsCalculators/Effect_Size_Calculator%2017.xls
## 心理統計学の基礎、p. 9 変化量データ。
y1 <- c(4, 3, -3, 4, 1, 4, 0, 6, 2, -6, 2, 6, 6, -2, 1, 1, 6, -2, 7, -1) # 男子
y2 <- c(9, 0, 6, 5, 5, 5, 2, 11, 5, 3, 4, 7, 4, 10, -2, 6, 1, 2, 4, 4) # 女子
n1 <- 20
n2 <- 20
m1 <- mean(y1)
m2 <- mean(y2)
(tres <- t.test(y1, y2, var.equal=T))
# 標本分散の場合
s1 <- var(y1)*((n1-1)/(n1))
s2 <- var(y2)*((n2-1)/(n2))
sps <- sqrt(((s1*n1+s2*n2)/(n1+n2-2)))
(m1-m2)/sps
# 不偏分散の場合
v1 <- var(y1)
v2 <- var(y2)
spe <- sqrt(((v1*(n1-1)+v2*(n2-1))/(n1+n2-2)))
(m1-m2)/spe
# つまるところ、プールされた標準偏差の分子は2群それぞれの平方和を足したもので、それを(n1+n2-2) = 自由度 で割る
# Peasonのrは対応のないt検定だけ
## t.testの結果オブジェクト、もしくはt値とサンプルサイズを与えて効果量を算出する
tef <- function(tres=NULL, tv, n1, n2) {
if (is.null(tres)) {
tv <- tv
df <- n1+n2-2
}
else {
tv <- tres$statistic
df <- tres$parameter
}
#dv1 <- (2*tv)/sqrt(df)
dv2 <- (2*tv)/sqrt(df+2)
r1 <- sqrt(tv^2/(tv^2+df))
#print("# G*Powerと一致するのはdv2")
#c(dv1=dv1, dv2=dv2, r=r1)
res <- c(dv2, r1)
names(res) <- c("Cohen's d", "Pearson's r")
res
}
tef(tres)
tef(tv=2.20, n1=50, n2=60)
# 芝・南風原 (1990). 行動科学における統計解析法 p251 より
n1v <- c(50,30,29,25,20,20,15,15)
n2v <- c(60,35,30,30,26,24,18,15)
tvv <- c(2.20,0.91,0.76,1.15,-0.06,0.93,1.64,0.38)
dvs <- list()
for (i in 1:8) {
n1i <- n1v[i]
n2i <- n2v[i]
tvi <- tvv[i]
x <- tef(tv=tvi, n1=n1i, n2=n2i)
dvs[[i]] <- x
}
dvs
# 効果量dを他の色んな効果量に変換するExcelシート
http://www.psychsystems.net/Manuals/StatsCalculators/Effect_Size_Calculator%2017.xls
# 心理統計学の基礎 共分散構造分析 p351-363 の例を実行
# ...してみようと思ったら以下のURLにそのものズバリが載っていた。感謝。自分用にコピペしておく
# http://home.hiroshima-u.ac.jp/ksatoh/documents/Rsem.txt
# http://home.hiroshima-u.ac.jp/ksatoh/download.htm
library(sem)
# 表10-8 観測変数間の相関係数 p354
mat <- matrix(c(
1.000,0.160,0.302,0.461,0.299,0.152,0.134,0.182,0.251,0.372,0.157,0.203,
0.160,1.000,0.341,0.400,0.404,0.320,0.403,0.374,0.285,0.100,0.291,-0.014,
0.302,0.341,1.000,0.372,0.552,0.476,0.467,0.572,0.316,0.408,0.393,0.369,
0.461,0.400,0.372,1.000,0.302,0.225,0.256,0.255,0.164,0.236,0.229,0.224,
0.299,0.404,0.552,0.302,1.000,0.708,0.623,0.776,0.361,0.294,0.472,0.342,
0.152,0.320,0.476,0.225,0.708,1.000,0.324,0.769,0.295,0.206,0.351,0.202,
0.134,0.403,0.467,0.256,0.623,0.324,1.000,0.724,0.260,0.071,0.204,0.152,
0.182,0.374,0.572,0.255,0.776,0.769,0.724,1.000,0.284,0.142,0.320,0.189,
0.251,0.285,0.316,0.164,0.361,0.295,0.260,0.284,1.000,0.295,0.290,0.418,
0.372,0.100,0.408,0.236,0.294,0.206,0.071,0.142,0.295,1.000,0.468,0.351,
0.157,0.291,0.393,0.229,0.472,0.351,0.204,0.320,0.290,0.468,1.000,0.385,
0.203,-0.014,0.369,0.224,0.342,0.202,0.152,0.189,0.418,0.351,0.385,1.000),
nr=12,
dimnames = list(paste("y", 1:12, sep=""), paste("y", 1:12, sep=""))
)
mat
# 日本語じゃなくて英語にしておく
model.agg <- specify.model()
mot -> y1, b11, NA
mot -> y2, b12, NA
mot -> y3, b13, NA
mot -> y4, b14, NA
int -> y5, NA, 1
int -> y6, b22, NA
int -> y7, b23, NA
int -> y8, b24, NA
agg -> y9, NA, 1
agg -> y10, b32, NA
agg -> y11, b33, NA
agg -> y12, b34, NA
mot -> int, g11, NA
mot -> agg, g12, NA
y1 <-> y1, e1, NA
y2 <-> y2, e2, NA
y3 <-> y3, e3, NA
y4 <-> y4, e4, NA
y5 <-> y5, e5, NA
y6 <-> y6, e6, NA
y7 <-> y7, e7, NA
y8 <-> y8, e8, NA
y9 <-> y9, e9, NA
y10 <-> y10, e10, NA
y11 <-> y11, e11, NA
y12 <-> y12, e12, NA
mot <-> mot, NA, 1
int <-> int, delta2, NA
agg <-> agg, delta3, NA
sem.agg <- sem(model.agg, mat, N=50)
summary(sem.agg)
std.coef(sem.agg,digit=4) # 標準解 図10-6と一致するのはこっち
# 図10-6 くらいは自分で書いてみよう
library(Rgraphviz)
path.diagram(sem.agg, ignore.double=FALSE, edge.labels="values", digits=3, file="out2")
# やはり、日本語フォントでは表示できない。Graphvizのバージョンが古いのかもしれない
# しかたないので協調性はagg, 母親価値はmot, 相互作用経験はintとした。でも文字の一部が消えたり、どうもうまくいかない。
# そのうえ、sem.aggの係数で描画するため、標準解を出力できない
# ...してみようと思ったら以下のURLにそのものズバリが載っていた。感謝。自分用にコピペしておく
# http://home.hiroshima-u.ac.jp/ksatoh/documents/Rsem.txt
# http://home.hiroshima-u.ac.jp/ksatoh/download.htm
library(sem)
# 表10-8 観測変数間の相関係数 p354
mat <- matrix(c(
1.000,0.160,0.302,0.461,0.299,0.152,0.134,0.182,0.251,0.372,0.157,0.203,
0.160,1.000,0.341,0.400,0.404,0.320,0.403,0.374,0.285,0.100,0.291,-0.014,
0.302,0.341,1.000,0.372,0.552,0.476,0.467,0.572,0.316,0.408,0.393,0.369,
0.461,0.400,0.372,1.000,0.302,0.225,0.256,0.255,0.164,0.236,0.229,0.224,
0.299,0.404,0.552,0.302,1.000,0.708,0.623,0.776,0.361,0.294,0.472,0.342,
0.152,0.320,0.476,0.225,0.708,1.000,0.324,0.769,0.295,0.206,0.351,0.202,
0.134,0.403,0.467,0.256,0.623,0.324,1.000,0.724,0.260,0.071,0.204,0.152,
0.182,0.374,0.572,0.255,0.776,0.769,0.724,1.000,0.284,0.142,0.320,0.189,
0.251,0.285,0.316,0.164,0.361,0.295,0.260,0.284,1.000,0.295,0.290,0.418,
0.372,0.100,0.408,0.236,0.294,0.206,0.071,0.142,0.295,1.000,0.468,0.351,
0.157,0.291,0.393,0.229,0.472,0.351,0.204,0.320,0.290,0.468,1.000,0.385,
0.203,-0.014,0.369,0.224,0.342,0.202,0.152,0.189,0.418,0.351,0.385,1.000),
nr=12,
dimnames = list(paste("y", 1:12, sep=""), paste("y", 1:12, sep=""))
)
mat
# 日本語じゃなくて英語にしておく
model.agg <- specify.model()
mot -> y1, b11, NA
mot -> y2, b12, NA
mot -> y3, b13, NA
mot -> y4, b14, NA
int -> y5, NA, 1
int -> y6, b22, NA
int -> y7, b23, NA
int -> y8, b24, NA
agg -> y9, NA, 1
agg -> y10, b32, NA
agg -> y11, b33, NA
agg -> y12, b34, NA
mot -> int, g11, NA
mot -> agg, g12, NA
y1 <-> y1, e1, NA
y2 <-> y2, e2, NA
y3 <-> y3, e3, NA
y4 <-> y4, e4, NA
y5 <-> y5, e5, NA
y6 <-> y6, e6, NA
y7 <-> y7, e7, NA
y8 <-> y8, e8, NA
y9 <-> y9, e9, NA
y10 <-> y10, e10, NA
y11 <-> y11, e11, NA
y12 <-> y12, e12, NA
mot <-> mot, NA, 1
int <-> int, delta2, NA
agg <-> agg, delta3, NA
sem.agg <- sem(model.agg, mat, N=50)
summary(sem.agg)
std.coef(sem.agg,digit=4) # 標準解 図10-6と一致するのはこっち
# 図10-6 くらいは自分で書いてみよう
library(Rgraphviz)
path.diagram(sem.agg, ignore.double=FALSE, edge.labels="values", digits=3, file="out2")
# やはり、日本語フォントでは表示できない。Graphvizのバージョンが古いのかもしれない
# しかたないので協調性はagg, 母親価値はmot, 相互作用経験はintとした。でも文字の一部が消えたり、どうもうまくいかない。
# そのうえ、sem.aggの係数で描画するため、標準解を出力できない
# 心理統計学の基礎, p347-350のもの。感謝
# 表10-1 p318 相関行列
mat <- matrix(c(
1.000,0.033,0.315,0.456,0.266,0.607,0.228,0.419,
0.033,1.000,0.637,0.250,0.528,0.195,0.522,0.420,
0.315,0.637,1.000,0.333,0.880,0.237,0.750,0.328,
0.456,0.250,0.333,1.000,0.362,0.432,0.398,0.449,
0.266,0.528,0.880,0.362,1.000,0.252,0.738,0.269,
0.607,0.195,0.237,0.432,0.252,1.000,0.335,0.463,
0.228,0.522,0.750,0.398,0.738,0.335,1.000,0.238,
0.419,0.420,0.328,0.449,0.269,0.463,0.238,1.000),
nrow=8, dimnames=list(c("onw", "yok", "gai", "sin", "sha", "kyo", "sek", "sun"), c("onw", "yok", "gai", "sin", "sha", "kyo", "sek", "sun"))
)
library(sem)
model <- specify.model()
f1 -> yok, b12, 1
f1 -> gai, b13, NA
f1 -> sha, b15, NA
f1 -> sek, b17, NA
f2 -> onw, b21, NA
f2 -> sin, b24, NA
f2 -> kyo, b26, NA
f2 -> sun, b28, NA
f1 <-> f2, g11, NA
f1 <-> f1, NA, 1
f2 <-> f2, NA,1
onw <-> onw, e1, NA
yok <-> yok, e2, NA
gai <-> gai, e3, NA
sin <-> sin, e4, NA
sha <-> sha, e5, NA
kyo <-> kyo, e6, NA
sek <-> sek, e7, NA
sun <-> sun, e8, NA
semres <- sem(model, mat, N=200)
# 表10-7, p394 のものと一致
summary(semres)
# 適合度指標基準の覚書
# 以下参考にしたサイト。感謝。David Kennyのをちゃんと読むのが一番いいだろう
# David Kenny
# 小塩先生のページ。感謝
# Rでsem入門。
# Indiana Universityとpdf
# 荘島先生のサイト
# GFI, AGFIは.90以上。サンプルサイズが大きいときは無視
# RMSEAは小さいほうがいい。0.05未満
# AIC, BICは相対基準。小さいほうがいい
# 表10-1 p318 相関行列
mat <- matrix(c(
1.000,0.033,0.315,0.456,0.266,0.607,0.228,0.419,
0.033,1.000,0.637,0.250,0.528,0.195,0.522,0.420,
0.315,0.637,1.000,0.333,0.880,0.237,0.750,0.328,
0.456,0.250,0.333,1.000,0.362,0.432,0.398,0.449,
0.266,0.528,0.880,0.362,1.000,0.252,0.738,0.269,
0.607,0.195,0.237,0.432,0.252,1.000,0.335,0.463,
0.228,0.522,0.750,0.398,0.738,0.335,1.000,0.238,
0.419,0.420,0.328,0.449,0.269,0.463,0.238,1.000),
nrow=8, dimnames=list(c("onw", "yok", "gai", "sin", "sha", "kyo", "sek", "sun"), c("onw", "yok", "gai", "sin", "sha", "kyo", "sek", "sun"))
)
library(sem)
model <- specify.model()
f1 -> yok, b12, 1
f1 -> gai, b13, NA
f1 -> sha, b15, NA
f1 -> sek, b17, NA
f2 -> onw, b21, NA
f2 -> sin, b24, NA
f2 -> kyo, b26, NA
f2 -> sun, b28, NA
f1 <-> f2, g11, NA
f1 <-> f1, NA, 1
f2 <-> f2, NA,1
onw <-> onw, e1, NA
yok <-> yok, e2, NA
gai <-> gai, e3, NA
sin <-> sin, e4, NA
sha <-> sha, e5, NA
kyo <-> kyo, e6, NA
sek <-> sek, e7, NA
sun <-> sun, e8, NA
semres <- sem(model, mat, N=200)
# 表10-7, p394 のものと一致
summary(semres)
# 適合度指標基準の覚書
# 以下参考にしたサイト。感謝。David Kennyのをちゃんと読むのが一番いいだろう
# David Kenny
# 小塩先生のページ。感謝
# Rでsem入門。
# Indiana Universityとpdf
# 荘島先生のサイト
# GFI, AGFIは.90以上。サンプルサイズが大きいときは無視
# RMSEAは小さいほうがいい。0.05未満
# AIC, BICは相対基準。小さいほうがいい
