#GLM

#logistic

install.packages("AER")

library(AER)

summary(Affairs)

table(Affairs$affairs)

fit.full <- glm(ynaffair ~ gender + age + yearsmarried + children +religiousness + education + occupation+rating,data=Affairs, family=binomial())

summary(fit.full)

fit.reduced <- glm(ynaffair ~ age + yearsmarried +religiousness +rating,data=Affairs, family=binomial())

summary(fit.reduced)

anova(fit.reduced, fit.full, test="Chisq")

testdata <- data.frame(rating=c(1, 2, 3, 4, 5), age=mean(Affairs$age),yearsmarried=mean(Affairs$yearsmarried),religiousness=mean(Affairs$religiousness))

testdata$prob <- predict(fit.reduced, newdata=testdata, type="response")

testdata <- data.frame(rating=mean(Affairs$rating),

                      age=seq(17, 57, 10),

                      yearsmarried=mean(Affairs$yearsmarried),

                      religiousness=mean(Affairs$religiousness))

testdata$prob <- predict(fit.reduced, newdata=testdata,  type="response")

#对过度离势进行检验

deviance(fit.reduced)/df.residual(fit.reduced)

#置换检验 coin

library(coin)

score <- c(40, 57, 45, 55, 58, 57, 64, 55, 62, 65)

treatment <- factor(c(rep("A",5), rep("B",5)))

mydata <- data.frame(treatment, score)

t.test(score~treatment, data=mydata, var.equal=TRUE)

oneway_test(score~treatment, data=mydata, distribution="exact")

library(MASS)

UScrime <- transform(UScrime, So = factor(So))

wilcox_test(Prob ~ So, data=UScrime, distribution="exact")

library(multcomp)

set.seed(1234)

oneway_test(response~trt, data=cholesterol,  

            distribution=approximate(nresample=9999))

#卡方检验

library(vcd)

Arthritis <- transform(Arthritis,Improved=as.factor(as.numeric(Improved)))

chisq_test(Treatment~Improved,data = Arthritis,distribution = approximate(B=9999))

#线性模型置换检验lmPerm

require(lmPerm)

fit<-lmp(weight~height,data=women,perm="Prob")

summary(fit)

psych包：principal() 功能较强

#主成分

install.packages("psych")

library(psych)

pc<-principal(USJudgeRatings[,-1],nfactors=1)

pc<-principal(USJudgeRatings[,-1],nfactors=11, rotate = "varimax")

score <- pc$scores#获取主成分得分

#层次聚类分析

row.names(nutrient)<-tolower(row.names(nutrient))

nutrient.scaled<-scale(nutrient)

d<-dist(nutrient.scaled)

fit.average<-hclust(d,method = "average")

plot(fit.average,hang=-1,cex=.8,main="Average Linkage Clustering")

clusters <- cutree(fit.average, k=5)

aggregate(as.data.frame(nutrient.scaled),by=list(cluster=clusters),median)

plot(fit.average, hang=-1, cex=.8,

    main="Average Linkage Clustering\n5 Cluster Solution")

rect.hclust(fit.average, k=5)

#分类算法

#数据准备

loc <- "http://archive.ics.uci.edu/ml/machine-learning-databases/"

ds <- "breast-cancer-wisconsin/breast-cancer-wisconsin.data"

url <- paste(loc, ds, sep="")

breast <- read.table(url, sep=",", header=FALSE, na.strings="?")

names(breast) <- c("ID", "clumpThickness", "sizeUniformity",

                  "shapeUniformity", "maginalAdhesion",

                  "singleEpithelialCellSize", "bareNuclei",

                  "blandChromatin", "normalNucleoli", "mitosis", "class")

df <- breast[-1]

df$class <- factor(df$class, levels=c(2,4),

                  labels=c("benign", "malignant"))

set.seed(1234)

#分为训练集和验证集

train <- sample(nrow(df), 0.7*nrow(df))

df.train <- df[train,]

df.validate <- df[-train,]

table(df.train$class)

table(df.validate$class)

#比较方法 逻辑回归

fit.logit <- glm(class~., data=df.train, family=binomial())

prob <- predict(fit.logit, df.validate, type="response")

logit.pred <- factor(prob > .5, levels=c(FALSE, TRUE),

                    labels=c("benign", "malignant"))

logit.perf <- table(df.validate$class, logit.pred,

                    dnn=c("Actual", "Predicted"))

logit.perf

#计算AUC

install.packages("pROC")

library(pROC)

require(pROC)

auc = roc(df.validate$class,predict(fit.logit, df.validate))

plot(auc)

#比较方法 决策树

library(rpart)

set.seed(1234)

dtree <- rpart(class ~ ., data=df.train, method="class", parms=list(split="information"))

dtree$cptable

dtree.pruned <- prune(dtree, cp=.0125)

dtree.pred <- predict(dtree.pruned,df.validate,type = "class")

install.packages("rpart.plot")

library(rpart.plot)

prp(dtree.pruned, type = 2, extra = 104,

    fallen.leaves = TRUE, main="Decision Tree")

#随机森林

library(randomForest)

set.seed(1234)

fit.forest <- randomForest(class~., data=df.train,

                          na.action=na.roughfix,

                          importance=TRUE)

fit.forest

Import <- importance(fit.forest, type=2)

dotchart(import[order(import),] ,main = "Importance of RF",pch=16)

forest.pred <- predict(fit.forest, df.validate)

forest.perf <- table(df.validate$class, forest.pred,

                    dnn=c("Actual", "Predicted"))

forest.perf

#预测准确性度量

performance <- function(table, n=2){

  if(!all(dim(table) == c(2,2)))

    stop("Must be a 2 x 2 table")

  tn = table[1,1]

  fp = table[1,2]

  fn = table[2,1]

  tp = table[2,2]

  sensitivity = tp/(tp+fn)

  specificity = tn/(tn+fp)

  ppp = tp/(tp+fp)

  npp = tn/(tn+fn)

  hitrate = (tp+tn)/(tp+tn+fp+fn)

  result <- paste("Sensitivity = ", round(sensitivity, n) ,

                  "\nSpecificity = ", round(specificity, n),

                  "\nPositive Predictive Value = ", round(ppp, n),

                  "\nNegative Predictive Value = ", round(npp, n),

                  "\nAccuracy = ", round(hitrate, n), "\n", sep="")

  cat(result)

}

#几种方法的比较

performance(logit.perf)

performance(dtree.perf)

performance(forest.perf)