r語言 做wald檢驗_r語言wald檢驗怎麼做 - CSDN

2020-11-21 CSDN技術社區

用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/")

 

 

 

 

相關焦點

  • R語言 小wald檢驗_lm檢驗 wald檢驗 - CSDN
    用R語言遇到的一些問題。經常看到rcs()函數,比如擬合回歸時:f <- cph(S ~ rcs(age,4) + sex, x=T, y=T)。
  • 似然比檢驗、Wald檢驗和拉格朗日檢驗的Stata實現 討論
    似然比檢驗(LR)、Wald檢驗、拉格朗日檢驗(LM)都基於最大似然估計(MLE),本文以logit模型為例討論三類檢驗的
  • r語言卡方檢驗和似然比檢驗_r語言似然比檢驗代碼 - CSDN
    本文對應《R語言實戰》第9章:方差分析;第10章:功效分析 ====================================================================方差分析:回歸分析是通過量化的預測變量來預測量化的響應變量,而解釋變量裡含有名義型或有序型因子變量時
  • r語言的p值檢驗 - CSDN
    輸入1: rdata = matrix(rnorm(1000* 6, 0, 3), 6) rvar = apply(rdata, 2, var) mean(rvar)結果1: 前文連結:醫學統計與R語言:多列分組正態性檢驗醫學統計與R語言:標準Z值一定服從標準正態分布?
  • adf檢驗r語言分析_r語言adf檢驗 - CSDN
    協整檢驗是為了檢驗非平穩序列的因果關係,協整檢驗是解決偽回歸為問題的重要方法。利用最小二乘法對回歸方程進行估計,從回歸方程中提取殘差進行檢驗。,適用於對多重一階單整I(1)序列進行協整檢驗。JJ檢驗有兩種:特徵值軌跡檢驗和最大特徵值檢驗。我們可以調用urca包中的ca.jo命令完成這兩種檢驗。
  • r語言卡方檢驗算法_r語言符號檢驗算法 - CSDN
    解:按題意,需檢驗                     H0: μ ≤ 225     H1: μ >  225      此問題屬於單邊檢驗問題      可以使用R語言t.test
  • r語言做白噪聲檢驗_r語言中如何做白噪聲檢驗 - CSDN
    我個人認為時間序列分析是一門挺重要的科目,如果做建模什麼的一定是知道的,或者處理數據的時候,很多數據都是和時間有關的,所以時間序列還是很值得學習的。    這次我申請了一個專欄,我會把文章放在專欄裡。截一張圖,做一個紀念。
  • r語言白噪聲檢驗眼_r語言白噪聲檢驗 - CSDN
    我個人認為時間序列分析是一門挺重要的科目,如果做建模什麼的一定是知道的,或者處理數據的時候,很多數據都是和時間有關的,所以時間序列還是很值得學習的。    這次我申請了一個專欄,我會把文章放在專欄裡。截一張圖,做一個紀念。
  • r 平穩性檢驗 語言_r語言平穩性檢驗方法 - CSDN
    協整檢驗是為了檢驗非平穩序列的因果關係,協整檢驗是解決偽回歸為問題的重要方法。利用最小二乘法對回歸方程進行估計,從回歸方程中提取殘差進行檢驗。,適用於對多重一階單整I(1)序列進行協整檢驗。JJ檢驗有兩種:特徵值軌跡檢驗和最大特徵值檢驗。我們可以調用urca包中的ca.jo命令完成這兩種檢驗。
  • r 秩和檢驗 - CSDN
    所述配對雙樣品的Wilcoxon檢驗一種的非參數檢驗,其可以被用於比較樣品的兩個獨立數據。 本文介紹如何在ř中計算兩個樣本的秩檢驗。
  • r語言檢驗 是否相關 - CSDN
    #分析:按題意,需檢驗#H0: μ ≤ 225 H1: μ > 225#此問題屬於單邊檢驗問題,可以使用R語言t.test#t.test(x,y=NULL,alternative=c(「two.sided」,「less」,「greater」),mu=0,paired=FALSE,var.equal=FALSE,conf.level=0.95
  • r語言 t檢驗 假設 - CSDN
    假設檢驗 -T檢驗 -F檢驗 -卡方檢驗 -正太性檢驗T檢驗2兩樣本的T檢驗 -有原始數據的獨立兩樣本T檢測 -有原始數據的配對T檢測 實例如下: Wage 數據中大學學歷的收入和中學一樣嗎
  • R語言:t檢驗
    BioNews,專注於報導生命科學領域相關新聞,長按下方二維碼即可關注"BioNews"(id : iBioNews)1. t分布(不同自由度)了解r語言幾個函數單樣本t檢驗(使用教材光碟血紅蛋白數據: 例03-05.sav)前提條件:取自正態分布的小樣本(<=60, 偏態用秩和檢驗);或者取自任意分布的大樣本(>60)# install.packages("memisc")library(memisc)## Loading required package: lattice
  • r語言 平穩性檢驗 - CSDN
    這一部分是時間序列預處理R語言的實現。目標是將課本和上課知識點整合。老師是用一節課講完的,本篇文章只做了平穩性檢驗~~~下一篇再寫純隨機性檢驗全部代碼yield <- c(15.2,16.9,15.3,14.9,15.7,15.1,16.7)par(mfrow=c(1,2))plot(yield)yield <- ts(yield,start = 1884)plot(yield)help(plot)plot(yield,
  • r語言重格蘭傑因果檢驗代碼_r語言格蘭傑因果檢驗代碼 - CSDN
    本博客在另外一篇文章中(參考文末給出的資料連結【1】)介紹過在R語言中進行格蘭傑因果檢驗的方法,本文將在解釋格蘭傑因果檢驗原理的基礎上,演示在Python中進行格蘭傑因果檢驗的具體步驟。二、Granger因果關係檢驗 變量X是否為變量Y的Granger原因,是可以檢驗的。 檢驗X是否為引起Y變化的Granger原因的過程如下:第一步:檢驗原假設「H0:X不是引起Y變化的Granger原因」。
  • 【R語言】相關性分析、相關係數的顯著性檢驗及可視化
    本篇文章介紹基於R語言的相關性分析、相關係數的顯著性檢驗及可視化,該教程為個人筆記,大家也可參考學習,不足之處也歡迎大家批評指正!相關性分析用於評估兩個或多個變量之間的關聯,能通過定量指標描述變量之間的強弱、直接或間接聯繫。
  • r語言檢驗序列相關 - CSDN
    ,則進行模式識別畫自相關圖和非自相關圖,根據兩圖的結尾性和拖尾性進行AR、MA、ARMA的模式識別對識別後模式中的位置參數進行參數估計arima()模型檢驗分為:①殘差的白噪聲檢驗;②過度擬合檢驗pt()模型檢驗通過則進行模型優化,否則重新進行模式識別模型優化中得到AIC和BIC值,進行模型的優化然後進行預測與控制2.
  • dw檢驗_dw檢驗的r語言代碼 - CSDN
    線性回歸—投資額(python、OLS最小二乘、殘差圖、DW檢驗)一、問題描述:    建立投資額模型,研究某地區實際投資額與國民生產值(GNP)及物價指數(PI)的關係,根據對未來GNP及PI的估計,預測未來投資額
  • r語言 似然比檢驗_對數似然比檢驗的r語言實現 - CSDN
    學習目標使用LRT提取結果,並與Wald檢驗進行比較從LRT顯著基因列表中識別共享表達譜似然比檢驗(LRT)結果探索DESeq2還提供了似然比檢驗作為跨兩個以上組別評估表達變化
  • r語言tseries - CSDN
    tsdiag(m1) #對估計進行診斷,判斷殘差是否為白噪聲summary(m1)r=m1$residuals #用r來保存殘差Box.test(r,type=」Ljung-Box」,lag=6, fitdf=1)#對殘差進行純隨機性檢驗,fitdf表示殘差減少的自由度AutocorTest(m1$resid) #加載FinTS包