source("http://www.matsuka.info/univ/course_folder/cutil.R")
dat<-read.csv("http://www.matsuka.info/data_folder/dktb312.csv")
dat$method=factor(dat$method, levels(dat$method)[c(1,3,4,2)])
dat.aov<-aov(result~method, data=dat)
new.alpha<-adj.alpha(5,0.05)
cu.bonF1F(dat.aov,dat,c(-3,1,1,1),new.alpha)
cu.bonF1F(dat.aov,dat,c(-1,1,0,0),new.alpha)
cu.bonF1F(dat.aov,dat,c(-1,0,1,0),new.alpha)
cu.bonF1F(dat.aov,dat,c(-1,0,0,1),new.alpha)
new.f<-adj.f(4,3,28,0.05)
cu.scheffe1F(aov.obj,dat,c(-1,0,1,0))
cu.scheffe1F(aov.obj,dat,c(-1,0,0,1))
dat<-read.csv("http://matsuka.info/univ/course_folder/datW03R.txt")
interaction.plot(dat$gender,
dat$affil,
dat$shoesize,
pch=c(20,20),
col=c("skyblue","orange"),
xlab="gender", ylab="shoesize",
lwd=3,type='b',cex=2,
trace.label="Affiliation")
dat.aov=aov(shoesize~gender*affil, data=dat)
dat.aov.sum=summary(dat.aov)
means<-tapply(dat$shoesize, list(dat$gender,dat$affil), mean)
SS_gen_CS<- 5*(means[2,1]^2 + means[1,1]^2 -0.5*sum(means[,1])^2) # SS_gender CS
SS_gen_PS<- 5*(means[2,2]^2 + means[1,2]^2 -0.5*sum(means[,2])^2) # SS_gender PS
dat.aov.sum=summary(dat.aov) # ANOVA table
MSe=dat.aov.sum[[1]][4,3] # MSE from ANOVA table or MSe=0.62
dfE=dat.aov.sum[[1]][4,1] # DF for error or dfE=16
dfG=1 # DF for gender
F_gen_CS=(SS_gen_CS/dfG)/MSe # F-value for gender effect given CS
F_gen_PS=(SS_gen_PS/dfG)/MSe # F-value for gender effect given PS
P_gen_CS=1-pf(F_gen_CS,1,dfE) # p-value for gender effect given CS
P_gen_PS=1-pf(F_gen_PS,1,dfE) # p-value for gender effect given PS
SS_affil_F<- 5*(means[1,1]^2+means[1,2]^2-0.5*sum(means[1,])^2) #SS_affil | F
SS_affil_M<- 5*(means[2,1]^2+means[2,2]^2-0.5*sum(means[2,])^2) #SS_affil | M
dfA=1 # DF for affil
F_affil_F=SS_affil_F/dfA/MSe # F-value for affiliation effect | F
F_affil_M=SS_affil_M/dfA/MSe # F-value for affiliation effect | M
P_affil_F=1-pf(F_affil_F,1,dfE) # p-value for affiliation effect | F
P_affil_M=1-pf(F_affil_M,1,dfE) # p-value for affiliation effect | M
tsme = CRF.tsme(dat.aov, dat)
dat<-read.csv("http://www.matsuka.info/data_folder/dktb321.txt")
interaction.plot(dat$duration,dat$method,dat$result, pch=c(20,20), col=c("skyblue","orange"), ylab="score", xlab="Duration",lwd=3,type='b',cex=2,trace.label="Method")
mod1=aov(result~method+duration,data=dat)
mod1.sum=print(summary(mod1))
mod2=aov(result~method*duration,data=dat)
mod2.sum=print(summary(mod2))
qv=qtukey(0.95,DFd+1,DFe)
hsd=qv*(sqrt(MSe/5))
means<-tapply(dat$result,list(dat$method,dat$duration),mean)
print(diffM<-outer(means[1,],means[1,],"-"))
abs(diffM)>hsd
CRF.tsme(mod2, dat)
Category Archives: データ解析
データ解析基礎論 W10
f=c(24.1,23.9,24.4,24.4,23.5)
m=c(25.6,26.1,25.8,25.9,26)
G.mean=mean(c(f,m))
ss.total=sum((c(f,m)-G.mean)^2)
ss.tr=sum((mean(f)-G.mean)^2)*5+sum((mean(m)-G.mean)^2)*5
ss.error=sum((f-mean(f))^2)+sum((m-mean(m))^2)
df.tr = 2 - 1
df.error = (5-1)*2
ms.tr = ss.tr /df.tr
ms.error = ss.error /df.error
dat<-data.frame(ssize=c(f,m),gender=c(rep("f",5),rep("m",5)))
dat.aov<-aov(ssize ~ gender,data=dat)
summary(dat.aov)
dat<-read.table("http://www.matsuka.info/data_folder/aov01.txt")
dat.aov<-aov(shoesize~club, dat)
summary(dat.aov)
qT<-qtukey(0.95, 3, 67)
HSD<-qT*sqrt((2.243*(1/23+1/24))/2)
means<-tapply(dat$shoesize,dat$club,mean)
abs(outer(means,means,"-"))
abs(outer(means,means,"-"))>HSD
TukeyHSD(dat.aov)
dat<-read.csv("http://www.matsuka.info/data_folder/dktb312.csv")
dat$method=factor(dat$method, levels(dat$method)[c(1,3,4,2)])
dat.aov<-aov(result~method, data=dat)
summary(dat.aov)
dat.means<-tapply(dat$result,dat$method,mean)
new.alpha = 1-(1-0.05)^(1/5)
cont=c(-3,1,1,1)
bunshi=sum(cont*dat.means)
bunbo=sqrt(5.29*(sum((cont^2)/8)))
t.value=bunshi/bunbo
2*(1-pt(t.value,28))
データ解析基礎論A W09
dat <- read.csv("http://www.matsuka.info/data_folder/tdkReg01.csv")
dat.lm1 <- lm(sales~price,data = dat)
summary(dat.lm1)
dat.lm2 <- lm(sales~price+material,data = dat)
summary(dat.lm2)
dat.lm01<-lm(sales~price, data=as.data.frame(scale(dat)))
plot(dat.lm01,which=1)
plot(dat.lm01,which=2)
norm.vars=rnorm(300)
qqnorm(norm.vars)
qqline(norm.vars,col='red',lwd=2)
unif.vars=runif(300)
qqnorm(unif.vars)
qqline(unif.vars,col='green',lwd=2)
par(mfrow=c(2,2))
plot(dat.lm01)
dat<-read.csv("http://www.matsuka.info/data_folder/dktb312.csv")
dat$method=factor(dat$method, levels(dat$method)[c(1,3,4,2)])
dat.lm<-lm(result~method,data=dat)
plot(dat.lm)
dat <- read.csv("http://www.matsuka.info/data_folder/tdkReg01.csv")
dat2 <- rbind(dat,c(0,0,100,250))
tail(dat2)
plot(dat2[1:50,],col=c(rep("black",50),"red"),pch=20,cex=3)
plot(dat2[1:50,],pch=20,cex=3)
dat.lm<-lm(sales~design,data=dat)
summary(dat.lm)
dat2.lm<-lm(sales~design,data=dat2)
summary(dat2.lm)
dat<-read.csv("http://www.matsuka.info/data_folder/tdkReg02.csv")
dat.lm<-lm(sales~., data=dat)
install.packages("DAAG")
library(DAAG)
vif(dat.lm)
dat<-read.csv("http://www.matsuka.info/data_folder/waa07_2.csv")
poly.lm<-lm(grade~study+study.sq, data= dat)
dat<-read.table("http://www.matsuka.info/data_folder/tdkPATH01.txt",header=TRUE)
plot(dat,pch=20,cex=2)
dat.lm1<-lm(absence~interest,dat)
dat.lm2<-lm(study~interest+knowledge,dat)
dat.lm3<-lm(grade~knowledge+study+absence,dat)
dat.STD<-as.data.frame(scale(dat))
summary(lm(absence~interest,dat.STD))
summary(lm(study~interest+knowledge,dat.STD))
summary(lm(grade~study+absence+knowledge,dat.STD))
データ解析基礎論 W08
dat <- read.csv("http://www.matsuka.info/data_folder/tdkReg01.csv")
dat.lm <- lm(sales~.,data = dat)
summary(dat.lm)
anova(dat.lm)
sse = anova(dat.lm)[[2]][4]
ssr = sum(anova(dat.lm)[[2]][1:3])
df.err = anova(dat.lm)[[1]][4]
df.model = sum(anova(dat.lm)[[1]][1:3])
F.val = (ssr/3) / (sse/46)
1 - pf(F.val, df.model, df.err)
dat<-read.csv("http://peach.l.chiba-u.ac.jp/course_folder/dktb312.csv")
tapply(dat$result,dat$method,mean)
dat2<-data.frame(result=dat$result,
c1=c(rep(0,8),rep(1,8),rep(0,16)),
c2=c(rep(0,16),rep(1,8),rep(0,8)),
c3=c(rep(0,24),rep(1,8)))
dat2.lm<-lm(result~.,data=dat2)
dat3<-data.frame(result=dat$result,
c1=c(rep(-3,8), rep(1,24)),
c2=c(rep(0,8),rep(-2,8),rep(1,16)),
c3=c(rep(0,16),rep(-1,8),rep(1,8)))
dat3.lm<-lm(result~c1+c2+c3,data=dat3
dat<-read.csv("http://www.matsuka.info/data_folder/dktb321.txt")
plot(dat$result~dat$duration,data=dat[dat$method=="method.X",])
result<-dat$result[dat$method=="method.X"]
CL=c(rep(-3,5),rep(-1,5),rep(1,5),rep(3,5))
CQ=c(rep(-1,5),rep(1,5),rep(1,5),rep(-1,5))
CC=c(rep(-3,5),rep(1,5),rep(-1,5),rep(3,5))
trend.lm<-lm(result~CL+CQ+CC)
summary(trend.lm)
dat.x = dat[dat$method=="method.X",]
trend.lm2 <- lm(result~duration, contrasts=list(duration = "contr.poly"),data=dat.x)
dat.x$result[16:20]=dat.x$result[16:20]-3
plot(tapply(dat.x$result,dat.x$duration,mean),
pch=20,col="red",lwd=3, type="o",cex=3,ylab="mean",xlab="time")
trend.lm2 <- lm(result~duration, contrasts=list(duration = "contr.poly"),data=dat.x)
subj.gender = c(rep("male",30),rep("female",30))
pic.gender = rep(c(rep("male",15),rep("female",15)),2)
# 15m-m, 15m-f, 15f-m, 15f-f
eye.fix = c(round(rnorm(15,50,5)),
round(rnorm(15,70,5)),
round(rnorm(15,65,5)),
round(rnorm(15,50,5)))
datE = data.frame(eye.fix = eye.fix, subj.gender= subj.gender,pic.gender=pic.gender)
interaction.plot(datE$subj.gender, datE$pic.gender, datE$eye.fix,
trace.label = "Gender of Stimuli", xlab = "Gender of Participants",
ylab = "Number of fixation",lwd = 3,type = "b",pch = c(17,20),cex =3,
col = c('skyblue','coral'),legend =T)
datE.lm <- lm(eye.fix~subj.gender*pic.gender,data=datE)
summary(datE.lm)
dat<-read.csv("http://www.matsuka.info/data_folder/tdkReg01.csv")
dat.lm01<-lm(sales~price, data=as.data.frame(scale(dat)))
plot(dat.lm01,which=1)
plot(dat.lm01,which=2)
par(mfrow=c(2,2))
plot(dat.lm01)
dat<-read.csv("http://peach.l.chiba-u.ac.jp/course_folder/dktb312.csv")
levels(dat$method)
dat$method=factor(dat$method, levels(dat$method)[c(1,3,4,2)])
dat.lm<-lm(result~method,data=dat)
par(mfrow=c(2,2))
plot(dat.lm)
dat<-read.csv("http://www.matsuka.info/data_folder/tdkReg02.csv")
plot(dat)
dat.lm<-lm(sales~., data=dat)
データ解析基礎論A W07
dat<-read.csv("http://www.matsuka.info/data_folder/hwsk8-17-6.csv")
dat.lm <- lm(ani~otouto, data=dat)
summary(dat.lm)
t.test(dat$ani, dat$otouto, var.equal=T)
dat2 <- data.frame(score=c(dat$ani,dat$otouto),order=c(rep("ani",10),rep("otouto",10)))
plot(dat2$score~as.numeric(dat2$order),pch=20,xlab="order",
ylab="score",xlim=c(0.5,2.5),cex=2,xaxt="n")
axis(1,c(1,2),c("ani","otouto"))
dat2.lm<-lm(score~order,data=dat2)
abline(dat2.lm,col='red',lwd=3)
dat.D = dat$ani - dat$otouto
boxplot(dat.D,col="skyblue",ylab="Difference")
t.test(dat.D)
plot(dat.D~rep(1,10),pch=20,xlab="",ylab="Difference",cex=3)
dat.D.lm<-lm(dat.D~1)
abline(dat.D.lm,col='red',lwd=3)
dat<-read.csv("http://www.matsuka.info/data_folder/tdkReg01.csv")
データ解析基礎論A 06
ssize = c(24,25,26,23.5,25,27,24,22,27.5,28)
t.test(ssize, mu=24)
dat<-read.csv("http://www.matsuka.info/data_folder/datWA01.txt")
t.test(dat$h[dat$gender=="M"],mu=171)
t.test(dat$h[dat$gender=="F"],mu=158)
A=c(12,19,10,10,14,18,15,11,16)
B=c(15,20,16,14,17,16,12,12,19)
d=A-B
t.test(d, mu=0)
X1=c(78,70,66,76,78,76,88,76)
X2=c(76,72,60,72,70,72,84,70)
t.value=(mean(X1)-mean(X2))/sqrt((var(X1)+var(X2))/8)
t.test(X1,X2,var.equal=T)
t.test(X1,X2)
dat<-read.csv("http://www.matsuka.info/data_folder//mobH.csv")
###
may = dat$height[dat$mob=="may"]
dec = dat$height[dat$mob=="dec"]
dat3 = data.frame(height = c(may,dec),
mob = c(rep('may',length(may)),rep('dec',length(dec))))
boxplot(height~mob, data=dat3,col=c('skyblue','orange'),
ylab = "height", xlab ="mob")
spring = dat$height[dat$mob =="march"|
dat$mob =="April"|
dat$mob =="may"]
fall = dat$height[dat$mob =="sep"|
dat$mob =="oct"|
dat$mob =="nov"]
dat4 = data.frame(height = c(spring,fall),
season = c(rep('spring',length(spring)),
rep('fall',length(fall))))
boxplot(height~season, data=dat4, col=c('skyblue','orange'),
ylab = "height", xlab ="season")
データ解析基礎論a W05
dat<-read.csv("http://www.matsuka.info/data_folder/datWA01.txt")
# Male
mean.M <-mean(dat$h[dat$gender=="M"])
stddev = 10
n.M = length(dat$h[dat$gender=="M"])
z.value=(mean.M-171)/(sqrt(stddev^2/n.M))
(1-pnorm(abs(z.value)))*2
# Female
mean.F <-mean(dat$h[dat$gender=="F"])
stddev = 5
n.F = length(dat$h[dat$gender=="F"])
z.value=(mean.F-158)/(sqrt(stddev^2/n.F))
(1-pnorm(abs(z.value)))*2
データ解析基礎論A 練習課題WA04
データ解析基礎論A 練習課題（提出不要）
練習課題WA04.1
マーケティングリサーチで商品Aと商品Bをレーティングしてもらい、そしてどちらか一方をレーティングに基づいて選択してもらった結果が以下のデータです。
S1 S2 S3 S4 S5 S6 S7 S8 S9
A 12 19 10 10 14 18 15 11 16
B 15 20 16 14 17 16 12 12 19
-----------------------------------
B B B B B A A B B
つまり、９人中７人が商品Bを選択しました。
商品Aと商品Bが同じ割合で選択されると仮定した場合、以下の確率を求めてください：
A: ９人中７人以上がBを選択する確率を計算してください。
B: ９人中７人以上が一方の商品を選択する確率を計算してください。
練習課題WA04.2
ある通りの平日午前10~11時の1時間の車の交通量は800台で標準偏差が40台だとします。また、この通りの同じ曜日の同じ時間帯の交通量は正規分布に従うと仮定して、以下の数値を求めてください；
A) 交通量が750台から850台になる確率
B）交通量が700台以下となる確率
C）平均値の800を中心に95%の範囲（下限、上限）
D）99%の上限（99%の確率でX台以下となるようなX）
練習課題WA04.3
一般成人の血液100cc中の赤血球の量の平均値が150ユニットで標準偏差が15とします。また、その分布は正規分布に従うと仮定します。運動部Xの10名の赤血球の量を計測したら、165ユニットでした。
A) 運動部Xの赤血球の量は一般成人と同等であると仮定した場合、部員の赤血球の平均値が165ユニット以上となる確率を求めてください。
B）運動部Xの赤血球の量は一般成人と同等であると仮定したが、Aで得られた確率がどの程度だとこの仮説が妥当だと思いますか？
2018 DAA W03 Data Viz & Descrptive statistics
dat<-read.table("http://www.matsuka.info/data_folder/aov01.txt")
par(mfrow=c(1,2))
hist(dat[dat$gender=='F',]$h, main="Dist. of Height for Female Participants",
xlab="Height", xlim=c(140,190), probability=T)
dens.F = density(dat[dat$gender=='F',]$h)
lines(dens.F, col='blue',lwd=2)
hist(dat[dat$gender=="M",]$h,
main="Dist. of Height for Male Participants",
xlab="Height", xlim=c(140,190), probability=T,ylim=c(0,0.08))
dens.M = density(dat[dat$gender=='M',]$h)
lines(dens.M, col='green', lwd=2)
par(mfrow=c(1,1))
plot(dens.F,col='blue',lwd=2, ylab='density', xlim=c(140,190),
main="Dist. of Height by gender",xlab='Height')
lines(dens.M,col='green',lwd=2)
legend("topleft", c('Female','Male'), col=c('blue','green'), cex=1.5,lwd=2)
text(x = 157.5, y = 0.04, 'Female', col='blue', cex=2)
text(x = 170, y = 0.04,'Male', col='green', cex=2)
plot(dat$shoesize, dat$h, main="Relationship b/w shoesize and height",
xlab = "shoesize", ylab="height", pch=19, col="red")
txt = paste("r =",round(cor(dat$shoesize,dat$h), 4))
text(22, 175, txt, cex = 1.5)
abline(h = mean(dat$h), col='blue');
abline(v = mean(dat$shoesize), col='green');
abline(lm(dat$h~dat$shoesize), lty=2, lwd=2)
plot(dat[dat$gender=='F',]$shoesize, dat[dat$gender=='F',]$h,
main="Relationship b/w shoesize and height", xlab='shoesize', ylab='height',
cex.lab=1.5, pch=19, col='blue', xlim=c(20,29), ylim=c(140,190))
points(dat[dat$gender=='M',]$shoesize,dat[dat$gender=='M',]$h,
pch = 15, col = 'green')
legend("topleft", c('Female','Male'), pch =c(19,15),
col = c('blue','green'), cex = 1.5)
dat<-read.csv("http://www.matsuka.info/data_folder/tdkReg01.csv")
plot(dat, pch=20, col='blue')
dat.pca<-read.table("http://www.matsuka.info/data_folder/tdkPCA01.txt")
plot(dat.pca, pch = rownames(dat.pca), cex = 1.7, col = 'blue')
2018 DAA W02 Data visualization
v1 = seq(-3,3,0.1)
v2 = v1^2
plot(x = v1, y = v2)
plot(v1, v2, col = 'red')
plot(v1, v2, col=“red”, pch = 20)
plot(v1, v2, col="red", pch = 20, cex = 3)
plot(v1, v2, type = "l")
plot(v1, v2,type="l",lty=4,lwd=3)
plot(v1, v2, main = "THIS IS THE TITLE",
xlab = "Label for X-axis",
ylab = "Label for Y-axis")
plot(v1, v2, main = "THIS IS THE TITLE", cex.lab = 1.5,
xlab = "Label for X-axis",ylab = "Label for Y-axis")
plot(v1, v2, main = "TITLE", xlab = "X here",ylab = "Y here",
xlim = c(-3.5, 3.5),
ylim = c(-0.5, 10))
plot(v1, v2, col = "blue", type = "o", lty = 2, pch = 19,
cex.lab = 1.5, lwd = 3, main = "Y=X*X", xlab = "X",
ylab="X*X", xlim=c(-3.5,3.5), ylim=c(-0.5, 10))
dat<- read.csv("http://www.matsuka.info/data_folder/datWA01.txt")
hist(dat$h)
hist(dat$h, breaks = 20, main = "Histogram of Height",
xlab = "Height", col = 'blue', xlim = c(140, 190))
dens<-density(dat$h); # 確率密度の算出
hist(dat$h, main = "Histogram of Height", xlab = "Height",
xlim = c(140,190), probability = T)
lines(dens, lwd = 2, col = "red", lty=2)
plot(v1, v2, col = "blue", type = "l", pch = 19, cex.lab = 1.5, lwd = 3, xlab = "X",
ylab="f(X)", xlim=c(-3.5,3.5), ylim=c(-0.5, 10))
lines(v1, v1^3, col='red',lwd = 3)
legend("bottomright", c("x^2","x^3"), col=c('blue','red'), lwd=2)
boxplot(dat$h,main="Boxplot of Height", ylab="Height", col='cyan', ylim=c(140,190))
boxplot(dat$h,main="Boxplot of Height", xlab="Height", col='orange', horizontal=T)
boxplot(dat$h ~ dat$gender,
main="Distribution of Height by Gender",
ylab="Gender", xlab="Height", col=c('blue','cyan'),
ylim=c(140,190), horizontal=T)
dat<-read.table("http://www.matsuka.info/data_folder/aov01.txt")
boxplot(dat$h ~ dat$gender + dat$affil,
main="Distribution of Height by Gender and Affiliation",
ylab="Gender x Affiliation", xlab="Height", col=c("blue","cyan","red","magenta"),
ylim=c(140,190),horizontal=T)
interaction.plot(dat$gender,
dat$affil,
dat$h,
pch=c(20,20),
col=c("skyblue","orange"),
xlab="gender", ylab="height",
lwd=3,type='b',cex=2,
trace.label="Affiliation")
par(mfrow=c(1,2))
hist(dat[dat$gender=="F",]$h, main="Dist. of Height for Female Participants", xlab="Height", xlim=c(140,190), probability=T)
dens.F = density(dat[dat$gender=='F',]$h)
lines(dens.F, col='blue',lwd=2)
hist(dat[dat$gender=="M",]$h, main="Dist. of Height for Male
Participants", xlab="Height", xlim=c(140,190), probability=T,ylim=c(0,0.08))
dens.M = density(dat[dat$gender=='M',]$h)
lines(dens.M, col='green', lwd=2)
par(mfrow=c(1,1))
plot(dens.F,col='blue',lwd=2, ylab='density', xlim=c(140,190), main="Dist. of Height by gender",xlab='Height')
lines(dens.M,col='green',lwd=2)
legend("topleft", c('Female','Male'), col=c('blue','green'), cex=1.2,lwd=2)
plot(dat$shoesize, dat$h, main="Relationship b/w shoesize and height",
xlab = "shoesize", ylab="height", pch=19, col="red")