А.Б.Шипунов

Материалы семинара по <strong>R</strong>

Материалы семинара по R

Шипунов А.Б. Материалы семинара по R [Электронный ресурс]. 2002. Режим доступа: http://herba.msu.ru/shipunov/software/r/r-01.htm
Shipunov A.B. R seminar's materials [Electronic resource]. 2002. Mode of access: http://herba.msu.ru/shipunov/software/r/r-01.htm

1  Семинар 1

(а) Работа - с правами "10Б" 
(б) Заведите в директории h:\r4work поддиректорию со своим именем
(в) Перепишите туда из p:\lists\bio\r файлы ws-isl2.dat и spiski.r

1.1  Выход

(0) q()

1.2  Арифметика

(1) sin(pi/2)
(2) sqrt(1024)
(3) log((16^2)/2)
(4) 10^(-5)
(5) round(log((16^2)/2),3)
(6) round(log((16^2)/2),1)

1.3  Присвоения

(7) a <- c(2,4,6,8)
(8) b <- 8:11
(9) a
(10) b
(11) z <- rep(3,4)
(12) z
(13) d <- seq(1,7,2)
(14) d

1.4  Помощь

(15) ?rep
(16) help(seq)

1.5  Объекты

(17) ls()
(18) is.vector(a)
(19) rm(a)
(20) is.vector(a)
(21) ls()
(22) e <- z
(23) e
(24) pjatyj.objekt <- e
(25) pjatyj.objekt
(26) shestoj objekt <- e
(27) a <- c(2,4,6,8)
(28) a
(29) m0 <- cbind(a,b,z,d,pjatyj.objekt,e)
(30) m0
(31) m1 <- log(m0)
(32) m1
(33) m2 <- m1 + m2
(34) m2 <- m1 + m0
(35) m2
(36) is.matrix(m0)
(37) is.vector(m0)
(38) ls()
(39) rm(list=ls())
(40) ls()

2  Семинар 2

Перепишите в рабочую директорию
из p:\lists\bio\r файлы:
fl.dat, sp.dat, dact02s.dat, ws-isl2.dat.

2.1  Чтение

(41) getwd()
(42) setwd("h:\\r4work\\моя_директория")
(43) getwd()
(44) ostr <- read.table
("ws-isl2.dat", head=TRUE, sep=";")
(45) is.matrix(ostr)
(46) is.data.frame(ostr)
(47) names(ostr)
(48) row.names(ostr)
(49) str(ostr)

2.2  Выборки

(50) ostr$ryz
(51) ostr[,32]
(52) ostr[5,]
(53) ostr[1:5,1:5]
(54) t(ostr[1:5,1:5])
(55) ostr.s <- ostr[,order(names(ostr))]
(56) names(ostr.s)
(57) rm(ostr.s)
(58) ostr.n <- names(ostr)
(59) is.vector(ostr.n)
(60) ostr <- as.matrix(ostr)
(61) is.data.frame(ostr)
(62) is.matrix(ostr)
(63) ostr[ostr > 0] <- 1
(64) ostr[1:5, 1:5]
(65) ryz <- ostr[,ostr.n == "ryz"]
(66) izb <- ostr[,ostr.n == "izb"]
(67) sp.ryz <- names(ryz[ryz > 0])
(68) sp.ryz
(69) sum <- ryz + izb
(70) ryz.i.izb <- names(sum)[sum > 1]
(71) ryz.i.izb
(72) length(ryz.i.izb)
(73) source("spiski.r")
(74) spisok(izb)
(75) oba.sp(ryz,izb)
(76) perv.only.sp(ryz,izb)
(77) perv.only.sp(izb,ryz)
(78) ls()

2.3  Запись

(79) save(file="29-09.rd",list=ls())
(80) write.table(ryz.i.izb,
file="ryz&izb.txt", quote=FALSE)
(81) rm(list=ls())
(82) ls()
(83) load("29-09.rd")
(84) ls()

2.4  Описательные статистики

(85) fl <- read.table("fl.dat", h=T, sep=";")
(86) sp <- read.table("sp.dat", sep=";")
(87) names(sp) <- c("VID","NAZ")
(88) fl.1 <- merge(fl, sp, by="VID")
(89) tapply(fl.1$LEPESTOK, fl.1$NAZ, summary)
(90) tapply(fl.1$LEPESTOK, fl.1$NAZ, median)
(91) fl.1$LJUBIT<- fl.1$LEPESTOK %% 2
(92) tapply(fl.1$LJUBIT, fl.1$NAZ, median)
(93) dact <- read.table("dact02s.dat",
h=T, sep=";", dec=",")
(94) str(dact)
(95) summary(dact)
(96) save(file="29-09.rd",list=ls())
(97) q()

3  Семинар 3

3.1  Начало работы

(98) setwd("h:\\r4work\\моя_директория")
(99) getwd()
(100) load("29-09.rd")
(101) ls()

3.2  Функции и перечни

(102) privet<-function() cat("Privet!\n")
(103) privet
(104) privet()
(105) testlist<-list(first=c(1:4),
second=list(sec.1=c(1,2), sec.2=c("A","B")),
third=privet)
(106) testlist
(107) testlist[1]
(108) testlist[[1]]
(109) testlist[[2]][[1]]
(110) testlist$second
(111) testlist$second$sub2
(112) testlist[[3]]()

3.3  Графики

(113) fl.2 <- table(fl.1$NAZ)/length(fl.1$NAZ)
(114) pie(fl.2)
(115) palette(rainbow(length(fl.2)))
(116) pie(fl.2,col=1:length(fl.2))
(117) barplot(rev(sort(fl.2))*100,
names.arg=1:(length(fl.2)))
(118) barplot(table(c("a","b","a","b","a","a")))
(119) names(dact)
(120) attach(dact)
(121) plot(LIP.L)
(122) plot(LIP.L,LIP.W)
(123) identify(LIP.L,LIP.W)
(124) text(locator(), "Sample\n text")
(125) par(mfrow=c(2,2))
(126) boxplot(MIDD.L); boxplot(LATER.L);
boxplot(LIP.L); boxplot(LEAF.L)
(127) par(mfrow=c(1,1))
(128) boxplot(dact[,-c(1:6,10:17)])
(129) boxplot(P.HIGH,LEAF.L,
names=c("Plant's high","Leaf length"))
(130) hist(LEAF.W, breaks=10)
(131) hist(LEAF.W, breaks=20)
(132) srtipchart(dact[,5:8])
(133) coplot(LEAF.L ~ LEAF.W | P.HIGH)
(134) stem(INFL.L)
(135) png(file="ris.png")
(136) coplot(LEAF.L ~ LEAF.W | P.HIGH)
(137) dev.off()
(138) dact.l <- aggregate(dact,
list(LOCALITY), median)
(139) dact.l
(140) pairs(dact.l)

3.4  Генерация данных

(141) y <- rnorm(1000, mean=0, sd=1)
(142) hist(y)
(143) qqnorm(y, main = "normal(0,1)"); qqline(y)
(144) rep(1:10, 100)
(145) f <- factor(rep(1:10, 100))
(146) boxplot(y ~ f)
(147) x <- rbinom(100, size=5, prob=.25)
(148) hist(x, probability=TRUE)
(149) x <- rexp(100, rate=1/2500)
(150) hist(x, probability=TRUE, col=gray(.9))
(151) curve(dexp(x, 1/2500), add=T)
(152) kubik <- function(n) sample(1:6, n, replace=T)
(153) kubik(5)
(154) kubik(12)
(155) kosti <- as.vector(outer(1:6, 1:6, paste))
(156) sample(kosti, 4, replace=T)
(157) sample(kosti, 5, replace=T)
(158) cards <- paste(rep(c(6:10,"V","D","K","T"),4),
c("Tr","Bu","Ch","Pi"))
(159) sample(cards,6)
(160) sample(cards,6)
(161) gl(5, 5, label=letters[1:5])
(162) x <- expand.grid(h=c(60,80),
w=c(100, 300), s=c("A", "B"))
(163) x
(164) unstack(x[,-1])

3.5  Тесты

(165) prop.test(42, 100, p=.5)
(166) prop.test(420, 1000, p=.5)
(167) x <- c(15, 10, 13, 7, 9, 8, 21, 9, 14, 8)
(168) y <- c(15, 14, 12, 8, 14, 7, 16, 10, 15, 12)
(169) t.test(x, y, alt="less", var.equal=TRUE)
(170) t.test(x, y, alt="less")
(171) x <- c(3, 0, 5, 2, 5, 5, 5, 4, 4, 5)
(172) y <- c(2, 1, 4, 1, 4, 3, 3, 2, 3, 5)
(173) t.test(x, y, paired=TRUE)
(174) qqnorm(LEAF.L)
(175) qqnorm(LEAF.W)
(176) quantile(LEAF.L, na.rm=T)
(177) wilcox.test(LEAF.L, conf.int=TRUE)
(178) t.test(LEAF.L, LEAF.W)
(179) table(LIP.DRW, LEAF.SP)
(180) chisq.test(table(LIP.DRW, LEAF.SP))
(181) fisher.test(table(LIP.DRW, LEAF.SP))
(182) summary(table(LIP.DRW, LEAF.SP))
(183) wilcox.test(LIP.DRW, LIP.COL)
(184) wilcox.test(LIP.DRW, LIP.COL)$p.value
(185) round(wilcox.test(LIP.DRW, LIP.COL)$p.value, 3)
(186) kruskal.test(LEAF.W, HERBAR)

3.6  Завершение работы

(187) ls()
(188) save(file="05-10.rd", list=ls())
(189) history()
(190) savehistory("his05-10.r")
(191) q()

4  Семинар 4

4.1  Начало работы

(а) Перепишите в рабочую директорию из p:\lists\bio\r файлы: nuphar99.dat, attrac00.dat, cerat00.dat, nuphar.r, shrtcor.r.

(б) Объяснение значений переменных содержится в заголовке каждого файла. Их можно прочитать, например, просмотрев этот файл в программе Windows Commander при помощи клавиши F3.

(192) setwd("h:\\r4work\\моя_директория")
(193) getwd()

4.2  Работа с данными цветках кубышек

(194) nuph <- read.table("nuphar99.dat", h=T,
sep=";")
(195) str(nuph)
(196) summary(nuph)
(197) stripchart(nuph)
(198) oldpar <- par(mfrow=c(2,2), bg="white")
(199) hist(nuph$SEPALS); hist(nuph$PETALS);
hist(nuph$ANTHERS); hist(nuph$PISTILS)
(200) for (i in 1:4) hist(nuph[[i]],
main=names(nuph[i]), xlab="")
(201) dev.print(png, file="his-nuph.png",
width=600, height=600)
(202) for (i in 1:4) plot(density(nuph[[i]],
na.rm=T), main=names(nuph[i]))
(203) dev.print(png, file="den-nuph.png",
width=600, height=600)
(204) for (i in 1:4) {qqnorm(nuph[i],
main=names(nuph[i])); qqline(nuph[i])}
(205) dev.print(png, file="qqn-nuph.png",
width=600, height=600)
(206) pairs(nuph)
(207) dev.print(png, file="pai-nuph.png",
width=600, height=600)
(208) cor(na.omit(nuph))
(209) nuph1 <- nuph
(210) nuph1$PISTILS[is.na(nuph1$PISTILS)] <-
median(nuph1$PISTILS, na.rm=T)
(211) str(nuph1)
(212) cor(nuph1)
(213) source("shrtcor.r")
(214) shrt.cor(cor(nuph1))
(215) sort(shrt.cor(cor(nuph1)))
(216) attach(nuph1)
(217) cor.test(ANTHERS,PISTILS)
(218) cor.test(ANTHERS,PISTILS, method="spearman")
(219) pet <- cut(PETALS, c(0,16,30),
labels=c("lpet", "hpet"))
(220) pist <- cut(PISTILS, c(0,14,30),
labels=c("lpist", "hpist"))
(221) table(pet,pist)
(222) summary(table(pet,pist))
(223) par(mfrow=c(1,2))
(224) hist(ANTHERS[pet=="lpet"]);
hist(ANTHERS[pet=="hpet"])
(225) wilcox.test(ANTHERS[pet=="lpet"],
ANTHERS[pet=="hpet"])
(226) wilcox.test(ANTHERS[pist=="lpist"],
ANTHERS[pist=="hpist"])
(227) hist(ANTHERS[pist=="lpist"]);
hist(ANTHERS[pist=="hpist"])
(228) par(mfrow=c(1,1))
(229) nuph.lm <- lm(PISTILS ~ ANTHERS)
(230) summary(nuph.lm)
(231) coef(nuph.lm)
(232) plot(ANTHERS,PISTILS,xlab="Number of anthers",
ylab="Number of pistils",main="Nuphar luteum")
(233) abline(nuph.lm, col=2, lty=1)
(234) text(145, 13.5, "PISTILS = 0.05*ANTHERS + 6.62",
cex=.8)
(235) dev.print(png, file="lm-nuph.png",
width=600, height=600)
(236) coplot(PISTILS ~ ANTHERS | PETALS)
(237) dev.print(png, file="cop-nuph.png",
width=600, height=600)

4.3  Автоматическая генерация отчета

(238) savehistory(file="mynuphar.r")

После этого файл можно отредактировать в любом текстовом редакторе для того, чтобы убрать ошибочные и лишние команды. Это можно сделать, например, нажав в Windows Commander кнопку F4.

(239) sink("nuphar.out")
(240) source("mynuphar.r", echo=T)
(241) sink()

4.4  Некоторые полезные команды

(242) edit(nuph)
(243) library(Hmisc)
(244) describe(nuph)
(245) rcorr(data.matrix(nuph))

4.5  Самостоятельная работа

Нужно выбрать любую из двух задач:

(а) Загрузить файл данных по морфометрии роголистника (cerat00.dat), выяснить зависимости между зафиксированными характеристиками. Учтите, что для построения таблиц нумерические данные надо преобразовывать в факторы.

(б) Загрузить файл данных по площади и количеству ``цветков'' (attrac00.dat) и выяснить зависимости между видом растения, типом опылителя, количеством цветков на растение, средней площадью цветка и общей площадью цветков на площадке. Учтите, что для подсчета корреляций надо заменять значения факторов на их коды.

4.6  Завершение работы

(246) ls()
(247) save(file="12-10.rd", list=ls())
(248) history()
(249) savehistory("his12-10.r")
(250) q()

5  Семинар 5

5.1  Начало работы

(а) Перепишите в рабочую директорию из p:\lists\bio\r файл schr02s.dat.

(б) Объяснение значений переменных содержится в заголовке файла.

5.2  Анализ признаков подорожников

5.3  Предварительный анализ

(251) setwd("h:\\r4work\\моя_директория")
(252) getwd()
(253) dir()
(254) file.show("schr02s.dat")
(255) pl <- read.table("schr02s.dat", h=T, sep=";")
(256) str(pl)
(257) summary(pl)
(258) cor(pl)
(259) source("shrtcor.r")
(260) sort(shrt.cor(cor(pl[,-c(1,2)])))

5.3.1  Анализ главных компонент

(261) library(mva)
(262) paste(1:10, names(pl))
(263) pl.pca<-princomp(scale(pl[,c(3:9)], scale=T,
center=T), cor=F)
(264) summary(pl.pca)
(265) loadings(pl.pca)
(266) pl.pcs <- predict(pl.pca)
(267) plot(pl.pcs[,1:2],xlab="PC1",ylab="PC2")
(268) palette(rainbow(max(pl[,1])))
(269) plot(pl.pcs[,1:2],type="n",xlab="PC1",
ylab="PC2")
(270) text(pl.pcs[,1:2],labels=pl[,10],col=pl[,1])
(271) plot(pl.pcs[,1:2],xlab="PC1",ylab="PC2")
(272) identify(pl.pcs[,1:2],labels=paste(pl[,1],
pl[,2],sep=":"),offset=-.5)
(273) plot(pl.pcs[,2:3],xlab="PC2",ylab="PC3")
(274) plot(pl.pcs[,1:3],xlab="PC1",ylab="PC3")
(275) palette("default")
(276) biplot(pl.pca,col=(1:2),cex=c(.5,.8))

5.3.2  Многомерное шкалирование

(277) ls()
(278) ls(pat="pl.p*")
(279) rm(list=ls(pat="pl.p*"))
(280) library(MASS)
(281) pl.dist <- dist(pl[,-c(1,2)])
(282) pl.iso <- isoMDS(pl.dist)
(283) eqscplot(pl.iso$points)
(284) identify(pl.iso$points)
(285) palette(rainbow(max(pl[,1])))
(286) eqscplot(pl.iso$points,type="n")
(287) text(pl.iso$points,labels=pl[,10],col=pl[,1])

5.3.3  Кластерный анализ

(288) pl.h <- hclust(pl.dist)
(289) pl.hw <- hclust(pl.dist, method="ward")
(290) plot(pl.h)
(291) plot(pl.hw)
(292) cutree(pl.h, k=5)
(293) id <- 1:106; attach(pl)
(294) id[STALK/SPIKE>5 & STALK+SPIKE<150 & UPPER<100]
<- "schrenkii"
(295) (id[!(STALK/SPIKE>5 & STALK+SPIKE<150 & UPPER<100)]
<- "maritima")
(296) table(id)
(297) cl <- cutree(pl.h, k=3)
(298) table(id,cl)
(299) summary(table(id,cl))
(300) cl.w <- cutree(pl.hw, k=3)
(301) summary(table(id,cl.w))
(302) summary(table(cl,cl.w))
(303) cl.5 <- cutree(pl.h, k=5)
(304) summary(table(cl.5,cl))
(305) summary(table(cl.5,id))
(306) cor.test(cl,id,method="spearman")

5.3.4  Деревья классификации

(307) library(rpart)
(308) (pl.r1 <- rpart(factor(cl) ~ NEIGH + factor(SUBSTR),
data=pl))
(309) plot(pl.r1,mar=.1)
(310) text(pl.r1)
(311) pl.r1.cl <- predict(pl.r1, type="class")
(312) table(cl,pl.r1.cl)
(313) sum(diag(table(cl,pl.r1.cl)))/length(cl)
(314) summary(table(cl,pl.r1.cl))
(315) cor.test(cl,pl.r1.cl,method="spearman")
(316) (pl.r2 <- rpart(factor(cl) ~ UPPER + STALK + SPIKE
+ FLOWER, data=pl))
(317) plot(pl.r2, mar=.1)
(318) text(pl.r2, cex=.8)
(319) (pl.r3 <- rpart(factor(id) ~ UPPER + STALK + SPIKE
+ FLOWER, data=pl))
(320) plot(pl.r3, mar=.1)
(321) text(pl.r3, cex=.8)

5.4  Завершение работы

(322) ls()
(323) save(file="19-10.rd", list=ls())
(324) savehistory("his19-10.r")
(325) q()




File translated from TEX by TTH, version 3.00.
On 18 Oct 2002, 19:34.