×
[PR]上記の広告は3ヶ月以上新規記事投稿のないブログに表示されています。新しい記事を書く事で広告が消えます。
psychパッケージのfa関数が不思議な挙動を見せるようなので
自分で調べてRjpwikiにレスしておいた。
せっかく調べたのでもう少し詳しく書いてみる
library(psych)
data(bfi)
x <- factor.pa(bfi, nfactors=3, rotate="promax")
# psych:::print.psych.faの中身を途中まで実行
# パラメータを設定しておく
digits = 2; all = FALSE; cut = NULL; sort = TRUE
#> psych:::print.psych.fa
#function (x, digits = 2, all = FALSE, cut = NULL, sort = FALSE,
# ...)
#{ ここまでは定義なので実行しないこと
if (!is.null(x$fn)) {
if (x$fn == "principal") {
cat("Principal Components Analysis")
}
else {
cat("Factor Analysis using method = ", x$fm)
}
}
cat("\nCall: ")
print(x$Call)
load <- x$loadings
if (is.null(cut))
cut <- 0
nitems <- dim(load)[1]
nfactors <- dim(load)[2]
loads <- data.frame(item = seq(1:nitems), cluster = rep(0,
nitems), unclass(load))
if (sort) {
loads$cluster <- apply(abs(load), 1, which.max)
ord <- sort(loads$cluster, index.return = TRUE)
loads[1:nitems, ] <- loads[ord$ix, ]
rownames(loads)[1:nitems] <- rownames(loads)[ord$ix]
items <- table(loads$cluster)
first <- 1
item <- loads$item
for (i in 1:length(items)) {
if (items[i] > 0) {
last <- first + items[i] - 1
ord <- sort(abs(loads[first:last, i + 2]), decreasing = TRUE,
index.return = TRUE)
loads[first:last, 3:(nfactors + 2)] <- load[item[ord$ix +
first - 1], ]
loads[first:last, 1] <- item[ord$ix + first -
1]
rownames(loads)[first:last] <- rownames(loads)[ord$ix +
first - 1]
first <- first + items[i]
}
}
}
if (cut > 0) {
ncol <- dim(loads)[2] - 2
rloads <- round(loads, digits)
fx <- format(rloads, digits = digits)
nc <- nchar(fx[1, 3], type = "c")
fx.1 <- fx[, 1, drop = FALSE]
fx.2 <- fx[, 3:(2 + ncol)]
load.2 <- as.matrix(loads[, 3:(ncol + 2)])
fx.2[abs(load.2) < cut] <- paste(rep(" ", nc), collapse = "")
fx <- data.frame(V = fx.1, fx.2)
if (dim(fx)[2] < 3)
colnames(fx) <- c("V", colnames(x$loadings))
if (nfactors > 1) {
if (is.null(x$Phi)) {
h2 <- rowSums(load.2^2)
}
else {
h2 <- diag(load.2 %*% x$Phi %*% t(load.2))
}
}
else {
h2 <- load.2^2
}
if (!is.null(x$uniquenesses)) {
u2 <- x$uniquenesses
}
else {
u2 <- (1 - h2)
}
vtotal <- sum(h2 + u2)
if (isTRUE(all.equal(vtotal, nitems))) {
cat("Standardized loadings based upon correlation matrix\n")
print(cbind(fx, h2, u2), quote = "FALSE", digits = digits)
}
else {
cat("Unstandardized loadings based upon covariance matrix\n")
print(cbind(fx, h2, u2, H2 = h2/(h2 + u2), U2 = u2/(h2 +
u2)), quote = "FALSE", digits = digits)
}
} else {
ncol <- dim(loads)[2] - 2
load.2 <- as.matrix(loads[, 3:(ncol + 2)])
if (nfactors > 1) {
if (is.null(x$Phi)) {
h2 <- rowSums(load.2^2)
}
else {
h2 <- diag(load.2 %*% x$Phi %*% t(load.2))
}
}
else {
h2 <- load.2^2
}
if (!is.null(x$uniquenesses)) {
u2 <- x$uniquenesses
}
else {
u2 <- (1 - h2)
}
vtotal <- sum(h2 + u2)
if (isTRUE(all.equal(vtotal, nitems))) {
cat("Standardized loadings based upon correlation matrix\n")
print(cbind(round(loads[-(1:2)], digits), h2, u2),
quote = "FALSE", digits = digits)
}
else {
cat("Unstandardized loadings based upon covariance matrix\n")
print(cbind(round(loads[-c(1:2)], digits), h2, u2,
H2 = h2/(h2 + u2), U2 = u2/(h2 + u2)), quote = "FALSE",
digits = digits)
}
}
if (is.null(x$Phi)) {
if (nfactors > 1) {
vx <- colSums(load.2^2)
}
else {
vx <- sum(load.2^2)
}
} else {
vx <- diag(x$Phi %*% t(load) %*% load)
}
names(vx) <- colnames(x$loadings)
varex <- rbind(`SS loadings` = vx)
varex <- rbind(varex, `Proportion Var` = vx/vtotal)
if (nfactors > 1)
varex <- rbind(varex, `Cumulative Var` = cumsum(vx/vtotal))
cat("\n")
print(round(varex, digits))
## この時点でh2はソートされているが、u2はソートされていない
h2
u2
# この後のコードでこの両者を使って標準化している (fx.2 ... sqrt(h2 + u2)) ため、誤りが生じます
#if (!isTRUE(all.equal(vtotal, nitems))) {
cat("\n Standardized loadings\n")
fx <- format(loads, digits = digits)
nc <- nchar(fx[1, 3], type = "c")
fx.1 <- fx[, 1, drop = FALSE]
fx.2 <- round(loads[, 3:(2 + ncol)]/sqrt(h2 + u2), digits)
########
# 誤った負荷量
fx.2
# 修正コード
u2 <- 1 - h2
fx.2 <- round(loads[, 3:(2 + ncol)]/sqrt(h2 + u2), digits)
fx.2
########
load.2 <- loads[, 3:(ncol + 2)]/sqrt(h2 + u2)
fx.2[abs(load.2) < cut] <- paste(rep(" ", nc), collapse = "")
fx <- data.frame(V = fx.1, fx.2)
## 以下本題
一エンドユーザーとしてみると、RがSPSSやSASより優れている点はフリーであること以外にない。Rを利用しているえらい人は (自分の本を売るためか) Rの使用をすすめる人が多いが、個人的には止めたほうがいいと思う。
それは使用に慣れるまでが大変で、容易に間違った結果を出力するからである。巷にあふれるRの入門書も誤った使用法を紹介しているものが山ほどある (厳密にいえば誤っているわけではなく、教科書の例題のようなごく限られた場合にのみ適用可能なコードを紹介している) …ような気がする。
特に、「世界中の人が使ってチェックするので，計算の信頼性もかなり高い。」といった宣伝文句はどうかと思う。baseに含まれている関数には当てはまるが、あまたあるパッケージの関数には当てはまらない。というか、下手に使うと間違った結果しか出力されず、信頼性はかなり低い。なぜなら、上記のようにフリーであるので誰にでも使用できるが、そのため誤ったコード例が量産されている (本サイトもかなり誤りがあるんです) 。エンドユーザーはコピペのように手軽に実行でき、しかも正しい結果の得られるものを求める。Rで自らのデータに合わせて適切にオプションを選択し、分析を実行していくのは容易ではない。
コピペで間違ったら自己責任でしょ、フリーなんだし、自分の研究なんだから…というのであれば、やはりSPSSやSASの方が優れているとしか思えない。
「Rで分析しました」という記述があったら、その結果を疑わしく見てしまうのは自分だけかな。「下手に使って、どっか間違ってんじゃないの？」と思ってしまう。
自分で調べてRjpwikiにレスしておいた。
せっかく調べたのでもう少し詳しく書いてみる
library(psych)
data(bfi)
x <- factor.pa(bfi, nfactors=3, rotate="promax")
# psych:::print.psych.faの中身を途中まで実行
# パラメータを設定しておく
digits = 2; all = FALSE; cut = NULL; sort = TRUE
#> psych:::print.psych.fa
#function (x, digits = 2, all = FALSE, cut = NULL, sort = FALSE,
# ...)
#{ ここまでは定義なので実行しないこと
if (!is.null(x$fn)) {
if (x$fn == "principal") {
cat("Principal Components Analysis")
}
else {
cat("Factor Analysis using method = ", x$fm)
}
}
cat("\nCall: ")
print(x$Call)
load <- x$loadings
if (is.null(cut))
cut <- 0
nitems <- dim(load)[1]
nfactors <- dim(load)[2]
loads <- data.frame(item = seq(1:nitems), cluster = rep(0,
nitems), unclass(load))
if (sort) {
loads$cluster <- apply(abs(load), 1, which.max)
ord <- sort(loads$cluster, index.return = TRUE)
loads[1:nitems, ] <- loads[ord$ix, ]
rownames(loads)[1:nitems] <- rownames(loads)[ord$ix]
items <- table(loads$cluster)
first <- 1
item <- loads$item
for (i in 1:length(items)) {
if (items[i] > 0) {
last <- first + items[i] - 1
ord <- sort(abs(loads[first:last, i + 2]), decreasing = TRUE,
index.return = TRUE)
loads[first:last, 3:(nfactors + 2)] <- load[item[ord$ix +
first - 1], ]
loads[first:last, 1] <- item[ord$ix + first -
1]
rownames(loads)[first:last] <- rownames(loads)[ord$ix +
first - 1]
first <- first + items[i]
}
}
}
if (cut > 0) {
ncol <- dim(loads)[2] - 2
rloads <- round(loads, digits)
fx <- format(rloads, digits = digits)
nc <- nchar(fx[1, 3], type = "c")
fx.1 <- fx[, 1, drop = FALSE]
fx.2 <- fx[, 3:(2 + ncol)]
load.2 <- as.matrix(loads[, 3:(ncol + 2)])
fx.2[abs(load.2) < cut] <- paste(rep(" ", nc), collapse = "")
fx <- data.frame(V = fx.1, fx.2)
if (dim(fx)[2] < 3)
colnames(fx) <- c("V", colnames(x$loadings))
if (nfactors > 1) {
if (is.null(x$Phi)) {
h2 <- rowSums(load.2^2)
}
else {
h2 <- diag(load.2 %*% x$Phi %*% t(load.2))
}
}
else {
h2 <- load.2^2
}
if (!is.null(x$uniquenesses)) {
u2 <- x$uniquenesses
}
else {
u2 <- (1 - h2)
}
vtotal <- sum(h2 + u2)
if (isTRUE(all.equal(vtotal, nitems))) {
cat("Standardized loadings based upon correlation matrix\n")
print(cbind(fx, h2, u2), quote = "FALSE", digits = digits)
}
else {
cat("Unstandardized loadings based upon covariance matrix\n")
print(cbind(fx, h2, u2, H2 = h2/(h2 + u2), U2 = u2/(h2 +
u2)), quote = "FALSE", digits = digits)
}
} else {
ncol <- dim(loads)[2] - 2
load.2 <- as.matrix(loads[, 3:(ncol + 2)])
if (nfactors > 1) {
if (is.null(x$Phi)) {
h2 <- rowSums(load.2^2)
}
else {
h2 <- diag(load.2 %*% x$Phi %*% t(load.2))
}
}
else {
h2 <- load.2^2
}
if (!is.null(x$uniquenesses)) {
u2 <- x$uniquenesses
}
else {
u2 <- (1 - h2)
}
vtotal <- sum(h2 + u2)
if (isTRUE(all.equal(vtotal, nitems))) {
cat("Standardized loadings based upon correlation matrix\n")
print(cbind(round(loads[-(1:2)], digits), h2, u2),
quote = "FALSE", digits = digits)
}
else {
cat("Unstandardized loadings based upon covariance matrix\n")
print(cbind(round(loads[-c(1:2)], digits), h2, u2,
H2 = h2/(h2 + u2), U2 = u2/(h2 + u2)), quote = "FALSE",
digits = digits)
}
}
if (is.null(x$Phi)) {
if (nfactors > 1) {
vx <- colSums(load.2^2)
}
else {
vx <- sum(load.2^2)
}
} else {
vx <- diag(x$Phi %*% t(load) %*% load)
}
names(vx) <- colnames(x$loadings)
varex <- rbind(`SS loadings` = vx)
varex <- rbind(varex, `Proportion Var` = vx/vtotal)
if (nfactors > 1)
varex <- rbind(varex, `Cumulative Var` = cumsum(vx/vtotal))
cat("\n")
print(round(varex, digits))
## この時点でh2はソートされているが、u2はソートされていない
h2
u2
# この後のコードでこの両者を使って標準化している (fx.2 ... sqrt(h2 + u2)) ため、誤りが生じます
#if (!isTRUE(all.equal(vtotal, nitems))) {
cat("\n Standardized loadings\n")
fx <- format(loads, digits = digits)
nc <- nchar(fx[1, 3], type = "c")
fx.1 <- fx[, 1, drop = FALSE]
fx.2 <- round(loads[, 3:(2 + ncol)]/sqrt(h2 + u2), digits)
########
# 誤った負荷量
fx.2
# 修正コード
u2 <- 1 - h2
fx.2 <- round(loads[, 3:(2 + ncol)]/sqrt(h2 + u2), digits)
fx.2
########
load.2 <- loads[, 3:(ncol + 2)]/sqrt(h2 + u2)
fx.2[abs(load.2) < cut] <- paste(rep(" ", nc), collapse = "")
fx <- data.frame(V = fx.1, fx.2)
## 以下本題
一エンドユーザーとしてみると、RがSPSSやSASより優れている点はフリーであること以外にない。Rを利用しているえらい人は (自分の本を売るためか) Rの使用をすすめる人が多いが、個人的には止めたほうがいいと思う。
それは使用に慣れるまでが大変で、容易に間違った結果を出力するからである。巷にあふれるRの入門書も誤った使用法を紹介しているものが山ほどある (厳密にいえば誤っているわけではなく、教科書の例題のようなごく限られた場合にのみ適用可能なコードを紹介している) …ような気がする。
特に、「世界中の人が使ってチェックするので，計算の信頼性もかなり高い。」といった宣伝文句はどうかと思う。baseに含まれている関数には当てはまるが、あまたあるパッケージの関数には当てはまらない。というか、下手に使うと間違った結果しか出力されず、信頼性はかなり低い。なぜなら、上記のようにフリーであるので誰にでも使用できるが、そのため誤ったコード例が量産されている (本サイトもかなり誤りがあるんです) 。エンドユーザーはコピペのように手軽に実行でき、しかも正しい結果の得られるものを求める。Rで自らのデータに合わせて適切にオプションを選択し、分析を実行していくのは容易ではない。
コピペで間違ったら自己責任でしょ、フリーなんだし、自分の研究なんだから…というのであれば、やはりSPSSやSASの方が優れているとしか思えない。
「Rで分析しました」という記述があったら、その結果を疑わしく見てしまうのは自分だけかな。「下手に使って、どっか間違ってんじゃないの？」と思ってしまう。
PR
# 交互作用重回帰の別の例。以下参考。感謝
# http://wilhelmhofmann.info/wordpress/?page_id=120
library(foreign)
dat <- read.spss("https://blog.cnobi.jp/v1/blog/user/89d80905c7038b4121822249e9062fba/1286213208", to.data.frame=T, max.value.labels=99)
head(dat)
# 記述統計
library(psych)
describe(dat)
x <- dat[c(6,2,3,4)]
attach(dat)
heightc <- scale(height, scale=F)
agec <- scale(age, scale=F)
summary(lm(lifesat~heightc*agec))
# simple slope
agech <- agec-sd(agec)
hmod <- lm(lifesat~heightc*agech)
summary(hmod)
agecl <- agec+sd(agec)
lmod <- lm(lifesat~heightc*agecl)
summary(lmod)
# プロット
x <- agec
y <- lifesat
plot(x,y, col=0)
abline(hmod, lty=1)
abline(lmod, lty=2)
legend("topright", legend=c("ageabove", "agebelow"), lty=1:2)
# 正しい標準化値。交互作用項を作成してから標準化するのは誤り
hsc <- scale(height)
asc <- scale(age)
lsc <- scale(lifesat)
summary(lm(lsc~hsc*asc))
# simple slopeの標準化値
asch <- asc-1
hmods <- lm(lsc~hsc*asch)
summary(hmods) # 高群
agecl <- asc+1
lmods <- lm(lsc~hsc*agecl)
summary(lmods) # 低群
# 標準化値のプロット
win.graph()
x <- asc
y <- lsc
plot(x,y, col=0)
abline(hmods, lty=1)
abline(lmods, lty=2)
legend("topright", legend=c("ageabove", "agebelow"), lty=1:2)
# そのうちやろう
# 統制変数
# 連続 x カテゴリカル
summary(lm(lifesat~heightc*gender))
# simple slope
mm <- factor(gender, levels=c("male", "female")) # 男性をリファレンスカテゴリにする。ダミーコードではmale=0, female=1になる。
ff <- factor(gender, levels=c("female", "male"))
summary(lm(lifesat~heightc*mm)) # 男性
summary(lm(lifesat~heightc*ff)) # 女性
# 標準化値
# The unstandardized solution from the output is the correct solution (Friedrich, 1982)!
summary(lm(lsc~hsc*gender))
# .507 = estimated difference in regression weights between groups
# simple slope. カテゴリカル変数は標準化しないこと
summary(lm(lsc~hsc*mm)) # 男性
summary(lm(lsc~hsc*ff)) # 女性
# プロット
plot(heightc, lifesat)
abline(lm(lifesat~heightc*mm), lty=1)
abline(lm(lifesat~heightc*ff), lty=2)
legend("bottomright", legend=c("male", "female"), lty=1:2)
# x軸とライン分けを変えたプロット
lmres <- lm(lifesat~heightc*gender) # 全体のモデル
coef(lmres)
inc <- coef(lmres)[1] # 全体の切片
cfidv <- coef(lmres)[2] # 独立変数の回帰係数
cfmod <- coef(lmres)[3] # 調整変数の回帰係数
cfint <- coef(lmres)[4] # 独立変数の回帰係数
sdidv <- sd(heightc)
# 全体の重回帰分析でfemaleをリファレンスカテゴリにしている場合
fh <- inc+sdidv*cfidv # 切片＋（平均＋1SD)*回帰係数
fl <- inc-sdidv*cfidv # 切片＋（平均-1SD)*回帰係数
mh <- inc+sdidv*cfidv + cfmod+((sdidv)*cfint) # 切片＋（平均＋1SD)*回帰係数 + 調整変数の回帰係数＋((平均＋1SD)*交互作用項の回帰係数)
ml <- inc-sdidv*cfidv + cfmod-((sdidv)*cfint) # 切片＋（平均-1SD)*回帰係数 + 調整変数の回帰係数＋((平均-1SD)*交互作用項の回帰係数)
mat <- matrix(c(fl, ml, fh, mh), nr=2, dimnames=list(c("female", "male"), c("heightlow", "heighthigh")))
mat
matplot(mat, type="l", xaxt="n")
axis(1,
at = 1:2,
labels=rownames(mat),
)
legend("bottomleft", legend=colnames(mat), lty=1:2, col=1:2)
# 棒グラフにしてみる
win.graph(); barplot(t(mat), beside=T, legend=c("heightlow", "heighthigh"), ylab="lifesat")
# http://wilhelmhofmann.info/wordpress/?page_id=120
library(foreign)
dat <- read.spss("https://blog.cnobi.jp/v1/blog/user/89d80905c7038b4121822249e9062fba/1286213208", to.data.frame=T, max.value.labels=99)
head(dat)
# 記述統計
library(psych)
describe(dat)
x <- dat[c(6,2,3,4)]
attach(dat)
heightc <- scale(height, scale=F)
agec <- scale(age, scale=F)
summary(lm(lifesat~heightc*agec))
# simple slope
agech <- agec-sd(agec)
hmod <- lm(lifesat~heightc*agech)
summary(hmod)
agecl <- agec+sd(agec)
lmod <- lm(lifesat~heightc*agecl)
summary(lmod)
# プロット
x <- agec
y <- lifesat
plot(x,y, col=0)
abline(hmod, lty=1)
abline(lmod, lty=2)
legend("topright", legend=c("ageabove", "agebelow"), lty=1:2)
# 正しい標準化値。交互作用項を作成してから標準化するのは誤り
hsc <- scale(height)
asc <- scale(age)
lsc <- scale(lifesat)
summary(lm(lsc~hsc*asc))
# simple slopeの標準化値
asch <- asc-1
hmods <- lm(lsc~hsc*asch)
summary(hmods) # 高群
agecl <- asc+1
lmods <- lm(lsc~hsc*agecl)
summary(lmods) # 低群
# 標準化値のプロット
win.graph()
x <- asc
y <- lsc
plot(x,y, col=0)
abline(hmods, lty=1)
abline(lmods, lty=2)
legend("topright", legend=c("ageabove", "agebelow"), lty=1:2)
# そのうちやろう
# 統制変数
# 連続 x カテゴリカル
summary(lm(lifesat~heightc*gender))
# simple slope
mm <- factor(gender, levels=c("male", "female")) # 男性をリファレンスカテゴリにする。ダミーコードではmale=0, female=1になる。
ff <- factor(gender, levels=c("female", "male"))
summary(lm(lifesat~heightc*mm)) # 男性
summary(lm(lifesat~heightc*ff)) # 女性
# 標準化値
# The unstandardized solution from the output is the correct solution (Friedrich, 1982)!
summary(lm(lsc~hsc*gender))
# .507 = estimated difference in regression weights between groups
# simple slope. カテゴリカル変数は標準化しないこと
summary(lm(lsc~hsc*mm)) # 男性
summary(lm(lsc~hsc*ff)) # 女性
# プロット
plot(heightc, lifesat)
abline(lm(lifesat~heightc*mm), lty=1)
abline(lm(lifesat~heightc*ff), lty=2)
legend("bottomright", legend=c("male", "female"), lty=1:2)
# x軸とライン分けを変えたプロット
lmres <- lm(lifesat~heightc*gender) # 全体のモデル
coef(lmres)
inc <- coef(lmres)[1] # 全体の切片
cfidv <- coef(lmres)[2] # 独立変数の回帰係数
cfmod <- coef(lmres)[3] # 調整変数の回帰係数
cfint <- coef(lmres)[4] # 独立変数の回帰係数
sdidv <- sd(heightc)
# 全体の重回帰分析でfemaleをリファレンスカテゴリにしている場合
fh <- inc+sdidv*cfidv # 切片＋（平均＋1SD)*回帰係数
fl <- inc-sdidv*cfidv # 切片＋（平均-1SD)*回帰係数
mh <- inc+sdidv*cfidv + cfmod+((sdidv)*cfint) # 切片＋（平均＋1SD)*回帰係数 + 調整変数の回帰係数＋((平均＋1SD)*交互作用項の回帰係数)
ml <- inc-sdidv*cfidv + cfmod-((sdidv)*cfint) # 切片＋（平均-1SD)*回帰係数 + 調整変数の回帰係数＋((平均-1SD)*交互作用項の回帰係数)
mat <- matrix(c(fl, ml, fh, mh), nr=2, dimnames=list(c("female", "male"), c("heightlow", "heighthigh")))
mat
matplot(mat, type="l", xaxt="n")
axis(1,
at = 1:2,
labels=rownames(mat),
)
legend("bottomleft", legend=colnames(mat), lty=1:2, col=1:2)
# 棒グラフにしてみる
win.graph(); barplot(t(mat), beside=T, legend=c("heightlow", "heighthigh"), ylab="lifesat")
# 参考URL。感謝
# http://psy.isc.chubu.ac.jp/~oshiolab/teaching_folder/datakaiseki_folder/07_folder/da07_02.html
cr <- c(1,1,1,0,1,0,0,1,0,0)
ex <- c(4,5,3,2,4,3,2,4,3,1)
ky <- c(3,4,4,5,4,2,3,5,2,2)
ke <- c(2,3,2,2,4,3,4,5,3,1)
dat <- data.frame(cr, ex, ky, ke)
青木先生の関数
source("http://aoki2.si.gunma-u.ac.jp/R/src/disc.R", encoding="euc-jp")
source("http://aoki2.si.gunma-u.ac.jp/R/src/geneig.R", encoding="euc-jp")
source("http://aoki2.si.gunma-u.ac.jp/R/src/candis.R", encoding="euc-jp")
resd <- disc(dat[2:4], dat[1]) # 線形判別の関数
resc <- candis(dat[2:4], dat[1]) # 正準判別の関数
# MASSパッケージのlda関数
library(MASS)
ldares <- lda(cr~., data=dat)
ldares # 基本的に教科書の結果とは正負が逆になる。以下も同様
# 切片は自分で計算
apply(ldares$means %*% ldares$scaling, 2, mean)
# 分解するとこんな感じ
gm <- ldares$mean
cf <- ldares$scaling
x <- gm %*% cf
apply(x, 2, mean)
# Wilks' lambda
summary(manova(cbind(ex, ky, ke)~cr, data=dat), test="Wilks")
# 複数の連続変数→カテゴリカル(2値) の一般線形モデルによる分析が判別分析で、カテゴリカル→複数の連続変数 の一般線形モデルの分析がmanovaにあたる
# 判別得点など
prres <- predict(ldares)
# 判別結果
prres$class
cr # 元の群
(x <- data.frame(prres$class, cr)) # 並べてみる
t(table(x)) # 判別結果と所属群
# 事後確率
prres$posterior
# 判別得点
prres$x
## 表15.1 (B) p150
data.frame("判別得点"=prres$x[,1], "判別結果"=prres$class)
# 交差妥当性
ldacv <- lda(cr~., data=dat, CV=T)
pcr <- predict(ldares)$class
cvpcr <- ldacv$class
table(pcr, cr)
# 判別率・誤判別率
cprp(table(pcr, cr))$pt # prop.tableでもいい？
table(cvpcr, cr) # 交差妥当性の確認
cprp(table(cvpcr, cr))$pt
# http://psy.isc.chubu.ac.jp/~oshiolab/teaching_folder/datakaiseki_folder/07_folder/da07_02.html
cr <- c(1,1,1,0,1,0,0,1,0,0)
ex <- c(4,5,3,2,4,3,2,4,3,1)
ky <- c(3,4,4,5,4,2,3,5,2,2)
ke <- c(2,3,2,2,4,3,4,5,3,1)
dat <- data.frame(cr, ex, ky, ke)
青木先生の関数
source("http://aoki2.si.gunma-u.ac.jp/R/src/disc.R", encoding="euc-jp")
source("http://aoki2.si.gunma-u.ac.jp/R/src/geneig.R", encoding="euc-jp")
source("http://aoki2.si.gunma-u.ac.jp/R/src/candis.R", encoding="euc-jp")
resd <- disc(dat[2:4], dat[1]) # 線形判別の関数
resc <- candis(dat[2:4], dat[1]) # 正準判別の関数
# MASSパッケージのlda関数
library(MASS)
ldares <- lda(cr~., data=dat)
ldares # 基本的に教科書の結果とは正負が逆になる。以下も同様
# 切片は自分で計算
apply(ldares$means %*% ldares$scaling, 2, mean)
# 分解するとこんな感じ
gm <- ldares$mean
cf <- ldares$scaling
x <- gm %*% cf
apply(x, 2, mean)
# Wilks' lambda
summary(manova(cbind(ex, ky, ke)~cr, data=dat), test="Wilks")
# 複数の連続変数→カテゴリカル(2値) の一般線形モデルによる分析が判別分析で、カテゴリカル→複数の連続変数 の一般線形モデルの分析がmanovaにあたる
# 判別得点など
prres <- predict(ldares)
# 判別結果
prres$class
cr # 元の群
(x <- data.frame(prres$class, cr)) # 並べてみる
t(table(x)) # 判別結果と所属群
# 事後確率
prres$posterior
# 判別得点
prres$x
## 表15.1 (B) p150
data.frame("判別得点"=prres$x[,1], "判別結果"=prres$class)
# 交差妥当性
ldacv <- lda(cr~., data=dat, CV=T)
pcr <- predict(ldares)$class
cvpcr <- ldacv$class
table(pcr, cr)
# 判別率・誤判別率
cprp(table(pcr, cr))$pt # prop.tableでもいい？
table(cvpcr, cr) # 交差妥当性の確認
cprp(table(cvpcr, cr))$pt
# 太郎丸 博 (2005). 人文・社会科学のためのカテゴリカル・データ解析入門 ナカニシヤ出版
## p135
mat <- matrix(c(36,13,5,10,2,8), nr=2, dimnames=list(hus=c("jh",
"cu"), wif=c("wjh", "wc", "wu")))
xtb<- as.table(mat)
xtb
library(MASS)
lgres <- loglm(formula=~hus+wif, data=xtb)
coef(lgres) # 係数
fitted(lgres) # 期待値
residuals(lgres) # 標準化残差
(xtb-fitted(lgres))/ sqrt(fitted(lgres)) # ピアソンタイプの残差？
xtb # 観測値
lgres # 検定
# 組み込みのloglin関数
linres <- loglin(xtb, margin=list(c(1), c(2)), param=T, fit=T)
linres$para
linres$fit
linres$lrt # 尤度比
linres$pearson # カイ二乗
linres$df # 自由度
# 検定
1-pchisq(linres$lrt, linres$df)
1-pchisq(linres$pearson, linres$df)
## p135
mat <- matrix(c(36,13,5,10,2,8), nr=2, dimnames=list(hus=c("jh",
"cu"), wif=c("wjh", "wc", "wu")))
xtb<- as.table(mat)
xtb
library(MASS)
lgres <- loglm(formula=~hus+wif, data=xtb)
coef(lgres) # 係数
fitted(lgres) # 期待値
residuals(lgres) # 標準化残差
(xtb-fitted(lgres))/ sqrt(fitted(lgres)) # ピアソンタイプの残差？
xtb # 観測値
lgres # 検定
# 組み込みのloglin関数
linres <- loglin(xtb, margin=list(c(1), c(2)), param=T, fit=T)
linres$para
linres$fit
linres$lrt # 尤度比
linres$pearson # カイ二乗
linres$df # 自由度
# 検定
1-pchisq(linres$lrt, linres$df)
1-pchisq(linres$pearson, linres$df)
基準関連妥当性
併存妥当性 不安の高い人は抑うつも高いはずだ的なやつ
予測妥当性 アウトカム変数の予測
構成概念妥当性
収束妥当性 同じような測定と相関する。他の不安傾向尺度との相関
弁別妥当性 違う測定と相関しない
内容妥当性 ぱっと見て内容がそれっぽいか
表面妥当性 パッと見て項目がそれっぽいか
因子妥当性 因子構造が予測あるいは先行研究どおりか
交差妥当性 別サンプル間での予測誤差が小さい
おおむね収束妥当性と弁別妥当性を満たしているとOKとなる
実は表面妥当性を満たしていると他が問われないことも多い。実験心理学の測定は操作的に定義されることもありほぼ全てがそうだろう
併存妥当性 不安の高い人は抑うつも高いはずだ的なやつ
予測妥当性 アウトカム変数の予測
構成概念妥当性
収束妥当性 同じような測定と相関する。他の不安傾向尺度との相関
弁別妥当性 違う測定と相関しない
内容妥当性 ぱっと見て内容がそれっぽいか
表面妥当性 パッと見て項目がそれっぽいか
因子妥当性 因子構造が予測あるいは先行研究どおりか
交差妥当性 別サンプル間での予測誤差が小さい
おおむね収束妥当性と弁別妥当性を満たしているとOKとなる
実は表面妥当性を満たしていると他が問われないことも多い。実験心理学の測定は操作的に定義されることもありほぼ全てがそうだろう
