用R語言遇到的一些問題。
經常看到rcs()函數,比如擬合回歸時:f <- cph(S ~ rcs(age,4) + sex, x=T, y=T)。
關於RCS的理解,可以參考:Restricted cubic splines
另外,丁香園裡有人給出這樣的看法:
rcs全稱是restricted cubic spline 即限制立方樣條函數。為什麼用它呢?我們在做回歸擬合數據時,經常對因變量和自變量的假定是:自變量和因變量呈線性關係,logistic回歸是自變量和logitP,cox是指自變量和一個h(t)變換的一個東東(記不清了,抱歉),這些假定就叫線性假定。我們做回歸的時候自變量包括連續變量和ordinal variable(等級有序變量)和因變量不成線性的時候,我可以轉化為分類,也可以變量log變換,還有就是spline技術,也就是引入多項式,但rcs又不是簡單的引入多項式,簡單的說就是引入3次方的項,這和你選擇的節點有關,這樣處理的結果就是可以和好的擬合數據有限制了引入項的自由度,所以Harrell (rms包的作者)推崇這種建模方法。說的很概括,具體的內容搜索rcs就可以了解。Harrell的論文也可以參考,他認為節點的位置不那麼重要,而節點數目很重要。
rms中的應用:
## Not run:options(knots=4, poly.degree=2)# To get the old behavior of rcspline.eval knot placement (which didnt' handle# clumping at the lowest or highest value of the predictor very well):# options(fractied = 1.0) # see rcspline.eval for detailscountry <- factor(country.codes)blood.pressure <- cbind(sbp=systolic.bp, dbp=diastolic.bp)fit <- lrm(Y ~ sqrt(x1)*rcs(x2) + rcs(x3,c(5,10,15)) + lsp(x4,c(10,20)) + country + blood.pressure + poly(age,2))# sqrt(x1) is an implicit asis variable, but limits of x1, not sqrt(x1)# are used for later plotting and effect estimation# x2 fitted with restricted cubic spline with 4 default knots# x3 fitted with r.c.s. with 3 specified knots# x4 fitted with linear spline with 2 specified knots# country is an implied catg variable# blood.pressure is an implied matrx variable# since poly is not an rms function (pol is), it creates a# matrx type variable with no automatic linearity testing# or plottingf1 <- lrm(y ~ rcs(x1) + rcs(x2) + rcs(x1) %ia% rcs(x2))# %ia% restricts interactions. Here it removes terms nonlinear in# both x1 and x2f2 <- lrm(y ~ rcs(x1) + rcs(x2) + x1 %ia% rcs(x2))# interaction linear in x1f3 <- lrm(y ~ rcs(x1) + rcs(x2) + x1 %ia% x2)# simple product interaction (doubly linear)# Use x1 %ia% x2 instead of x1:x2 because x1 %ia% x2 triggers# anova to pool x1*x2 term into x1 terms to test total effect# of x1## End(Not run)
SVM的調參,關於e1071包,好像如果把數據「尺度化」(scaling)後,使用默認的參數就能訓練出較好的模型。
R語言裡,回歸的參數,如果傳formula,比如Y~X,那麼這裡的X不應該是dataframe或matrix,而應該用向量比如Y~x1+x2。如果向量太多,那麼可以這樣傳兩個參數:formula和data,比如glm(Y~., data=X)。這是在Logistics回歸中被坑過,在predict時解決了將newdata用dataframe形式傳入的問題後,總是報錯說變量數不對(一方面可能fit 時有問「glm傳參數data可能沒傳好」,另一方面newdata的dataframe可能需要轉置一下)。
R語言的向量、矩陣、數據框、數組、列表等等,有點煩人,和我之前的C語言系列思維差別太大。
R的e1701中都在svm中,僅當y變量是factor類型時構建SVC,否則構建SVR
ggplot做出來的圖,如何去除默認的背景和網格?(摘自:移除ggplot2的網格和背景色)
# 方法1myplot<-myplot+ theme(panel.grid.major =element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.line = element_line(colour = "black"))# 方法2myplot<-myplot+theme_bw()+ theme(panel.border = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), axis.line = element_line(colour = "black"))
生存分析
cph與coxph是點不同的,雖然得到的係數一樣,p值一樣,但是wald統計量有些區別(進一步進行anova分析時可以看出區別)。cph來自rms包,coxph來自survival包。
cph的一些進一步的說明:
「Modification of Therneau’s coxph function to fit the Cox model and its extension, the AndersenGill model. The latter allows for interval time-dependent covariables, time-dependent strata, and repeated events. The Survival method for an object created by cph returns an S function for computing estimates of the survival function. The Quantile method for cph returns an S function. for computing quantiles of survival time (median, by default). The Mean method returns a function for computing the mean survival time. This function issues a warning if the last follow-up time is uncensored, unless a restricted mean is explicitly requested.」
繪製Nomogram舉例1:
require(rms) ##調用rms包# 建立數據集y = c(0, 1, 0, 0, 0, 1, 1, 1, 0, 0, 0, 1, 1, 0, 0, 1, 0, 0, 0, 1, 1, 0, 1, 1, 0, 1, 1, 1, 0, 1, 0, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 1, 0, 0, 0, 1, 1, 1, 1, 0, 1, 1, 0, 0, 0, 1, 1, 1, 0, 1, 1, 1, 1, 0, 0, 1, 1, 1, 0, 0, 0, 1, 0, 1, 0, 1, 0, 1)age = c(28, 42, 46, 45, 34, 44, 48, 45, 38, 45, 49, 45, 41, 46, 49, 46, 44, 48, 52, 48, 45, 50, 53, 57, 46, 52, 54, 57, 47, 52, 55, 59, 50, 54, 57, 60, 51, 55, 46, 63, 51, 59, 48, 35, 53, 59, 57, 37, 55, 32, 60, 43, 59, 37, 30, 47, 60, 38, 34, 48, 32, 38, 36, 49, 33, 42, 38, 58, 35, 43, 39, 59, 39, 43, 42, 60, 40, 44)sex = c(0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1)ECG = c(0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 2, 1, 0, 0, 2, 2, 0, 0, 2, 2, 0, 1, 2, 2, 0, 1, 0, 2, 0, 1, 0, 2, 1, 1, 0, 2, 1, 1, 0, 2, 1, 1, 0, 2, 1, 1, 0, 2, 1, 1)ECG<-as.factor(ECG)# 設定nomogram的參數ddist <- datadist(age, sex, ECG)options(datadist='ddist')# logistic回歸f <- lrm(y ~ age + sex + ECG)# nomogramnom <- nomogram(f, fun=plogis, fun.at=c(.001, .01, .05, seq(.1,.9, by=.1), .95, .99, .999), lp=F, funlabel="Risk")plot(nom)
繪製nomagram舉例2
library(foreign) # 讀取spss格式數據library(survival)library(rms) library(VIM) # 包中aggr()函數,判斷數據缺失情況par(family = "STXihei") # 指定繪製圖片中的字體# 數據來自張文彤《SPSS統計分析高級教程》,侵刪!url <- 'http://online.hyconresearch.com:4096/spss20/spss20advdata.zip'# 書中數據下載地址file <- basename(url) # 獲取文件名 spss20advdata.zipif (!file.exists(file)) download.file(url, file) # 判斷工作目錄下是否有spss20advdata.zip文件,如果不存在則執行下載數據命令unzip(file) # 解壓縮zip文件pancer <- read.spss('part4/pancer.sav') # 讀取文件pancer <- as.data.frame(pancer) # Cox回歸所需數據類型為數據框格式,將其轉換為數據框格式aggr(pancer,prop=T,numbers=T) # 判斷pancer各個變量的數據缺失情況,出現紅色即有缺失,繪製列線圖不允許缺失值存在pancer$censor <- ifelse(pancer$censor=='死亡',1,0) # Cox回歸結局變量需為數值變量pancer$Gender <- as.factor(ifelse(pancer$sex=='男',"Male","Female"))# 更改變量名稱以及變量取值pancer$Gendr <- relevel(pancer$Gender,ref='Female') # 設置參考組dd<-datadist(pancer) options(datadist='dd') coxm <- cph(Surv(time,censor==1)~age+Gender+trt+bui+ch+p+stage,x=T,y=T,data=pancer,surv=T) # 建立Cox回歸方程surv <- Survival(coxm) # 建立生存函數surv1 <- function(x)surv(1*3,lp=x) # lp: linear predictor, 用這個函數去surv2 <- function(x)surv(1*6,lp=x) surv3 <- function(x)surv(1*12,lp=x) # maxscale 參數指定最高分數,一般設置為100或者10分# fun.at 設置生存率的刻度# xfrac 設置數值軸與最左邊標籤的距離,可以調節下數值觀察下圖片變化情況plot(nomogram(coxm,fun=list(surv1,surv2,surv3),lp= F, funlabel=c('3-Month Survival','6-Month survival','12-Month survival'),maxscale=100, fun.at=c('0.9','0.85','0.80','0.70','0.6','0.5','0.4','0.3','0.2','0.1')),xfrac=.45)med <- Quantile(coxm)surv4 <- function(x) med(x)plot(nomogram(f, fun=surv4, fun.at=c(12,6,3,1,0),lp=F, funlabel="Median Survival Time"))
繪製Nomogram舉例3:
require(rms)require(Hmisc) ##需要下載安裝require(survival) ##R默認自帶的用於做生存分析的包# 建立數據集(使用rms包example的代碼,未改動)n <- 1000set.seed(731)age <- 50 + 12*rnorm(n)# label(age) <- "Age"sex <- factor(sample(c('Male','Female'), n, rep=TRUE, prob=c(.6, .4)))cens <- 15*runif(n)h <- .02*exp(.04*(age-50)+.8*(sex=='Female'))dt <- -log(runif(n))/h# label(dt) <- 'Follow-up Time'e <- ifelse(dt <= cens,1,0)dt <- pmin(dt, cens)# units(dt) <- "Year"# 設定nomogram的參數ddist <- datadist(age, sex)options(datadist='ddist')# Cox回歸S <- Surv(dt,e)f <- cph(S ~ rcs(age,4) + sex, x=T, y=T, surv=T)med <- Quantile(f)# nomogramnom <- nomogram(f, fun=function(x) med(x), fun.at=c(15,12,11,9,6,5),lp=F, funlabel="Median Survival Time")plot(nom) ##繪製Nomgram圖
先貼代碼,以後整理下。再就是,可以研究下rms包文檔裡面的nomogram例子。
-----------------------------------Fisher exact test------------------------------------
Fisher exact test 不僅適用於2*2表,在R*C表中也能使用Fisher精確檢驗。
有文獻對此進行過說明(Statistical notes for clinical researchers: Chi-squared test and Fisher's exact test)。
-----------------------------------非正態數據的正態化------------------------------------
http://blog.sina.com.cn/s/blog_13ec735f50102x59u.html
https://blog.csdn.net/qq_40587575/article/details/84349900
偏度和峰度的計算:
http://blog.sina.com.cn/s/blog_4d2fda500102wk0r.html
---------------------------------------------------------------------------------------------------
Bioconductor包的安裝——老方法
在BioC的官方網站上,存放著Bioc包的安裝腳本,biocLite.R,每次需要安裝BioC的包之前,我們運行該腳本。source是運行代碼腳本的命令,可以是本地文件,可以是網絡上的文件。需要你有流暢的網絡連結。
source(「http://bioconductor.org/biocLite.R」)
biocLite()是安裝函數,相當於R中的常用的install.package命令,如果不傳遞需要安裝的包的名稱,則安裝一些基本組件。
Bioconductor包的安裝——新方法
if (!requireNamespace("BiocManager", quietly = TRUE)) install.packages("BiocManager")BiocManager::install()
這種方法是將安裝函數打包在BiocManager這個包裡,而這個包又放在了CRAN上,相比老方法更為安全(直接掛網上的代碼容易被黑客做手腳)。
---------------------------------------------------------------------------------------------------
安裝包的方法:
除了install.packages()和BiocManager::install()這兩種經典方法之外,還有如下的方法。
方法1 找到安裝包所在的url,install.packages(packageurl, repos=NULL, type="source")。
方法2 github安裝,需要devtools工具包,install_github()。
方法3 將安裝包下載到本地,然後 install.packages(file, repos = NULL) 。
但是本地安裝一般都容易出錯,因為經常需要安裝依賴包。
---------------------------------------------------------------------------------------------------
library和require都可以載入包,但二者存在區別。
在一個函數中,如果一個包不存在,執行到library將會停止執行,require則會繼續執行。
---------------------------------------------------------------------------------------------------
關於依賴包的更新
#Update old packages: #Update all/some/none? [a/s/n]
這種情況其實可以選擇n的,放心去嘗試吧。
---------------------------------------------------------------------------------------------------
setdiff函數是有方向的,是第一個減去第二個參數。
---------------------------------------------------------------------------------------------------
Rstudio的使用習慣糾正
1、一個正式的項目,請新建一個project,並且階段性收工的時候應當push到雲端。
2、一個project內,應當有如下幾個子目錄:data存放數據、R存放代碼、plots存放圖片、以及各其他子分析模塊存放對應文件。
3、Workspace數據應當保存。
4、每個項目可以有臨時代碼文件,用於嘗試,請將相應代碼保存到temp子目錄中,並以temp.R命名。
5、可以有全局的臨時代碼(與項目無關),可以在電腦中建立一個temp目錄,用於保存temp代碼。
6、每一個project最好至少生成一份rmarkdown報告。
7、一定記得將代碼push到線上倉庫,將數據push到網盤中;需要展示的結果可使用github的page功能。
8、Rmarkdown的高效使用方法:
8.0 雖然Rmarkdown沒有Jupyter notebook那麼直觀好看,但是對於echo=FALSE這個設定的存在很滿意,以及可以輸出office系列文檔,可以做各種slide,很不錯。
8.1 請不要將cache設置為TRUE,緩存機制容易讓自己的代碼運行出錯(在沒能完全了解清楚cache的原理時,不要使用)。
8.2 請將耗時操作運算(一般運行超過1min的都可歸為耗時操作)的代碼標記為eval=FALSE,並且將其運算結果緩存到本地(以表格文件、圖片形式、R格式數據形式保存),eval的代碼加載數據即可;一個合理的Rmarkdown運行及渲染的時間不應當超過1min。
8.3 代碼實行分塊管理,即大的模塊之間使用注釋分割符號進行劃分(大塊間使用等號系列分割、小塊之間使用減號系列分割);無論代碼是否eval,相同功能模塊的代碼都應該寫在同一個塊內,當然,塊內可分為eval區和非eval區。另外,第一個塊內需要設定工作目錄、全局使用的工具包和代碼(如dplyr、ggplot2、FanCodeV1等)。
8.4 選項設置:全局請使用cache=FALSE,message=FALSE;局部代碼根據情況使用eval和warning;局部適當設置fig系列選項。
8.5 代碼可以先在臨時script文件中整理好再導入Rmarkdown,但運行調試也可以適當使用Rmarkdown,畢竟數據查看很方便。
---------------------------------------------------------------------------------------------------
設置全局鏡像:
options(repos="https://mirrors.tuna.tsinghua.edu.cn/CRAN/")