R in action -- chapter 7

mycars <- c("mpg",'hp','wt')
head(mtcars[mycars])
summary(mtcars[mycars])

mystats <- function(x,na.omit=FALSE){   #偏度峰度
  if (na.omit)
    x <- x[!is.na(x)]
  m <- mean(x)
  n <- length(x)
  s <- sd(x)
  skew <- sum((x-m)^3/s^3)/n
  kurt <- sum((x-m)^4/s^4)/n - 3
  return(c(n=n,mean=m,stdev=s,skew=skew,kurtosis=kurt))
}
sapply(mtcars[mycars], mystats)
sapply(mtcars[mycars], mystats,na.omit=T)

install.packages('Hmisc')
library(Hmisc)
describe(mtcars[mycars])

library(pastecs)
# stat.desc(x,basic=T,desc=T,norm=F,p=0.95)
stat.desc(mtcars[mycars])

library(psych)
describe(mtcars[mycars])

head(mtcars)
aggregate(mtcars[mycars],by=list(am=mtcars$am),mean) #分组 描述
aggregate(mtcars[mycars],by=list(am=mtcars$am),sd)

dstats <- function(x)sapply(x,mystats)
a <- by(mtcars[mycars],mtcars$am,dstats)


library(doBy)
summaryBy(mpg+hp+wt~am,data=mtcars,FUN=mystats)


library(psych)
describeBy(mtcars[mycars],list(am=mtcars$am))


library(vcd)
head(Arthritis)

table(Arthritis$Treatment,Arthritis$Improved)
with(Arthritis,table(Treatment,Improved))
mytable <- xtabs(~Treatment+Improved,data = Arthritis)
with(Arthritis,xtabs(~Treatment+Improved,data = Arthritis))

margin.table(mytable,2) # sum by row
prop.table(mytable,2)  #proportion by column
prop.table(mytable)  #proportion by total

addmargins(mytable)
addmargins(mytable,1)
addmargins(prop.table(mytable,2),1)


library(gmodels)
a <- CrossTable(Arthritis$Treatment,Arthritis$Improved) ##SAS format

mytable <- xtabs(~Treatment+Improved,data = Arthritis)
chisq.test(mytable)

mytable <- xtabs(~Sex+Improved,data = Arthritis)
chisq.test(mytable)
fisher.test(mytable)

mytable <- xtabs(~Treatment+Improved+Sex,data = Arthritis)
mantelhaen.test(mytable)

mytable <- xtabs(~Sex+Improved,data = Arthritis)
assocstats(mytable)

head(state.x77)
states <- state.x77[,1:6]
cov(states)
cor(states)
cor(states,method = 'spearman')

x <- states[,c('Population','Income','Illiteracy','HS Grad')]
y <- states[,c('Life Exp','Murder')]
cor(x,y)

#partial correlation
library(ggm)
colnames(states)
pcor(c(1,5,2,3,6),cov(states))  # convarice between first variables(1 and 5)

cor.test(states[,3],states[,5])

library(psych)
corr.test(states,use='complete')

library(MASS)
head(UScrime)
t.test(Prob~So, data = UScrime)

sapply(UScrime[c('U1','U2')], function(x)c(mean=mean(x), sd=sd(x)))

with(UScrime,t.test(U1,U2,paired = T))

#nonparameter
with(UScrime,by(Prob,So,mean))
wilcox.test(Prob~So,data = UScrime)

sapply(UScrime[c('U1','U2')], median)
with(UScrime,wilcox.test(U1,U2,paired = T))

states <- data.frame(state.region,state.x77)
kruskal.test(Illiteracy~state.region,data=states)

source("http://www.statemethods.net/RiA/wmx.txt")
wmc(Illiteracy~state.region,state.x77)