p.24～
＜Two-sided, unpaired t-test＞
t検定＝データ分布の平均の有意差を検定する
データが正規分布の時に使う
> data(sleep)
> sleep
extra group ID
1 0.7 1 1
2 -1.6 1 2
3 -0.2 1 3
> dim(sleep)
[1] 20 3
> sleep[,1]
[1] 0.7 -1.6 -0.2 -1.2 -0.1 3.4 3.7 0.8 0.0 2.0 1.9 0.8 1.1
[14] 0.1 -0.1 4.4 5.5 1.6 4.6 3.4
> sleep[,2]
[1] 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2
Levels: 1 2
> test <- t.test(sleep[,1] ~ sleep[,2])
test
Welch Two Sample t-test
data: sleep[, 1] by sleep[, 2]
t = -1.8608, df = 17.776, p-value = 0.07939
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-3.3654832 0.2054832
sample estimates:
mean in group 1 mean in group 2
0.75 2.33
＜Chi-squared test＞
χ２検定＝期待値からのずれ(有意差)を検定
> cont <- matrix(c(14,33,7,3), ncol=2)
> cont
[,1] [,2]
[1,] 14 7
[2,] 33 3
> colnames(cont) <- c("Sedan","Convertible")
> rownames(cont) <- c("Male","Female")
> cont
Sedan Convertible
Male 14 7
Female 33 3
> test_chi2 <- chisq.test(as.table(cont))
Warning message:
In chisq.test(as.table(cont)) : カイ自乗近似は不正確かもしれません
> test_chi2
Pearson's Chi-squared test with Yates' continuity correction
data: as.table(cont)
X-squared = 4.1324, df = 1, p-value = 0.04207
＜Wilcoxon signed-rank test＞
符号順位和検定 非正規分布（ノンパラメトリック）の代表値検定
(Wilcoxon rankl sum test(U検定)の方がt検定の非正規分布相当でよく使う)
> x <- c(1.83, 0.50, 1.62, 2.48, 1.68, 1.88, 1.55, 3.06, 1.30)
> y <- c(0.878, 0.647, 0.598, 2.05, 1.06, 1.29, 1.06, 3.14, 1.29)
> test <- wilcox.test(x,y,paired=TRUE,alternative="greater")
> str(test)
List of 7
$ statistic : Named num 40
..- attr(*, "names")= chr "V"
$ parameter : NULL
$ p.value : num 0.0195
$ null.value : Named num 0
..- attr(*, "names")= chr "location shift"
$ alternative: chr "greater"
$ method : chr "Wilcoxon signed rank test"
$ data.name : chr "x and y"
- attr(*, "class")= chr "htest"
> test$p.value
[1] 0.01953125
==================================================================
<Visualization data>
> s1 <- iris[,1]
> p1 <- iris[,4]
> plot(x=p1,y=s1,xlab="Petal length",ylab="Sepal length", col="black", main="Variation of sepal length with petal length")
> boxplot(Sepal.Length~Species, data=iris, ylab="sepal length",xlab="Species",main="Sepal length for different species")
> genex <- c(rnorm(100,1,0.1),rnorm(100,2,0.1),rnorm(50,3,0.1))
> plot(x=genex,xlim=c(0,6),type='l',main="line diagram")
> plot(x=genex,xlim=c(0,20),type='l',main="line diagram")
> x <- rnorm(1000,3,0.02)
> hist(x)
> lines(density(x),col="red")
========p.29==================================================
Working with Pubmed in R
ライブラリ準備
> install.packages("RISmed")
Installing package into ‘C:/Users/ekaminuma/Documents/R/win-library/3.1’
(as ‘lib’ is unspecified)
trying URL 'http://cran.rstudio.com/bin/windows/contrib/3.1/RISmed_2.1.5.zip'
Content type 'application/zip' length 313627 bytes (306 KB)
opened URL
downloaded 306 KB
package ‘RISmed’ successfully unpacked and MD5 sums checked
The downloaded binary packages are in
C:\Users\ekaminuma\AppData\Local\Temp\Rtmp2BIe2S\downloaded_packages
> library(RISmed)
> data(myeloma)
> ls()
[1] "const" "cont" "genex" "myeloma" "p1"
[6] "s1" "sleep" "test" "test_chi2" "x"
[11] "x1" "y"
> str(myeloma)
Formal class 'Medline' [package "RISmed"] with 59 slots
..@ Query : chr "\"multiple myeloma\"[MeSH Terms] AND 2012/05/08[EDAT] : 2013/05/08[EDAT]"
..@ PMID : chr [1:10] "23648714" "23648667" "23648347" "23648290" ...
..@ YearReceived : num [1:10] NA 2013 NA 2013 2013 ...
> Author(myeloma)
[[1]]
LastName ForeName Initials order
1 Kobayashi Tsutomu T 1
2 Kuroda Junya J 2
3 Fuchida Shin-ichi S 3
> ArticleTitle(myeloma)
[1] "The response to second-line induction with bortezomib and dexamethasone is predictive of long-term outcomes prior to high-dose chemotherapy with autologous stem cell transplantation for multiple myeloma."
[2] "Long-term outcome with lenalidomide and dexamethasone therapy for newly diagnosed multiple myeloma."
[3] "[Retrospective analysis on therapeutic effect of autologous hematopoietic stem cell transplantation in multiple myeloma patients]."
> Title(myeloma)
[1] "Internal medicine (Tokyo, Japan)"
[2] "Leukemia"
[3] "Zhonghua yi xue za zhi"
> PMID(myeloma)
[1] "23648714" "23648667" "23648347" "23648290" "23647456"
[6] "23647318" "23647020" "23646139" "23645688" "23645169"
-----------[p.31]----------------------------------------------------------------------------------
NCBI EUtilによるDB検索
> cancer <- EUtilsSummary("cancer[ti]",type="esearch", db="pubmed")
> str(cancer)
Formal class 'EUtilsSummary' [package "RISmed"] with 6 slots
..@ db : chr "pubmed"
..@ count : num 664808
..@ retmax : num 1000
..@ retstart : num 0
..@ PMID : chr [1:1000] "26002954" "26002887" "26002788" "26002780" ...
..@ querytranslation: chr "cancer[ti]"
> sra_plant <- EUtilsSummary("plant", type="esearch",db="sra")
> str(sra_plant)
Formal class 'EUtilsSummary' [package "RISmed"] with 6 slots
..@ db : chr "sra"
..@ count : num 46447
..@ retmax : num 1000
..@ retstart : num 0
..@ PMID : chr [1:1000] "1493794" "1493793" "1493792" "1493791" ...
..@ querytranslation: chr "plant[All Fields]"
---------[p.32]----------------------------------------------------------------------------
> sra_plant@PMID[1:10]
[1] "1493794" "1493793" "1493792" "1493791" "1493480" "1493479"
[7] "1493478" "1493477" "1493476" "1493475"
> cancer.ris <- EUtilsGet(cancer,type="efetch",db="pubmed")
> class(cancer.ris)
[1] "Medline"
attr(,"package")
[1] "RISmed"
> str(cancer.ris)
Formal class 'Medline' [package "RISmed"] with 59 slots
..@ Query : chr "cancer[ti]"
..@ PMID : chr [1:1000] "26002954" "26002887" "26002788" "26002780" ...
..@ YearReceived : num [1:1000] NA NA 2015 2014 NA ...
..@ MonthReceived : num [1:1000] NA NA 2 11 NA NA 3 NA 1 NA ...
..@ DayReceived : num [1:1000] NA NA 11 26 NA NA 4 NA 25 NA ...
..@ HourReceived : num [1:1000] NA NA NA NA NA NA NA NA NA NA ...
..@ MinuteReceived : num [1:1000] NA NA NA NA NA NA NA NA NA NA ...
..@ YearAccepted : num [1:1000] NA NA 2015 2015 2015 ...
..@ MonthAccepted : num [1:1000
