使用的是R包glmnet 函数,有些时候由于方法不收敛可以增加迭代次数,使得模型收敛。
log_res<-glm(pheno~PRS,family=binomial(link='logit'),control=list(maxit=150),data=train_data)
summary(log_res)
logistic.display(log_res)
交叉验证:1000次10倍交叉验证
install.packages('caret')
install.packages('lattice')
library('caret')
set.seed(7)
require("RColorBrewer")
require(caret)
colss=brewer.pal(8,"Set1")[1:10]
display.brewer.all(type = "qual")
for(j in 1:1000){
pdf(paste0(j,"_perm.pdf",sep=""),width=6,height=6)
a<-dev.cur()
png(paste0(j,"_perm.png",sep=""))
dev.control("enable")
folds<-createFolds(y=train_data$pheno,k=5)
for (i in 1:5)
{
print(i)
fold_test<-train_data[folds[[i]],]
fold_train<-train_data[-folds[[i]],]
print ("进行中")
##fold_pre<-glm(pheno~age+gender+PRS,family=binomial(link='logit'),data=fold_train)
##fold_pre<-glm(pheno~PRS,family=binomial(link='logit'),data=fold_train)
summary(fold_pre)
fold_predict <- predict.glm(fold_pre,type='response',newdata=fold_test)
fold_predicts = ifelse(fold_predict > 0.5,1,0)
fold_test$predict = fold_predicts
fold_real<-fold_test$pheno
modelroc<-roc(fold_real,fold_predict)
print(modelroc$auc)
if(i==1)
{
#plot.roc(modelroc, percent = TRUE, main="Smoothing",print.auc =TRUE)
plot(modelroc,col=colss[i],
auc.polygon = TRUE,max.auc.polygon=TRUE,
max.auc.polygon.col="white",auc.polygon.col = "white",
legacy.axes=TRUE,lwd=2,identity=F,main="ROC curve")
}
else
{
plot(modelroc, add=TRUE,col=colss[i],
max.auc.polygon.col="white",
auc.polygon.col = "white",legacy.axes=TRUE,lwd=2,identity=F)
}
legend(0.5,0.1+0.05*i,paste("AUC=",signif(modelroc$auc,3)),col=colss[i],lty=2,cex=1.5,bty="n")
}
dev.copy(which=a)
dev.off()
dev.off()
}
```
网友评论