Rを通じて統計学を学ぶ備忘録ブログ

SPSSからRに移行したい私のような人向けのR解説ブログ兼学習用備忘録。

ちょっとした備忘録のページ。

#確率密度関数=xをとる確率を返す
#xの確率は?
dnorm(0 ,0,1)
[1] 0.3989423

#累積分布関数=xをとるまでの累積確率
#xの値を取るまでの累積確率は?
pnorm(0)
[1] 0.5

#分位点=下側確率がxとなるZの値が求まる
#xの確率を取るときのZ値は?
qnorm(0.025)
[1] -1.959964
qnorm(0.975)
[1] 1.959964

pnorm(-1.96) + (1-pnorm(1.96))
[1] 0.04999579
x <- c(4, 1, 2, 3, 1, 2, 4, 2, 2, 1, 3, 1, 3, 2)
Eng <- c("A", "B", "C", "D")
y <- rep(NA, length(x))
for(i in 1:4){
  y[x==i] <- Eng[i]
}
data.frame(x = x, y = y)
set.seed(71)
dat <- matrix(floor(runif(30, min = 0,max = 10)), ncol = 3)
dat <- as.data.frame(dat)
colnames(dat) <- c("a", "b", "c")
dat

> dat
   a b c
1  3 1 6
2  5 7 2
3  3 1 7
4  2 7 2
5  3 8 1
6  9 7 4
7  6 1 4
8  8 8 3
9  3 8 2
10 4 9 8


#a列で重複している行を削除
dat[!duplicated(dat$a), ]

   a b c
1  3 1 6
2  5 7 2
4  2 7 2
6  9 7 4
7  6 1 4
8  8 8 3
10 4 9 8

#b,c列の組み合わせで重複している行を削除
dat[!duplicated(paste(dat$b, dat$c, sep = ",")), ]

   a b c
1  3 1 6
2  5 7 2
3  3 1 7
5  3 8 1
6  9 7 4
7  6 1 4
8  8 8 3
9  3 8 2
10 4 9 8
options(scipen = 5)

##想定クリック数1000回、クリック率1%の時の二項分布。
d <- as.data.frame(dbinom(x = 0:30, size = 1000, prob = 0.01))
d$x <- 0:30
colnames(d) <- c("y","x")
ggplot(data = d, aes(x = factor(x), y = y))+ geom_bar(stat = "identity")

##平均λ=1.43のポアソン分布。7時間で10cvの場合。
d <- as.data.frame(dpois(x = 0:10, lambda = 1.428))
d$x <- 0:10
colnames(d) <- c("y","x")
ggplot(data = d, aes(x = factor(x), y = y))+ geom_bar(stat = "identity")

参照元http://www.mm-lab.jp/analysis/data-analysis-start-with-excel-5th/

【オンライン広告の認知効果】

x <- 1:1000
y <- 5.980469*x^0.265342 + rnorm(x, 0, 5)
d <- data.frame(x,y)
plot(x,y)
fit <- lm(log(y) ~ log(x))
summary(fit)

result.con <- data.frame(x = d$x, conf = exp(predict(fit, interval="confidence", level = 0.90)))

ggplot(d, aes(x = x, y = y)) + geom_point() +
  geom_ribbon(data = result.con, aes(ymin=result.con$conf.lwr, ymax = result.con$conf.upr), fill = "darkgreen", alpha=0.2) +
  geom_line(data = result.con, aes(x = x, y = conf.fit)) 

参照元http://www.videor.co.jp/press/2015/150903.pdf

【重回帰分析の信頼区間作図】

m1071 <- lm(ACHIEVE ~ STUDY + FACULTY + STUDY:FACULTY, data=df1)

mean(df1$FACULTY)
dummy0 <- 0

min(df1$STUDY) #xの範囲を知るため
max(df1$STUDY) #xの範囲を知るため
xindex <- seq(4, 10, 1) #動かす変数の範囲


#動かす変数以外に平均値を放り込む
conf <- predict(m1071, 
                data.frame(STUDY = xindex, 
                           FACULTY = dummy0),
                interval="confidence") 

conf <- as.data.frame(conf)
conf$STUDY <- seq(4, 10, 1)

ggplot() + 
  geom_line(data = conf, aes(STUDY, fit)) +
  geom_ribbon(data = conf, aes(STUDY, ymin=lwr, ymax=upr), alpha=1/2) +
  geom_point(data = df1, aes(STUDY, ACHIEVE))

【重回帰分析の交互作用の±1SD】

m1071 <- lm(ACHIEVE ~ STUDY + INTEREST + STUDY:INTEREST, data=df1)

mean <- mean(df1$INTEREST)
mean_upp_1sd <- mean(df1$INTEREST) + sd(df1$INTEREST)
mean_low_1sd <- mean(df1$INTEREST) - sd(df1$INTEREST)

min(df1$STUDY) #xの範囲を知るため
max(df1$STUDY) #xの範囲を知るため
xindex <- seq(4, 10, 1) #動かす変数の範囲


#動かす変数以外に平均値を放り込む
conf <- predict(m1071, 
                data.frame(STUDY = xindex, 
                           INTEREST = mean),
                interval="confidence") 

conf <- as.data.frame(conf)
conf$STUDY <- seq(4, 10, 1)
##----------------------------
conf_upp <- predict(m1071, 
                  data.frame(STUDY = xindex, 
                             INTEREST = mean_upp_1sd),
                  interval="confidence") 
conf_upp <- as.data.frame(conf_upp)
conf_upp$STUDY <- seq(4, 10, 1)
##----------------------------
conf_low <- predict(m1071, 
                    data.frame(STUDY = xindex, 
                               INTEREST = mean_low_1sd),
                    interval="confidence") 
conf_low <- as.data.frame(conf_low)
conf_low$STUDY <- seq(4, 10, 1)
##----------------------------

ggplot() + 
  geom_line(data = conf, aes(STUDY, fit)) +
  geom_line(data = conf_upp, aes(STUDY, fit), colour="red",linetype="dashed") +
  geom_line(data = conf_low, aes(STUDY, fit), colour="blue",linetype="dashed") +
  geom_point(data = df1, aes(STUDY, ACHIEVE))

predcit関数

n <- 5
a <- 2
b <- 2

x <- rnorm(n, 10, 1)
e <- rnorm(n, 0, 1)
y <- a + b*x + e

plot(x,y)
model <- lm(y~x)
model

test <- data.frame(x = seq(min(x), max(x), 0.5))
A <- predict(object = model,
             newdata = newdata,
             se.fit = TRUE,
             interval = "confidence",
             levels = 0.95)

##predict関数による信頼区間
lines(test[,1], A$fit[,1])
lines(test[,1], A$fit[,2], col = "red", lty = "dotdash")
lines(test[,1], A$fit[,3], col = "red", lty = "dotdash")

##自作による信頼区間
test <- data.frame(x = seq(min(x), max(x), 0.5), predcit = A$fit)
test$self.lwr <- A$fit[,1] - A$se.fit * qt(0.975, A$df)
test$self.upr <- A$fit[,1] + A$se.fit * qt(0.975, A$df)
lines(test[,1], test[,5], col = "blue", lty = "dotted")
lines(test[,1], test[,6], col = "blue", lty = "dotted")

##予測区間
A$fit[,1] - sqrt(A$se.fit^2 + A$residual.scale^2)*qt(0.975, A$df)
A$fit[,1] + sqrt(A$se.fit^2 + A$residual.scale^2)*qt(0.975, A$df)
A <- predict(object = model,
             newdata = newdata,
             se.fit = TRUE,
             interval = "predicion",
             levels = 0.95)


#割合49.5%と50.5%に固定し、サンプルサイズだけを変動させた例
test1 <- prop.test(10*c(0.495,0.505), c(10,10))
test2 <- prop.test(100*c(0.495,0.505), c(100,100))
test3 <- prop.test(1000*c(0.495,0.505), c(1000,1000))
test4 <- prop.test(10000*c(0.495,0.505), c(10000,10000))
test5 <- prop.test(100000*c(0.495,0.505), c(100000,100000))
test6 <- prop.test(100000*c(0.495,0.505), c(1000000,1000000))

pval <- c(test1$p.value,
          test2$p.value,
          test3$p.value,
          test4$p.value,
          test5$p.value,
          test6$p.value)

d <- data.frame(pval = pval,
                test = paste("Q",1:length(pval), sep=""))

ggplot(d, aes(test, pval)) + geom_point() + geom_line(aes(group = 1)) +
  geom_hline(yintercept = c(0.10, 0.05, 0.001), linetype = "dashed") +
  annotate("text", label = c("p-value < 0.10", "p-value < 0.05", "p-value < 0.001"),
           x = 1, y = c(0.12, 0.07, 0.02), size = 7)

d <- read.csv("ana.csv", header = TRUE)
d$upr <- qbinom(0.975, size = d$access, prob = d$cvr)/d$access
d$lwr <- qbinom(0.025, size = d$access, prob = d$cvr)/d$access

ggplot(d, aes(month, cvr, fill = format)) + geom_bar(stat = "identity") + 
  geom_errorbar(data = d ,aes(ymin = lwr, ymax = upr, width = 0.1)) + facet_grid(.~format)

ggplot(d, aes(month, cvr, fill = format)) +
  geom_ribbon(data = d ,aes(ymin = lwr, ymax = upr, alpha = .5)) + geom_line() 

ggplot(d, aes(month, cvr, fill = format)) + geom_line() + 
  geom_line(aes(month, upr), linetype = "dashed") +
  geom_line(aes(month, upr), linetype = "dashed")

ggplot(d, aes(month, cvr, col = format)) + geom_line() + geom_point(size = 2) + 
  geom_errorbar(data = d ,aes(ymin = lwr, ymax = upr, width = .1))

http://sc1.cc.kochi-u.ac.jp/~murakami/cgi-bin/FSW/fswiki.cgi?page=%B8%FA%B2%CC%CE%CC%A1%A6%B8%A1%C4%EA%CE%CF%A1%A6%A5%B5%A5%F3%A5%D7%A5%EB%A5%B5%A5%A4%A5%BA
検定力分析では,(1) 実験を実施した後に,サンプルサイズ,効果量,有意水準(α)から,検定力(1-β)を確認する方法と,
(2) 実験を実施する前に,これまでの先行研究からわかっている(推測される)効果量,4 有意水準(α),目指している検定力(1-β)からサンプルサイズを決定する目的のものが多い 
#サンプルサイズ設計問題〜検出力0.8を保つ〜
#有意水準(α)#検出力(1-β)
#広告 A の経験的 CVR
#広告 B の期待される CVR
#α=0.05、検出力 1-β = 0.8、広告Aの経験的CVRは0.2%=0.002、広告BのCVR0.3%=0.02だったらいいな、と思っている。

power.prop.test(p1 = 0.02, 
                p2 = 0.04,  
                sig.level = 0.05, 
                power = 0.8,
                alternative = "two.sided")
pwr.2p.test(ES.h(0.01, 0.02), power = 0.8, sig.level = 0.05, alternative = "two.sided")

pwr.t.test(d=0.2, power = 0.8, sig.level= 0.01)

#サンプルサイズ設計問題〜誤差の観点から〜
#delta=抑えたい誤差の範囲
#alpha=有意水準
#p=母比率。わからない時は0.5

size_calc <-function(delta,p,alpha){
  ceiling(((qnorm(alpha/2) ^ 2) *p*(1-p))/delta^2)
}
size_calc(delta = 0.01, p = 0.5, alpha = 0.05)