R语言 实例操作2

for循环、assign()函数、paste()函数、sapply()函数、自编简单函数

setwd("D:/FDU/mine/mydata")
options("scipen"=100, "digits"=4)
dat = read.csv( file='D:/FDU/mine/mydata/data.csv', header = T)
dat <- within(dat, {
  WAVE <- as.factor(WAVE)
  GENDER <- as.factor(GENDER)
  edu <- as.factor(edu)
  smoke <- as.factor(smoke)
  drink <- as.factor(drink)
})
###该数据库是针对同一人群七年随访的数据，20个变量。
################将income和index变量分组的函数
fun1 <- function(x){
  q <- quantile(x$income, probs = c(1/3, 2/3), na.rm = T, names = F) ###计算三分位数
  p <- quantile(x$index, probs = c(1/3, 2/3), na.rm = T, names = F)
  x <- within(x, {
    incomel <- NA
    incomel[income < q[1]] <- "low income"
    incomel[income >= q[1] & income < q[2]] <- "middle income"
    incomel[income >= q[2]] <- "high income"
    indexl <- NA
    indexl[index < p[1]] <- "low urbanization"
    indexl[index >= p[1] & index < p[2]] <- "middle urbanization"
    indexl[index >= p[2]] <- "high urbanization"
    incomel <- as.factor(incomel)
    indexl <- as.factor(indexl)
  })
}

##将七年的数据拆成七个数据框，后续还会用到。对income和index分别按照三分位数分组后再合并
for(j in c(1993, 1997, 2000, 2004, 2006, 2009, 2011)){
  assign(paste("y", j, sep = ""), fun1(subset(dat, WAVE == j)))
}

datn <- rbind(y1993, y1997)
datn <- rbind(datn, y2000)
datn <- rbind(datn, y2004)
datn <- rbind(datn, y2006)
datn <- rbind(datn, y2009)
datn <- rbind(datn, y2011)


#############################连续性变量描述性分析
####连续性变量为数据框的第5, 8, 9, 10, 15, 17, 19列,共7个变量。

##对x数据框的7个变量计算中位数
fun2 <- function(x){
  var <- c(5, 8, 9, 10, 15, 17, 19)
  media <- sapply(x[, var[1:7]], median, na.rm = T)
  return(round(media, 2))
}
##对x数据框的7个变量计算四分位数
fun3 <- function(x){
  var <- c(5, 8, 9, 10, 15, 17, 19)
  Q <- sapply(x[, var[1:7]], quantile, probs = c(0.25, 0.75), names = F, na.rm = T)
  return(round(Q, 2))
}

##输出形式########**缺失值NA怎么让它显示为空格而不是NA？？？paste函数的用法**
fun4 <- function(x){
  rslt <- paste(fun2(x), "(", fun3(x)[1, ], ",", fun3(x)[2, ], ")")
  return(rslt)
}
col <- c(fun4(y1993), fun4(y1997), fun4(y2000), fun4(y2004), fun4(y2006), fun4(y2009), fun4(y2011))

##Kruskal-Wallis test及p值输出形式
fun5 <- function(x){
  kru <- kruskal.test(datn[, x] ~ datn$WAVE)
  if (kru[[3]] < 0.0001) print(paste("<", 0.001, sep = ""))
  else print(round(kru[[3]], 3))
}
col2 <- c(fun5(8), fun5(9), fun5(10), fun5(15), fun5(17), fun5(19))

##输出矩阵
m1 <- matrix(nrow = 7, ncol = 8)
m1[, 1:7] <- matrix(col, byrow = F) 
m1[2:7, 8] <- matrix(col2, byrow = F)


##############################分类变量描述性分析，6个变量

##计算七年第i列分类变量的频数和比例，每个变量的结果以矩阵形式输出
for(i in c(4, 7, 13, 14, 21, 22)){
  assign(paste("n", i, sep = ""), table(datn[, i], datn[, 3]))
  assign(paste("prop", i, sep = ""), prop.table(table(datn[, i], datn[, 3]), 2))
}

##对6个分类变量进行卡方检验， 结果为列表
list <- list(n4, n7, n13, n14, n21, n22)
chi <- sapply(list, chisq.test)

##输出形式
fun6 <- function(x, y){
  rslt <- paste(x, "(", round(y, digits = 2), ")")
  return(rslt)
}
fun7 <- function(x){
  if (chi[[3, x]] < 0.0001) print(paste("<", 0.001, sep = ""))
  else print(round(chi[[3, x]], 3))
}

##矩阵输出
m2 <- matrix(nrow = 21, ncol = 8)
m2[1, 1:7] <- c(5457, 7577, 7715, 7603, 7494, 7380, 6456)
m2[2:3, 1:7] <- matrix(fun6(n4, prop4), byrow = F) 
m2[5:7, 1:7] <- matrix(fun6(n7, prop7), byrow = F)
m2[9:10, 1:7] <- matrix(fun6(n13, prop13), byrow = F)
m2[12:13, 1:7] <- matrix(fun6(n14, prop14), byrow = F)
m2[15:17, 1:7] <- matrix(fun6(n21, prop21), byrow = F)
m2[19:21, 1:7] <- matrix(fun6(n22, prop22), byrow = F)
m2[4, 8] <- fun7(2)
m2[8, 8] <- fun7(3)
m2[11, 8] <- fun7(4)
m2[14, 8] <- fun7(5)
m2[18, 8] <- fun7(6)
m <- rbind(m2, m1)