R语言笔记3_回归分析（EDA OLS Power）

R语言笔记3——回归分析（EDA OLS Power）

探索性数据分析EDA

# Exploratory data analysis
#install.packages("faraway")
library(faraway)
head(pima)
summary(pima)

## sanity checks(合理性检查)
sort(pima$diastolic) #sort()函数是对向量进行从小到大的排序
pima$diastolic[pima$diastolic==0] <- NA
pima$glucose[pima$glucose==0] <- NA
pima$triceps[pima$triceps==0] <- NA
pima$insulin[pima$insulin==0] <- NA
pima$bmi[pima$bmi==0] <- NA

## factors - categorical variables
head(pima$test)
pima$test <- factor(pima$test)
summary(pima$test)
levels(pima$test) <- c('negative','positive')
summary(pima$test)

## Histogram
hist(pima$diastolic)

## Kernel density estimate
plot(density(pima$diastolic, na.rm = TRUE)) #删除缺失值

## Scatterplot
plot(diabetes~diastolic,pima)

## side-by-side boxplot
boxplot(diabetes~test,pima)

用scatterplot3d画三维散点图

#install.packages("scatterplot3d")
library(scatterplot3d)
data(iris)
head(iris)
length(unique(iris$Species))

# Basic 3d graphics
scatterplot3d(iris[,1:3])
colors <- c("#999999", "#E69F00", "#56B4E9")
colors <- colors[as.numeric(iris$Species)]
scatterplot3d(iris[,1:3], pch = 16, color = colors)#pch表示点的样式，pch:15~20点为实心
scatterplot3d(iris[,1:3], pch = 16, type = "h", color = colors)

最小二乘法OLS

响应变量：the lean（Y）

解释变量：year（X）

# 最小二乘法的例子
## data input
towerpisa <- list()
towerpisa$year <- c(75:87,116)
towerpisa$lean <- c(642,644,656,667,673,688,696,698,713,717,725,742,757)
print(towerpisa)

## scatter plot
attach(towerpisa)#为直接引用数据集中的元素
plot(year[1:13], lean[1:13], xlab = "Year", ylab = "Lean", las=1)#las表示刻度数字与坐标轴的方向
lines(lowess(year,lean), col=2)#lowess 是一种局部加权回归散点图平滑化技术

## 线性回归分析
bc <- lm(lean ~ year) # R语言线性回归函数 lm()
print(summary(bc))
e = bc$residuals
sse = 0
for( i in 1:length(e)){
  sse = sse + e[i]^2
} #SSE
mse = sse/11 #df = 13-2 =11
sqrt(mse) # residual standard error这么计算的结果与回归lm得出的一样

## residual plot
plot(year[1:13], bc$residuals, las=1, ylab = "Residual", xlab = "Year")
abline(h=0)

## QQ plot #Q-Q 图通过将两个概率分布的相同分位数点的值映射为 x 和 y 轴。
qqnorm(bc$residuals, datax = T)
qqline(bc$residuals, datax = T)

功效计算

# Power for regression

## only these values need to be changed
n = 25
sig2 = 2500
ssx = 19800
beta1 = 1.5
alpha = 0.05

## function of calculating power
find_power <- function(n,sig2,ssx,beta1,alpha){
  sig2b1 = sig2/ssx
  df = n-2
  delta = beta1/sqrt(sig2b1)
  c = qt(1 - alpha/2, df) #非中心t分布的分位数函数
  power = 1 - pt(c,df,delta) + pt(-c,df,delta) #pt(c,df,delta)=P(T<c),T~t(df,delta)
  #print(power)
  return(power)
}

##example
power = find_power(n,sig2,ssx,beta1,alpha)
beta1 = seq(-2,2,by=0.05)
power_vec = rep(0,length(beta1))
for (i in 1:length(beta1)) {
  power_vec[i] = find_power(n,sig2,ssx,beta1[i],alpha)  
}

plot(power_vec ~ beta1, type="l",ylab="Power",main="Power for the slope in simple linear regression")