4-2. 통계분석 with Customized R
# Independent Sample t-test
myttesti <- function(df,nt) {
#**************************************************
josa <- function(word) {
k <- substr(word,nchar(word),nchar(word))
if (((k >="가") & (k <= "힝")) & (((utf8ToInt(k) - utf8ToInt("가")) %% 28) > 0)) {
return (paste0(word,"은"))
} else {
return (paste0(word,"는"))
}
}
#**************************************************
df <- na.omit(data.frame(df))
vname <- colnames(df)
y <- as.matrix(as.numeric(df[,1]))
x <- factor(as.matrix(df[,2]))
cname <- levels(x)
m=aggregate(y ~ x, df, FUN="mean")
s=aggregate(y ~ x, df, FUN="sd")
n=aggregate(!is.na(y) ~ x, df, FUN="sum")
ndf=nrow(df) - 2
sp2 = ((n[1,2]-1)*s[1,2]^2 + (n[2,2]-1)*s[2,2]^2) / ndf
t = (m[1,2]-m[2,2]) / sqrt (sp2/ n[1,2] + sp2/ n[2,2])
p <- (1 - pt(abs(t), ndf))*2
if (nt == 1) {p = p/2.0}
out <- matrix(NaN,3,5)
out[1,1] <- m[1,2]
out[2,1] <- m[2,2]
out[3,1] <- mean(y)
out[1,2] <- s[1,2]
out[2,2] <- s[2,2]
out[3,2] <- sd(y)
out[1,3] <- n[1,2]
out[2,3] <- n[2,2]
out[3,3] <- nrow(df)
out[1,4] <- t
out[1,5] <- p
rownames(out) <- c(cname,"Total")
colnames(out) <- c("Mean","SD","N","t","p")
hypo1 <- sprintf("%s의 %s %s의 %s보다 높을 것이다.",cname[1],josa(vname[2]),cname[2],vname[2])
hypo2 <- sprintf("%s %s에 따라 차이가 있을 것이다.",josa(vname[1]),vname[2])
res0 <- sprintf("%s %s에 따라 통계적으로 차이가 없다",josa(vname[2]),vname[1])
res21 <- sprintf("%s %s에 따라 통계적으로 p < 0.05 수준에서 의미있는 차이가 있다.",josa(vname[1]),vname[2])
res22 <- sprintf("%s %s에 따라 통계적으로 p < 0.01 수준에서 의미있는 차이가 있다.",josa(vname[1]),vname[2])
res23 <- sprintf("%s %s에 따라 통계적으로 p < 0.001 수준에서 의미있는 차이가 있다.",josa(vname[1]),vname[2])
res11 <- sprintf("%s의 %s %s의 %s보다 통계적으로 p < 0.05 수준에서 의미있게 높았다.",cname[1],josa(vname[1]),cname[2],vname[1])
res12 <- sprintf("%s의 %s %s의 %s보다 통계적으로 p < 0.01 수준에서 의미있게 높았다.",cname[1],josa(vname[1]),cname[2],vname[1])
res13 <- sprintf("%s의 %s %s의 %s보다 통계적으로 p < 0.001 수준에서 의미있게 높았다.",cname[1],josa(vname[1]),cname[2],vname[1])
if (nt == 1) {
hypo=hypo1
tail <- "One-Tailed test"
} else {
hypo=hypo2
tail <- "Two-Tailed test"
}
if (p >= 0.05) {
res=res0
} else if (nt == 1 & (p < 0.05 & p >= 0.01)) {
res=res11
} else if (nt == 1 & (p < 0.01 & p >= 0.001)) {
res=res12
} else if (nt == 1 & p < 0.001) {
res=res13
} else if (nt == 2 & (p < 0.05 & p >= 0.01)) {
res=res21
} else if (nt == 2 & (p < 0.01 & p >= 0.001)) {
res=res22
} else if (nt == 2 & p < 0.001) {
res=res23
}
res <- paste0(res, " (t = ",round(t,3),", df = ",ndf,", p = ",round(p,5),")")
result <- list(tail=tail, hypothesis=hypo, out=out, conclusion=res)
return(result)
}
library(openxlsx)
data <- read.xlsx("ttest.xlsx")
# two-tailed test
myttesti(data,2)
## $tail
## [1] "Two-Tailed test"
##
## $hypothesis
## [1] "근무평정점수는 성별에 따라 차이가 있을 것이다."
##
## $out
## Mean SD N t p
## 1 80.28571 9.411239 7 3.880231 0.002562346
## 2 58.83333 10.534072 6 NaN NaN
## Total 70.38462 14.643192 13 NaN NaN
##
## $conclusion
## [1] "근무평정점수는 성별에 따라 통계적으로 p < 0.01 수준에서 의미있는 차이가 있다. (t = 3.88, df = 11, p = 0.00256)"
# 일원분산분석 (One-way ANOVA)
myanova <- function(df) {
#**************************************************
josa <- function(word) {
k <- substr(word,nchar(word),nchar(word))
if (((k >="가") & (k <= "힝")) & (((utf8ToInt(k) - utf8ToInt("가")) %% 28) > 0)) {
return (paste0(word,"은"))
} else {
return (paste0(word,"는"))
}
}
#**************************************************
df <- na.omit(df)
y <- as.matrix(df[,1])
x <- factor(as.matrix(df[,2]))
vname <- colnames(df)
m=aggregate(y ~ x, df, FUN="mean")
s=aggregate(y ~ x, df, FUN="sd")
n=aggregate(!is.na(y) ~ x, df, FUN="sum")
tmean <- mean(y)
tsd <- sd(y)
tn <- nrow(df)
SST <- sum((y-tmean)^2)
SSB <- sum((m[,2] - tmean)^2 * n[,2])
mm=cbind(m[,2],s[,2],n[,2])
out <- matrix(0,3,5)
out[1,1] <- SSB
out[2,1] <- SST-SSB
out[3,1] <- SST
out[1,2] <- nrow(mm)-1
out[2,2] <- sum(n[,2]) - nrow(mm)
out[3,2] <- out[1,2] + out[2,2]
out[1,3] <- out[1,1] / out[1,2]
out[2,3] <- out[2,1] / out[2,2]
out[1,4] <- out[1,3] / out[2,3]
out[1,5] = p = pf(out[1,4],out[1,2],out[2,2],lower.tail = FALSE)
colnames(out) <- c("SS","df","MS","F","p")
rownames(out) <- c("Between","Within","Total")
mm <- rbind(mm,c(tmean,tsd,tn))
stat <- mm
stat <- cbind(mm,matrix(0,nrow(mm),2))
stat[1,4] <- out[1,4]
stat[1,5] <- out[1,5]
rownames(stat) <- c(levels(factor(x)),"Total")
colnames(stat) <- c("Mean","SD","N","F","p")
hypo <- sprintf("%s %s에 따라 차이가 있을 것이다.",josa(vname[1]),vname[2])
res0 <- sprintf("%s %s에 따라 통계적으로 차이가 없다",josa(vname[1]),vname[2])
res1 <- sprintf("%s %s에 따라 통계적으로 p < 0.05 수준에서 의미있는 차이가 있다.",josa(vname[1]),vname[2])
res2 <- sprintf("%s %s에 따라 통계적으로 p < 0.01 수준에서 의미있는 차이가 있다.",josa(vname[1]),vname[2])
res3 <- sprintf("%s %s에 따라 통계적으로 p < 0.001 수준에서 의미있는 차이가 있다.",josa(vname[1]),vname[2])
if (p >= 0.05) {
res=res0
} else if (p < 0.05 & p >= 0.01) {
res=res1
} else if (p < 0.01 & p >= 0.001) {
res=res2
} else if (p < 0.001) {
res=res3
}
res <- paste0(res, " (F = ",round(out[1,4],3),", df1 = ",out[1,2],", df2 = ",out[2,2],", p = ",round(out[1,5],5),")")
result <- list(hypothesis=hypo, out=out, stat=stat, conclusion=res)
return(result)
}
library(openxlsx)
df=read.xlsx("anova1.xlsx")
# One-way ANOVA
myanova(df)
## $hypothesis
## [1] "근무만족도는 상사의.유형에 따라 차이가 있을 것이다."
##
## $out
## SS df MS F p
## Between 40.44444 2 20.222222 4.0625 0.03891091
## Within 74.66667 15 4.977778 0.0000 0.00000000
## Total 115.11111 17 0.000000 0.0000 0.00000000
##
## $stat
## Mean SD N F p
## 1 8.000000 1.095445 6 4.0625 0.03891091
## 2 8.333333 1.032796 6 0.0000 0.00000000
## 3 11.333333 3.559026 6 0.0000 0.00000000
## Total 9.222222 2.602161 18 0.0000 0.00000000
##
## $conclusion
## [1] "근무만족도는 상사의.유형에 따라 통계적으로 p < 0.05 수준에서 의미있는 차이가 있다. (F = 4.063, df1 = 2, df2 = 15, p = 0.03891)"
# 회귀분석 (regression analysis)
library(openxlsx)
myreg <- function(df) {
#*******************************************************************
#*** Required Functions for Regression Analysis: vif,ones,colSds ***
#*******************************************************************
vif <- function(xx) {
vif <- matrix(NaN,ncol(xx)+1,1)
for (i in (1:ncol(xx))) {
zy <- as.matrix(xx[,i])
zx <- as.matrix(xx[,-i])
beta <- solve(t(zx) %*% zx) %*% (t(zx) %*% zy)
ssq2 <- (t(zx %*% beta) %*% zy) / (t(zy) %*% zy)
vif[i+1,1] <- 1 / (1 - ssq2)
}
return(vif)
}
#*******************************************************************
ones <- function(r,c) {
x <- matrix(1,r,c)
return(x)
}
#*******************************************************************
colSds <- function(x){
return(apply(x,2,sd))
}
#*******************************************************************
df <- as.matrix(na.omit(df))
y <- df[,1]
x <- df[,-1]
zdf <- as.matrix((df-ones(nrow(df),1)%*%colMeans(df))/ones(nrow(df),1)%*%colSds(df))
zy <- zdf[,1]
zx <- zdf[,-1]
x <- cbind(ones(nrow(df),1),x)
b <- solve(t(x) %*% x) %*% (t(x) %*% y)
cname <- colnames(x)
meanx <- colMeans(x)
beta <- solve(t(zx)%*%zx) %*% t(zx)%*%zy
sigma <- (t(y-x%*%b) %*% (y-x%*%b)) / (nrow(df) - ncol(df))
ss <- solve(t(x) %*% x) * as.vector(sigma)
se <- matrix(sqrt(diag(ss)))
t <- b/se
dff <- nrow(df) - ncol(df)
p <- matrix(0,ncol(df))
for (i in c(1:ncol(df))) {
p[i] <- (1 - pt(abs(t[i]), dff))*2
}
regtable <- cbind(b,se,rbind(NaN,beta),t,p,vif(zx))
sst <- sum((y - ones(nrow(df),1)*mean(y))^2)
sse <- sum((y-x%*%b)^2)
out <- matrix(NaN,3,5)
out[1,1] <- sst - sse
out[2,1] <- sse
out[3,1] <- sst
out[1,2] <- ncol(df) -1
out[2,2] <- nrow(df) - ncol(df)
out[3,2] <- nrow(df) - 1
out[1,3] <- out[1,1] / out[1,2]
out[2,3] <- out[2,1] / out[2,2]
out[1,4] <- out[1,3] / out[2,3]
out[1,5] <- pf(out[1,4],out[1,2],out[2,2],lower.tail=F)
ssq2 <- (t(zx %*% beta) %*% zy) / (t(zy) %*% zy)
R2 <- ssq2
adjR2 <- 1 - (1-ssq2) * out[3,2]/out[2,2]
colnames(out) <- c("SS","df","MS","F","p")
rownames(out) <- c("Regression","Residual","Total")
cname[1] <- '(Intercept)'
rownames(regtable) <- cname
colnames(regtable) <- c("B","SE","BETA","t","p","VIF")
result <- list(R=round(sqrt(ssq2),3),Rsquared=round(R2,3),
AdjustedRSquared=round(adjR2,3),ANOVATable=round(out,3),
RegressionTable=round(regtable,3))
# plot(df[,2],df[,1])
# abline(b[1],b[2],col="red")
return(result)
}
library(openxlsx)
df <- read.xlsx("regress.xlsx")
# 회귀분석 (regression analysis)
myreg(df)
## $R
## [,1]
## [1,] 0.884
##
## $Rsquared
## [,1]
## [1,] 0.782
##
## $AdjustedRSquared
## [,1]
## [1,] 0.747
##
## $ANOVATable
## SS df MS F p
## Regression 7801.785 4 1950.446 22.415 0
## Residual 2175.415 25 87.017 NaN NaN
## Total 9977.200 29 NaN NaN NaN
##
## $RegressionTable
## B SE BETA t p VIF
## (Intercept) 22.043 13.533 NaN 1.629 0.116 NaN
## 대인관계 -0.008 0.098 -0.009 -0.083 0.935 1.300
## 자아개념 0.528 0.130 0.631 4.066 0.000 2.760
## 근무평정 -0.139 0.152 -0.090 -0.913 0.370 1.126
## SES점수 0.414 0.171 0.339 2.415 0.023 2.253