ロジスティック回帰分析
ロジスティック回帰分析（2値）
#データの読み込み library(readxl) dat1 <- read_excel("C:/Users/2ken/Desktop/Umeda3.xlsx") View(dat1) dat1=na.omit(dat1) #基礎集計 library(psych) describe(dat1) #データ型の確認と変更 str(dat1) dat1$graf=as.factor(dat1$graf) dat1$subject=as.factor(dat1$subject) #ロジスティック回帰 res1=glm(answer~graf+subject,data = dat1,family = "binomial") summary(res1) exp(res1$coefficients) #決定係数 library(BaylorEdPsych) PseudoR2(res1) #ロジスティック回帰の適合度の検定 library(ResourceSelection) hoslem.test(res1$y, res1$fitted.values)
メタ分析（変量効果モデル）
データの読み込み
library(metafor) library (MAc) library (MAd) dat=read.csv("csv_all2.csv") summary(dat)
変量効果モデル
res=rma(yi,vi,data=dat) res forest(res)
調整変数ごとの平均効果量
macat(yi,vi,mod = subject,method ="random",data = dat) macat(yi,vi,mod = grade,method ="random",data = dat) macat(yi,vi,mod = century,method ="random",data = dat) macat(yi,vi,mod = Measure,method ="random",data = dat) macat(yi,vi,mod = purpose,method ="random",data = dat) macat(yi,vi,mod = technology,method ="random",data = dat)
funnel plot
#funnel plot funnel(res) regtest(res,model = "rma")
fail-safe N
#fail-safe N FSN=(res$zval*11)^2/2.706-11 FSN fsn (y = res$yi, v = res$vi) plot(x=dat$year,y=dat$yi)
trim and fill
#trim and fill taf2=trimfill(res) funnel(taf2)
p.rep
#p.rep p.rep.old=pnorm(pnorm(1-0.0008)/sqrt(2)) p.rep.old #Lilleen(2005) t=res$b/res$se p.rep=pt(abs(t)/sqrt(2),df=9) p.rep #Lecoutre(2010)
多重比較（ボンフェローニ）
群間比較（多重比較，ボンフェローニ）
#学校ごとの平均算出 aggregate(X1.all~school,data=dat4,FUN = mean) #多重比較 pairwise.t.test(dat4$X1.all, dat4$school, p.adj = "bonf")
変数関係の視覚化（ペアプロット）
GGallyパッケージ
library(ggplot2) library(GGally) dat6=subset(dat4,select = c(school,think1_pre,think1_post,eap_pre,eap_post,X4_8,Q2.8)) colnames(dat6) <- c("学校","思考プレ","思考ポスト","態度プレ","態度ポスト","理科好きプレ","理科好きポスト") ggpairs(dat6, aes_string(colour="学校", alpha=0.5))
psychパッケージ
pairs.panels(dat)
段階反応モデルにおける特性値θのベイズ推定
困難度・識別力パラメータの読み込み
#Q7,9,10,11,12,13,14 雲財中村（2018）に基づくEAP推定 Extrmt1=c(-2.09,-1.55,-1.99,-1.86,-2.35,-1.48,-1.38) Extrmt2=c(-0.79,-0.55,-0.90,-0.78,-1.35,-0.49,-0.42) Extrmt3=c(0.21,0.43,0.07,0.36,-0.43,0.60,0.64) Extrmt4=c(1.51,1.37,0.96,1.47,0.74,1.47,1.54) Dscrmn=c(1.35,2.49,1.95,1.72,1.38,2.46,2.71) para2=data.frame(Extrmt1=Extrmt1,Extrmt2=Extrmt2,Extrmt3=Extrmt3,Extrmt4=Extrmt4,Dscrmn=Dscrmn)
特性値θの推定
library(irtoys) library(ltm) grm.theta(dat2, a=para2[,4],bc=para2[,c(1:4)], D=1.0, method = "ML") #EAP推定用の関数を定義 grm.eap <- function(X, a, bc, D=1.7, mu=0, sigma=1, n=50){ X<-as.matrix(X); a<-as.vector(a); bc<-as.matrix(bc) ns <- nrow(X); ni<- ncol(X); nc <- ncol(bc) qp = normal.qu(n,mu=mu,sigma=sigma)$quad.points qw = normal.qu(n,mu=mu,sigma=sigma)$quad.weights o=sapply(1:ns, function(i) eap.inside(x=X[i, ],qp,qw,a,bc,D,nc,ni,n),USE.NAMES=F) rownames(o) = c("est", "sem") return(t(o)) } eap.inside <- function (x,qp,qw,a,bc,D,nc,ni,n) { grm.LL <- function(theta){ LL01<-numeric(ni) p01<-1/(1+exp(-D*(a %*% t(rep(1,nc)))*(matrix(theta,ni,nc)-bc))) #ni*nc p02<- cbind( rep(1,ni), p01, rep(0,ni)) #ni*(nc+2) for (j in 1:ni){ LL01[j]<-p02[j,x[j]+1]-p02[j,x[j]+2]} LL02<-sum(log(LL01)) return(LL02) } ll<-numeric(n) for (i in 1:n) { ll[i]<-grm.LL(qp[i]) } wl = exp(ll) * qw swl = sum(wl) xx = sum(wl * qp)/swl dev = qp - xx sem = sqrt(sum(wl * dev * dev)/swl) return(c(xx, sem)) } cc3<-grm.eap(dat2,a=para2[,5],bc=para2[,c(1:4)],D=1.0,mu=0,sigma=1,n=50) #ロジスティック計量 dat3=cbind(dat2,cc3) write.table(dat3, file="eap.txt") #書き出し
項目反応理論（段階反応モデル，ベイズ推定）
データの読み込み
library(readr) SNFC_1 <- read_csv("SNFC_6.csv") View(SNFC_1) dat <- data.frame(SNFC_1[,8:22]) #特定の項目だけ取り出し dat=dat-1 dat2=subset(dat,select=c("Q1","Q2","Q3","Q4","Q6","Q7","Q8","Q9","Q10","Q11","Q12","Q13","Q14","Q15"))
IRTによる尺度構成
局所独立の検討
library(mirt) x <- mirt(dat2, 1) residuals(x) Theta <- fscores(x) round(residuals(x, type = 'Q3', Theta=Theta),2)
単調性の確認
library(mokken) mono.list<-check.monotonicity(dat2) summary(mono.list) (plot(mono.list,curves="ISRF", ci=TRUE,alpha = .05)) coefH(dat2)
段階反応モデルの実行
基本情報
library(ltm) descript(dat2,chi.squared = T)
特性反応曲線
par(mfrow=c(1,1)) par(ps=24) par(mar=c(5, 5, 5, 2)) plot(fitirt, legend=T, cx="left",main = "",xlab="特性値",ylab="確率",lwd=3)
テスト情報関数
par(mfrow=c(1,1)) par(ps=24) par(mar=c(5, 5, 5, 2)) plot(fitirt, type = "IIC", items = 0, plot = T,main="",lwd=3,xlab="特性値",ylab="情報量") par(mfrow=c(1,1)) vals <- plot(fitirt, type = "IIC", items = 0, plot = T) plot(vals[, "z"], 1 / sqrt(vals[, "test.info"]), type = "l", lwd = 3,xlab="特性値",ylab="標準誤差", main = "")
項目情報関数
par(mfrow=c(1,1)) #グラフ出力デフォルトに戻す plot(fitirt, type = "IIC", plot = T, legend=T, cx="left") information(fitirt,c(-3.1,3.1))
個人の能力値
fitt=factor.scores(fitirt, resp.pattern=dat2) SNFC_1$theta=fitt$score.dat$z1 View(SNFC_1)
因子得点と能力値θのプロット
plot(SNFC_1$SNFC,SNFC_1$theta) cor(SNFC_1[,23:29],use="complete.obs") #リストワイズ削除
水平テストの構築
#データの修正 typeA dat4=subset(dat,select=c("Q8","Q4","Q2","Q1","Q15","Q6","Q3")) #段階反応モデルの実行 (fitirtA=grm(dat4,Hessian = T)) #テスト情報関数 par(mfrow=c(1,1)) #グラフ出力デフォルトに戻す plot(fitirtA, type = "IIC", items = 0, plot = T,main="テスト情報曲線",xlim=c(-4,4), ylim=c(0,14), ylab="", lty=2) #データの修正 typeB dat5=subset(dat,select=c("Q14","Q13","Q9","Q11","Q7","Q10","Q12")) #段階反応モデルの実行 (fitirtB=grm(dat5,Hessian = T)) #テスト情報関数 par(new=T) plot(fitirtB, type = "IIC", items = 0, plot = T,main="",xlim=c(-4,4), ylim=c(0,14))
段階反応モデル verベイズ
呼び出し
require(rstan) options(max.print=400000000)
stanコード格納
stancode <- "
data{
int ni;
int nj;
int nc;
real D;
int<lower=1,upper=6> y[ni,nj];
}
parameters{
vector<lower=0,upper=5>[nj] a;
ordered[nc-1] ba[nj];
vector<lower=-4,upper=4>[ni] theta;
}
transformed parameters{
real b[nj,nc];
vector<lower=0,upper=1>[nc-1] pa[ni,nj];
simplex[nc] p[ni,nj];
for (j in 1:nj){
for (c in 1:nc){
if (c ==1){
b[j,c]=ba[j,c];
}else if (c ==nc){
b[j,c]=ba[j,c-1];
}else{
b[j,c]=(ba[j,c-1]+ba[j,c])/2;
}
}
}
for (i in 1:ni){
for (j in 1:nj){
for (c in 1:nc-1){
pa[i,j,c]= 1/(1+exp(-D*a[j]*(theta[i]-ba[j,c])));
}
}
}
for (i in 1:ni){
for (j in 1:nj){
for(c in 1:nc){
if (c==1){
p[i,j,c]=1-pa[i,j,c];
}else if(c==nc){
p[i,j,c]= pa[i,j,c-1];
}else{
p[i,j,c]= pa[i,j,c-1]-pa[i,j,c];
}
}
}
}
}
model{
for (i in 1:ni){
theta[i]~normal(0,1);
for (j in 1:nj){
y[i,j]~categorical(p[i,j]);
}
}
for (j in 1:nj){
a[j]~lognormal(0,sqrt(0.5));
for (c in 1:nc-1){
ba[j,c]~normal(0,2);
}
}
}
generated quantities{
real bg[nj,nc];
bg=b;
}このstanコードは豊田(2017)を引用
実行
u_ori<-dat2[1:14]+1 ni<-1508 nj<-14 nc<-5 D<-1 data_ori<-list(y=u_ori,nj=nj,ni=ni,nc=nc,D=D) par<-c("theta","ba","a","b") war<-1000 ite<-21000 see<-1234 dig<-3 cha<-5 fit_ori<-stan(file=stancode,model_name=stancode,data=data_ori,pars=par,verbose=F,seed=see,chains=cha,warmup=war,iter=ite) print(fit_ori,pars=par,digits_summary=dig) sink("./result1.txt") #保存 print(fit_ori,pars=par,digits_summary=dig) sink()
２群の平均値の差の検討（t検定，ベイズ推定）
２群の平均値の差の検討（性差の検討を例に）
データの格納
x<-SNFC_1$theta[SNFC_1$sex=="2"] #第1群の指定 y<-SNFC_1$theta[SNFC_1$sex=="1"] #第2群の指定 #平均値の高い方をyにしておくと後で解釈しやすい nx <- length(x) ny <- length(y) data <- list(N1=nx,N2=ny,x1=x,x2=y,c=0)
対応のないt検定
var.test(x, y) #等分散性の検定.p < 0.05 => 等分散でない．p >= 0.05 => 等分散である t.test(x,y,var.equal=T) #t検定.等分散でない時はvar.equal=F(Welch検定)
対応のないt検定 ver.ベイズ
stanコード格納
pairedF <- '
data{
int<lower=0> N1;
int<lower=0> N2;
real x1[N1];
real x2[N2];
real c;
}
parameters{
real mu1;
real mu2;
real<lower=0> sigma1;
real<lower=0> sigma2;
}
model{
x1 ~ normal(mu1,sigma1);
x2 ~ normal(mu2,sigma2);
}
generated quantities{
real delta;
real d_overC;
real cohen_d;
delta = mu2 - mu1;
d_overC = (delta>c? 1:0);
cohen_d = delta /((sigma1^2*(N1-1)+sigma2^2*(N2-1)) /((N1-1)+(N2-1)))^0.5;
}
'require(rstan) rstan_options(auto_write = TRUE) options(mc.cores = parallel::detectCores()) fit<-stan(model_code=pairedF,data=data,chains=5,warmup=1000,iter=21000) print(fit,digits_summary=3)
