Учитывая, что народ проявил недюжинный интерес к проекту R ( http://molbiol.ru/forums/index.php?showtopic=102358 ), предлагаю здесь организовать русскоязычную консультацию. Принимаются любые вопросы: от "Зачем это все вообще надо, разве не хватает Excel?" до "Как посчитать то-то и то-то?". Надеюсь, что мне помогут отвечать. Цель -- составить русскоязычный FAQ.
Имейте в виду, что существует немало полезных англоязычных ресурсов по R:
1) R FAQ (official) -- http://cran.r-project.org/faqs.html
2) R Help mail list -- https://stat.ethz.ch/pipermail/r-help (non-searchable, updated continuously); http://tolstoy.newcastle.edu.au/~rking/R (searchable, updated once a day)
3) R Help search engine -- http://finzi.psych.upenn.edu/nmz.html
4) R Tips -- http://pj.freefaculty.org/R/Rtips.html
Это лишь небольшая часть.
Поскольку (так я думаю) в русскоязычной среде пользователей R начинающих большинство, давайте сделаем основной акцент на самых базовых вещах, типа установки, импортирования/экспортирования данных, вставления графиков в презентации и пр. Более сложные вопросы тоже принимаются. Я, например, готов отвечать за многомерную статистику.

Я работаю с R уже года два (правда, не очень активно), и тоже готов помогать по мере сил. Работа с R часто состоит из поиска (ну или комбинаторного перебора ) нужного "заклинания", поэтому очень важно иметь возможность оперативно с кем-то посоветоваться, чтобы не наступать на одни и те же грабли.

Вот например, недавно часа 1.5 потратил на выяснение причин, почему не читается текстовый файл. Оказывается дело было в том, что пропущенные значения в файле кодировались как "#N/A", а значок # распознается R как комментарий до конца строки...

У меня такой вопрос: нужно загружать большую матрицу в память. Работает ли R не только с double (8 byte, что и стоит по умолчанию) но и с single (4 байта), а иначе даже 4 Гб памяти не хватает... Возможно ли как-то сконвертировать double в single, и вообще буду благодарен за любые советы по экономии памяти. Может, можно как-то на диск свопить?..

Вместе искать веселее
Вот про double/single:
===
R has no single precision data type. All real numbers are stored in double precision format.
http://finzi.psych.upenn.edu/R/library/base/html/double.html
===
Посмотрите еще здесь -- http://finzi.psych.upenn.edu/R/tmp/Rhelp02...hive/71569.html ,
советуют пользоваться базами данных (может, и вправду?)
R ставит рекорды по скорости вычислений, потому что у него все в оперативке. Так он был изначально задуман. Естественно, у этого есть оборотная сторона...
Вот нашел Вам два пакета:
===
data.table: Just like a data.frame but without rownames, up to 10 times faste
This package does very little. The only reason for its existence is that the white book specifies that data.frame must have rownames. This package defines a new class data.table which operates just like a data.frame, but uses up to 10 times less memory, and can be up to 10 times faster to create (and copy). It also takes the opportunity to allow subset() and with() like expressions inside the []. Most of the code is copied from base functions with the code manipulating row.names removed.
http://cran.arsmachinandi.it/bin/windows/c...a.table_1.0.zip
===
g.data: Delayed-Data Packages
Create and maintain delayed-data packages (DDP's). Data stored in a DDP are available on demand, but do not take up memory until requested. You attach a DDP with g.data.attach(), then read from it and assign to it in a manner similar to S-Plus, except that you must run g.data.save() to actually commit to disk.
http://cran.at.r-project.org/bin/windows/c.../g.data_1.6.zip
===
HTH

Идея с R конечно интересная. Раз уж можно задавать свои вопросы, то я тоже вставлю свои пять копеек. Я думаю, что для того чтобы R действительно пошел в массы, то неплохо было бы иметь описание как на нем производить простейший анализ. Попробовал обработать в R один эксперимент, который ранее в был обсчитан в STATISTICA.

Ввожу данные:
LINE = c('A','A','A','A','A','A','A','A','A','B','B','B','B','B','B','B','B','B',
'P','P','P','P','P','P','P','P','P')
DRUG = c('X','X','X','D','D','D','M','M','M','X','X','X','D','D','D','M','M','M',
'X','X','X','D','D','D','M','M','M')
DOSE = c(3.276,3.226,3.226,5.380,5.388,5.571,6.603,6.452,6.539,3.161,
3.068,3.694,5.170,5.194,5.782,6.857,6.677,7.944,0.596,0.255,
0.286,2.743,2.595,2.774,3.211,4.033,4.377)

Делаю датасет:

civ = data.frame(LINE, DRUG, DOSE)

Далее строю график, просто чтобы посмотреть на данные:

interaction.plot(civ$DRUG,civ$LINE,civ$DOSE)

Делаю ANOVA (без интеракций)

g <- lm(DOSE ~ DRUG+LINE, civ)
anova(g)

А вот вопросы которуе возникли у меня к этому времени:

1. Теперь хочу проверить assumptions. Как в R вывести Normal Probability Plot of Raw Residuals? Такой как в прикрепленном файле.
Нашел qqnorm(g$res), но это не то...

2. Ага и как провести тем Levene в R? Для проверки homogeneity of variances?

3. Есть ли в R другие post-hoc тесты дял ANOVA, кроме Tukey и Sheffe? Например Hochberg GT2, Dunnet T3, Games-Howell (есть в SPSS)?

4. Есть ли аналог Tukey для групп с разным N?

Заранее спасибо за ответы

1. plot(g)
2. levene.test <- function(y, group) {
meds <- tapply(y, group, median, na.rm=TRUE)
resp <- abs(y - meds[group])
table <- anova(lm(resp ~ group))
rownames(table)[2] <- " "
cat("Levene's Test for Homogeneity of Variance\n\n")
table[,c(1,4,5)]
}
3. Нет, но народ советует использовать пакет multtest (например, http://finzi.psych.upenn.edu/R/library/mul....rawp2adjp.html ). Посмотрите еще thread здесь -- http://finzi.psych.upenn.edu/R/tmp/Rhelp02...hive/65720.html И еще вот чего нашел: http://faculty.washington.edu/~jstorey/qvalue
4. ?TukeyHSD (R 2.2.1): "Technically the intervals constructed in this way would only apply to balanced designs where there are the same number of observations made at each level of the factor. This function incorporates an adjustment for sample size that produces sensible intervals for mildly unbalanced designs."

Неужели вот так, "руками", и вводили данные?! Люди могут такое об R подумать...

Хм, у меня эта функция выводит четыре графика, и среди них ни одного P-P Plot, есть только Q-Q Plot Версия пакета - самая последняя - 2.3.1. Что я неправильно делаю что-то может?


Ввел приведенную выше функцию.

А как ею пользоваться в приведенном мною выше примере? Положим я хочу проверить homogeneity of variances для фактора DRUG и LINE. На levene.test(g,LINE) ругается


Вроде как там есть Hochberg, но не уверен что именно тот что есть в SPSS. Тем более что вот в этой ссылке:

написано что эти тесты таки не реализованы

Жалко что этот multtest такое скудное описание имеет - там не расшифровывают. какой тест для каких случаев они предлагают и могут ли они быть полноценной заменой для описанных выше


Это все понятно. Знать бы еще где кончается этот "mildly unbalanced designs" и начинается "severe unbalanced designs". Вот именно для таких случаев и пригодился бы Hochberg.


Ну да, руками и вводил А разве есть способ импортировать файлы STATISTICA?

Так Вам шашечки или ехать? QQ Plot показывает то же самое, и даже эффективнее: http://finzi.psych.upenn.edu/R/tmp/Rhelp02...hive/57078.html

Когда кто-то "ругается", обязательно приводите "ругань"! Иначе отвечать сложно. Но в Вашем случае все просто:
> levene.test(DOSE,factor(LINE))
> levene.test(DOSE,factor(DRUG))

В SPSS-ном виде.

Вот Вам статья про него: http://www.bepress.com/cgi/viewcontent.cgi...text=ucbbiostat

Либо другие multiple testing процедуры. Q-Value, по-моему, весьма интересен -- поскольку он Bayesian, все assumptions ему должны быть по барабану.

Способа нет, потому что STATISTICA серьезные люди не пользуются еще с той давней истории
(см., например, http://www.math.yorku.ca/Who/Faculty/Monet...-stat/0030.html , http://www.math.yorku.ca/Who/Faculty/Monet...-stat/0038.html ).
Зато в R есть способы читать из буфера обмена. Под виндами это совсем просто:
> x <- read.delim("clipboard")
Под всеми платформами можно копировать по одной колонке:
> x <- scan()
Затем вставка, затем 2 раза ENTER.

Почитал я про этот QQ plot. Различия между ними есть, но как я понял его тоже можно для проверки нормальности применять. Для примера взял опять таки те самые данные:

LINE = c('A','A','A','A','A','A','A','A','A','B','B','B','B','B','B','B','B','B','P','P','P','P','P','P','P','P','P')
DRUG = c('X','X','X','D','D','D','M','M','M','X','X','X','D','D','D','M','M','M','X','X','X','D','D','D','M','M','M')
DOSE = c(3.276,3.226,3.226,5.380,5.388,5.571,6.603,6.452,6.539,3.161,3.068,3.694,5.170,5.194,5.782,6.857,6.677,7.944,0.596,0.255,0.286,2.743,2.595,2.774,3.211,4.033,4.377)
civ = data.frame(LINE, DRUG, DOSE)

Провел ANOVA:

g <- lm(DOSE ~ DRUG+LINE, civ)
anova(g)

Теперь хочу проверить нормальность в Shapiro-Wilk normality test:

shapiro.test(g$res)

Shapiro-Wilk normality test

data: g$res
W = 0.9343, p-value = 0.08822


и на QQ plot. Для этого строю его командой:

qqnorm(g$res)
Теперь вопрос - как мне его оценить на нормальность? В PP plot в STATISTICA я знал как это сделать. А тут как?

попробовал вывести линию соответствующую нормальному распределению:
qqline(rnorm(1000))
Получилось вроде что residuals не распределены нормально. (на прикрепленном рисунке)

С другой стороны попробовал комбинацию (которую вроде-бы советуют)
qqnorm(g$res)
qqline(g$res)

Получился график на котором практически все значения лежат на прямой линии, т.е. residuals распределены нормально (второй график).

Так какой из них правильный? Разъясните мне вооще пожалуйста стратегию проверки assumptions в ANOVA в R? Как Вы, например, их делаете?


Попробовал - заработало, спасибо. Но что-то он почему то выдает результат абсолютно отличный от STATISTICA. Вот что выдала STATISTICA:

для LINE:
MS effect ! MS error ! F ! p
0,044755 ! 0,634804 ! 0,070502 ! 0,932119

для DRUG:
MS effect ! MS error ! F ! p
0,032082 0,346457 0,092599 0,911883

А R выдал вот это:
> levene.test(DOSE,factor(LINE))
Levene's Test for Homogeneity of Variance

Df F value Pr(>F)
group 2 0.139 0.871
24

> levene.test(DOSE,factor(DRUG))
Levene's Test for Homogeneity of Variance

Df F value Pr(>F)
group 2 0.0252 0.9751
24

А по идее результат должен одинаковый быть. Где тут собака порылась?

Данные - приведенные в моем примере выше. Процедура Levene для R - написанная Вами выше.

Pr в output R, это я так понимаю "p"?


По post-hoc пока вопросов задавать не будем, еще дойдем до него, давайте сначала с assumptions разберемся.

Жду с нетерпением Ваших ответов по assumptions.

Вот спасибо, это заработало, действительно удобно

1. Про QQ Plot.
Я уже писал, что строить надо командой
> plot(g)
Если Вам нужны не все пять, нужно сделать так:
> plot(g, which=2)
Чтобы получить справку по этой команде, надо сделать так:
> ?plot.lm
потому что R -- по-настоящему объект-ориентированный язык.
2. Про R vs. STATISTICA
Можно задать про это вопрос в R Help, а можно автору кода "Levene Test" -- John Fox ( http://socserv.mcmaster.ca/jfox ). Но боюсь, что без особых резутьтатов -- про отношение к STATISTICA среди статистиков я уже писал, а самое главное -- Вы никогда не узнаете, что действительно использует STATISTICA, тогда как в случае с R это очень просто.
Кстати, если уж про тесты:
> ?bartlett.test
и все тесты, которые там в "See Also".
ANOVA пользуюсь нечасто, штатных средств R мне вполне хватает.

Хорошо, справку почитал, но не нашел там ответа на свой вопрос Имеет ли значение градус наклона линии на QQ plot или нет? Т.е. если все точки (или по крайней мере подавляющее большинство) располагаются на линии или очень близко к ней, то данные можно считать нормальными? Или наде еще смотреть на угол наклона линии на этом графике и насколько он соответсвует линии нормального распределения (т.е. y=x, 45 градусов, если шкалы по x и y равны)?


Разобрался сам. Дело в том, что STATISTICA считает Levene test на основе means (доступа к ее исходным кодам у меня, понятно нет, но установил это опытным путем). А John Fox реализовал levene.test в R на основе medians. Отсюда и некоторая разница в p values.

А проверил это так, что в процедуре John Fox'а заменил median на mean:

levene.test <- function(y, group) {
meds <- tapply(y, group, mean, na.rm=TRUE)
resp <- abs(y - meds[group])
table <- anova(lm(resp ~ group))
rownames(table)[2] <- " "
cat("Levene's Test for Homogeneity of Variance\n\n")
table[,c(1,4,5)]
}

и тест в R начал выдавать результат идентичный со STATISTICA

Почитал интернет по этому вопросу и оказалось, что второй вариант (реализованный в R, через medians) - точнее (например тут: http://tolstoy.newcastle.edu.au/R/help/03b/1901.html). Не то чтобы я не доверял R, но эта разница в p values меня тревожила. Теперь, когда причина разницы между STATISTICA и R найдена, все стало на свои места.

Вот еще инфа о тестах Homogeneity of Variances, которую нашел в своих конспектах, переписанная с какой-то книжки: "Hartley's F-max test was examined when discussing the two-sample t-test. It is simply the max variance divided by the min variance. Bartlett's test is computationally a bit more prolonged. A nice description and worked example can be found in Sokal & Rohlf (1995; Box 13.2). Both Hartley's and Bartlett's are sensitive to departures from normality, so this needs to be determined first. Further, Hartley's test requires equal sample sizes. The Scheffé-Box test is less sensitive to departures from normality and can be used for unequal sample sizes, but requires the data to be acquired in a stratified group fashion. Perhaps the best overall test (because of its insensitivity to sample size and normality) to examine the homogeneity of variance assumption is the Modified Levene Equal Variance test. Here, all variates are redefined by subtracting the median of each subgroup and running a one-way ANOVA on theses redefined variates. If you fail to reject the null hypothesis, conclude that variances are equal."

Чисто познавательно, интересно какой из вариантов Levene реализован в SPSS... Но что-то мне не хочется демо ставить только ради этого

А вообще в ходе поисков обнаружил что в пакете ctest есть тест, который некоторыми считается еще более robust чем Levene:

Fligner-Killeen Test of Homogeneity of Variances. The Fligner-Killeen (median) test has been determined in a simulation study as one of the many tests for homogeneity of variances which is most robust against departures from normality, see Conover, Johnson & Johnson (1981). It is a k-sample simple linear rank which uses the ranks of the absolute values of the centered samples and weights a(i) = qnorm((1 + i/(n+1))/2). The version implemented here uses median centering in each of the samples (F-K:med X^2 in the reference).

Последнюю реплику не совсем понял. А разве ANOVA не является штатным средством R? Я думал это один из самых часто используемых параметрических методов? А чем сравниваете несколько групп, если такая нужда есть? Или сразу регрессию проводите? Но там ведь assumptions тоже надо сравнивать...


Кстати, Вы может подскажите мне, как заставить R выводить p values в десятичных дробях, а не в научном формате? Т.е. 0.000002238 вместо 2.238e-06 ?

Да

Нет

Отлично! И спасибо за инфу о тестах.

Это я неудачно выразился. Конечно, несколько групп сравниваю с помощью ANOVA.

http://finzi.psych.upenn.edu/R/Rhelp02a/archive/15472.html

А какими тестами/графиками assumptions проверяете?

Прочитал, спасибо.
Только там форматировать предлагают каждый раз - это немного неудобно для меня. Я бы предпочел глобально это один раз переопределить для всех репортов. Решение такое нашел в доках:
options(scipen=20)
правда длинные они тогда иногда числа получаются, но все равно, для моих eyeballs оно так нагляднее и приятнее

У меня еще несколько вопросов накопилось.

1. Вы строете графики в R или GNUPlot? Какой из них мощнее?

2. Как в R на график впихнуть p values с какого либо теста?
Вот например сделал я ANOVA, вывел QQ plot командой
plot(g, which=2)

сделал тест shapiro-wilk на нормальность:
shapiro.test(g$res)

выдал он мне такой результат:
Shapiro-Wilk normality test

data: g$res
W = 0.9343, p-value = 0.08822

Как его заставить еще на моем QQPlot еще типа такой надписи написать:
"Shapiro-Wilk normality test: p=0.08822"

Причем как дополнительную надпись, а не вместо какого-нибудь титула?

3. Как можно заставить R каждый новый график открывать в новом окошке, а не стирать старый?

QQ Plot, Levene test

Gnuplot давно забросил. R, конечно же, мощнее. Вот, для забавы, посмотрите график:

и код к нему http://www.stat.auckland.ac.nz/~paul/RGrap...examples-once.R
В R не очень много интерактивности в графиках (один из примеров ниже), но есть замечательный GGobi ( http://www.ggobi.org ), который сию интерактивность добавляет.

Извиняюсь за выражение,
> text(locator(), paste("Shapiro-Wilk normality test: p=", round(shapiro.test(g$res)$p.value, 5), sep=""), pos=4)
Щелкните мышкой куда хотите вставить, а затем Esc.

Перед графиком
> x11()
либо
> windows()
(менее универсально), на Mac --
> quartz()

А чем пользуетесь для постройки SVG графиков (если строите)? gridSVG? Хотелось бы чего нибудь scalable, чтобы не пихать в документ статичные пикчерсы в PNG. Хотябы какая-та замена динамическим графикам STATISTICA, которые можно в любой момент поправить/изменить. Хотя ни MSOffice, ни OpenOffice все равно SVG не понимают...



Господи, как только люди не извращаются Пришлось для этого даже R Graphics установить.


Скачал, поставил, он начал орать на отсутствие какой-то GTKшной либы, которой почему то в дистрибутиве не было. Пришлось ему аналогичную из ГИМПа скормить

А есть к нему какой-нибудь туториал для чайников? Такой который можно безболезненно прожевать?


Я предпочитаю автоматизировать - даю ему сразу координаты. Только заметил интересную особенность - он надпись пихает не по абсолютным координатам, а по координатам осей - т.е. внутри графика. И не в какую не хочет его помещать снаружи собственно области графика. Как это побороть можно? Так чтобы надпись дать например вверху, рядом с титулом?


Мои знания по R благодаря Вам расширяются просто со пугающей скоростью. Аж самому страшно А ведь помню когда его первый раз поставил - поковырялся, да и снес от греха подальше... А оно оказывается совсем и ничего

Графики в SVG не строю, хотя формат очень перспективен, вот FireFox с 1.5 имеет нативную его поддержку, Inkscape тот же развивается семимильными шагами... Но ни одно приложение, кроме Adobe Illustrator, не имеет пока полной и безглючной его поддержки
Я строю PDF -- на Mac OS X это особенно хорошо, потому что его понимают _все_ приложения (ведь Quartz -- это DisplayPDF), да и под виндами с PDF вполне ничего. Редактировать, правда, опять же в AI. Можно EPS, редакторов больше (Mayura, к примеру), но больше и трабла. Если не слезать с виндов, можно WMF (редакторов полно). R умеет еще picTeX (редактировать ручками, поскольку TeX) и xfig (редактировать при помощи jfig). Все названное -- scalable.

Дык, это ж Paul Murrell, который с нуля сделал Trellis graphics в R!

У них на сайте есть даже учебные фильмы для чайников, поищИте.

Естественно, особенно если Вы предпочитаете автоматизировать! Поменяются точки, и Ваша надпись закроет что-нибудь ценное...

?mtext

Работать с R несколько муторно, но никуда не деться т.к. многое есть только здесь. Начал потихоньку осваивать геостатистический анализ (вариограммы, кригинг). Пробовал в 3-х пакетах посчитать доверительные интервалы для малых выборок с помощью бутстрэпа, но на одной выборке все пакеты почему-то сбоили. В R получилось. Здесь бутстрэп есть в нескольких библиотеках, но успешно посчитал толко в модуле boot.ci библиотеки boot.
> Plantago
Не знаете ли, есть в R возможность провести факторный анализ с какой-нибудь ресэмплинг-техникой: бутстрэпом или перекрестной оценкой (n-1 fold cross-validation)? Не главные компоненты, а именно факторный анализ (конктретный метод менее важен - нужны вращения). Дело в том, что интересные результаты получил на очень малых для R-техники выборках (10-15 случаев при 7-10 переменных) и теперь нужна хотя бы формальная оценка надежности решения.

В пакете sem есть bootstrap для structural equation models, что близко к FA. Вы можете выдать свой объект за sem-объект или сделать по аналогии, посмотрев код пакета.
Пакет boot предназначен для применения bootstrap к каким угодно задачам, поэтому Вы можете сами написать (думаю, не очень сложный) код под свою задачу. Вот здесь ( http://cran.at.r-project.org/doc/Rnews/Rnews_2002-3.pdf ) есть статья A. Canty о том, как это делать. Поскольку R все автоматизирует, любые resampling-техники достаточно легко пишутся руками. Я вот сам писал для кластерного анализа, но потом перешел на более продвинутый пакет pvclust. Главная трудность там была -- разобраться в структуре hclust-объекта. Объекты PCA/FA имеют куда более простую структуру, так что особых проблем быть не должно.
Посмотрите еще сюда:
http://tolstoy.newcastle.edu.au/R/help/05/02/11520.html
Объяснено на примере PCA, но, думаю, с FA будет похоже.

Большое спасибо, попробую разобраться!

Спасибо.

У меня к Вам еще вопросы

1. Скажите, а как можно при постройке interaction.plot указать свой порядок сортировки значений?

Данные ввожу в таком виде:

LINE = c('A','A','A','A','A','A','A','A','A','
B','B','B','B','B','B','B','B','B',
'P','P','P','P','P','P','P','P','P')
DRUG = c('X','X','X','D','D','D','M','M','M',
'X','X','X','D','D','D','M','M','M',
'X','X','X','D','D','D','M','M','M')
DOSE_ORIG = c(0.0656,0.0624,0.0624,0.5381,0.5420,0.6512,1.8270,1.5720,
1.7140,0.0585,0.0533,0.0997,0.4359,0.4468,0.8040,2.3550,1.9670,6.9890,
0.0045,0.0032,0.0033,0.0385,0.0332,0.0397,0.0615,0.1399,0.1973)
civ = data.frame(LINE, DRUG, DOSE_ORIG)


Потом строю график функцией:

interaction.plot(civ$DRUG,civ$LINE,civ$DOSE_ORIG)

Но он почему-то фактор DRUG (по оси X), выводит в последовательности D, M, X. Хотя при вводе данных, как видно выше я вводил в последовательности X, D, M Как его заставить plot строить тоже для оси x в последовательности X, D, M? Пример вывода графика R прикрепляю к посту.

3. Как заставить его смещать разные линии по оси x на небольшую величину, чтобы они не накладывались друг на друга?

4. Как заставить его отображать ось y в логарифмической скале?

5. Как сделать чтобы на графике отображались также SD, а не только means?

6. Как его заставить выводить этот plot не в lines, a в bars (как это показано на прикрепленном графике из STATISTICA)?

7. Как заставить R весь вывод (причем не только text output (результаты тестов), но также построенные графики) перенаправлять в какой-либо файл сразу? Очень желательно документ MSOffice/OpenOffice чтобы можно было подредактировать, скомпоновать и распечатать или на худой конец хотя бы в PDF?

Заранее спасибо за ответы

Извините, что отвечаю поздно и мало -- много работы...
1) > civ$DRUG <- factor(civ$DRUG, levels=c("X","D","M"))
> interaction.plot(civ$DRUG,civ$LINE,civ$DOSE_ORIG)
3) Не понял, объясните подробнее.
4) > interaction.plot(civ$DRUG,civ$LINE,civ$DOSE_ORIG, log="y")
5) Только sd: > interaction.plot(civ$DRUG,civ$LINE,civ$DOSE_ORIG, fun=sd)
Если хотите в одном флаконе, то надо что-то паять при помощи plot.new() или lines() -- чтобы R не открывал новый девайс, а писал прямо на старый.
6) А зачем? Удобнее?
Так не нравится: > interaction.plot(civ$DRUG,civ$LINE,civ$DOSE_ORIG, type="p", pch=21:24) ?
interaction.plot() -- вещь простая, но заточенная под линии/точки, ежели желаете столбики, надо просто его сымитировать через barplot()
7) Есть sink(), есть pdf(), есть win.metafile() и другие девайсы. Можно их комбинировать. R, запущенный без графического вывода (по-моему, опция --vanilla) складывает картинки в один многостраничный PostScript-файл. Есть Sweave(), но для этого надо быть хоть сколько-нибудь знакомым с TeX/LaTeX (хотя пристойные конверторы в OpenOffice существуют). В Вашем случае, наверное, лучше всего пакет r2html ( http://finzi.psych.upenn.edu/R/library/R2HTML/DESCRIPTION ) Можно еще Rpad ( http://www.rpad.org/Rpad ), но он, по-моему, уж слишком наворочен.

Кстати, не обратили внимание -- столбики на Вашем втором графике кажутся наклоненными друг к другу из-за зрительного эффекта

Да ничего страшного, спасибо что вообще отвечаете, а то больше сюда никто не приходит


Ага, именно это!


Ну так как на графике прикрепленном к этому письму. Видите, там линии немного смещены относительно друг друга и за счет этого их SD whiskers не накладываются друг на друга.


Да, это тоже именно то что надо


Да, для меня нагляднее и привычнее. Нашел в галерее как раз тот тип столбиков что мне нужен, установил этот пакет (gregmisc) с функцией barplot2. Но никак не могу прикрутить чтобы он мне построил нечто подобное с моими данными. Не могу я понять, чего куда там посылать, неужели все это как-то проще нельзя сделать - просто дать ему в какой колонке grouping factor(s) а в какой dependent variables, чтобы он мне посчитал mean + SD и выдал столбики с error bars.

Вот такой график мне и нужен:
http://addictedtor.free.fr/graphiques/RGra...ry.php?graph=54

Там по x - civ$DRUG хочу дать (причем в каждом civ$DRUG он еще разделяется на три столбика civ$LINE), а по y - соответственно mean для civ$DOSE_ORIG +/- SD.

Вот как сделать столбики:
===
c.mean <- tapply(civ$DOSE_ORIG, list(civ$LINE, civ$DRUG), mean)
c.sd <- tapply(civ$DOSE_ORIG, list(civ$LINE, civ$DRUG), sd)
c.up <- c.mean + c.sd
c.down <- c.mean - c.sd
barplot2(c.mean, beside=T, legend.text=c("A","B","P"), log="y", ci.l=c.down, ci.u=c.up, plot.ci=T, col=rainbow(3))
===
Я уж Ваши данные в скрипт загнал, чтобы быстрее было

Кажется, понял Ваш вопрос про разнесение точек!
Это то, что нужно?
===
x <- c(rep(3,3), rep(6,3), rep(9,3))
y <- c(rep(5,3), rep(8,3), rep(11,3))
plot(x, y)
plot(jitter(x), jitter(y)))
library(gregmisc)
plot(space(x,y))
===
Либо jitter(), либо space() из библиотеки gregmisc

Hello!

Ja tolko nachinau rabotat s R i u menja est vopros. Pravda, ne uverena, chto eto napryamuyu otnositsya k R/

Podskajite, kak cozdavat v R paketi (packages) dlya Bioconductor. Chto, krome R coda, dlya etogo neobhodimo?

Pri proverke sozdannogo paketa (s pomowju R CMD check) - voznikaet owibka, cvyazannaya s tem, 4to paket ne mojet bit installirovan.

Esli kto-to s etim stalkivalsja - pomogite!

Spasibo.

Сталкивался немного. Чтобы ответить на Ваши вопросы, мне нужна дополнительная информация. Объясните, пожалуйста, для чего Вы создаете пакет? У Вас есть значительный объем нового кода, который Вы желаете распространять? Или какая-то другая причина?

Kod uge napisan, i oformlen f funkcii, Teper' hotim sdelat' iz etogo paket. Sobsrvenno ia uge poniala v 4em bila oshibka, tak 4to prodolgaju rabotat' Izvinite za bespokoystvo.

Отлично, что Вы нашли причину! Вообще говоря, по разработке пакетов есть классическое руководство http://cran.r-project.org/doc/manuals/R-exts.html , а если Вы работаете под винды, то очень полезно вот это: http://cran.r-project.org/doc/contrib/Wang-WinBook.pdf

Plantago, большущее Вам спасибо опять. Я как раз не знал как мне из нескольких векторов сделать таблицу. Теперь с Вашей помощью все стало ясно. Я немного его изменил, под свои нужды, на основе Вашего примера и примеров из сети. Вот что получилось (на barplot3.png):

===

LINE = c('A','A','A','A','A','A','A','A','A','B','B','B','B','B','B','B','B','B',
'P','P','P','P','P','P','P','P','P')
DRUG = c('X','X','X','D','D','D','M','M','M',
'X','X','X','D','D','D','M','M','M',
'X','X','X','D','D','D','M','M','M')
DOSE_ORIG = c(0.0656,0.0624,0.0624,0.5381,0.5420,0.6512,1.8270,1.5720,
1.7140,0.0585,0.0533,0.0997,0.4359,0.4468,0.8040,2.3550,
1.9670,6.9890,0.0045,0.0032,0.0033,0.0385,0.0332,0.0397,0.0615,0.1399,0.1973)

civ = data.frame(LINE, DRUG, DOSE_ORIG)

civ$DRUG <- factor(civ$DRUG, levels=c("X","D","M"))
civ$LINE <- factor(civ$LINE, levels=c("P","A","B"))

# Make bar plot with SD

c.mean <- tapply(civ$DOSE_ORIG, list(civ$LINE, civ$DRUG), mean)
c.sd <- tapply(civ$DOSE_ORIG, list(civ$LINE, civ$DRUG), sd)
c.up <- c.mean + c.sd
c.down <- c.mean - c.sd

windows()
pcolors <- c("lightblue", "mistyrose", "lavender")
mp <- barplot2(c.mean, beside = TRUE, col = pcolors,
log="y", ylim = c(0.001, 10), ylab = expression(mu*"g/ml"),
main = "Sample title", font.main = 4, sub = "Sample graph",
cex.names = 1.5, ci.l=c.down, ci.u=c.up, plot.ci=TRUE,
plot.grid = TRUE)
smartlegend(x="left",y="top", inset = 0.03, rownames(c.mean), fill = pcolors)
box()
===

Ляпота!


Да, Вы правильно уловили что именно мне нужно. Именно функция space мне бы подошла больше чем jitter(), так как она каждое значение разносит на постоянную величину. Но как заставить ее работать с категоризированными данными, например когда у меня по оси х – три кода лекаства (буквы), а не значения. Так, например, как в этом примере:

===
LINE = c('A','A','A','A','A','A','A','A','A','B','B','B','B','B','B','B','B','B','P','P','P','P','P','P','P','P','P')
DRUG = c('X','X','X','D','D','D','M','M','M',
'X','X','X','D','D','D','M','M','M','X','X','X','D','D','D','M','M','M')
DOSE_ORIG = c(0.0656,0.0624,0.0624,0.5381,0.5420,0.6512,1.8270,1.5720,1.7140,0.0585,
0.0533,0.0997,0.4359,0.4468,0.8040,2.3550,1.9670,6.9890,0.0045,0.0032,0.0033,0.0385,0.0332,0.0397,0.0615,0.1399,0.1973)

# Make dataset

civ = data.frame(LINE, DRUG, DOSE_ORIG)

# Make correct order of the factors

civ$DRUG <- factor(civ$DRUG, levels=c("X","D","M"))
civ$LINE <- factor(civ$LINE, levels=c("P","A","B"))

# Display plot to see data graphically

interaction.plot(civ$DRUG,civ$LINE,civ$DOSE_ORIG,log="y")
===

Там, например, линии накладываются друг на друга, по x, а я хотел бы чтобы они немножко расходились, как в вашем примере со space. Можно ли как нибудь эту функцию срастить с interaction.plot? Так чтобы получилось нечто похожее как на картинке int.jpg (внизу)


У меня появился к Вам еще другой вопрос, немного другого плана:

Есть некоторый ряд данных. Представляют собой IC50, т.е. концентрации некоторого вещества, которые ингибируют пролиферацию 50% клеток. Данные получены из несколько тестов, например 0.801, 0.902, 0.850, 0.760, 0.810 microg/ml. Надо посчитать среднее и SD. У нас в лаборатории из этого считают геометрическое среднее, мотивируя что дозы IC50 расчитываются по графику на котором по x – логарифмы дозы (например для 100, 10, 1, 0.1 microg/ml будет 2, 1, 0, -1), а по y – проценты ингибирования от 0 до 100 (линейная шкала). Но у меня возникли определенные сомнения, может надо все-таки обычное среднее + SD считать? Насколько я знаю геометрическое среднее считается для процентов, долей или для данных, которые распределены сильно ненормально, например: 1, 2, 5, 6, 1500. А эти наши данные хотя и отчитываются по логарифмической шкале, но после преобразования опять представляют собою нормальные концентрации. Как вы считаете?

Ого, как у Вас здорово все выходит!
С учетом этого не буду предлагать готового решения. Попробуйте сделать так: выньте из R код interaction.plot() -- для этого достаточно ввечти эту команду без скобок, а потом засуньте внутрь туда куда-нибудь space().
Второй вопрос для меня практически ясен: то, как работать с данными, зависит прежде всего не от того, КАК они были получены, а от того, ЧТО они собою представляют. С этой точки зрения правильно, конечно, вычислять обычные средние.

Плантаго, не мучайте меня

Попробовал так как вы сказали:
interaction.plot(space(civ$DRUG,civ$LINE),civ$DOSE_ORIG,log="y")

На это он ругнулся:

Error in Summary.factor(..., na.rm = na.rm) :
max not meaningful for factors

Я так подозреваю потому что civ$DRUG и civ$LINE - не численные факторы, а категориальные. Есть какая-нибудь space для категориальных данных может?


Да, целая дискуссия сегодня развернулась. Пришли к выводу, что наверно в наших условиях обычное среднее и стандартное отклонение все-таки более подходит, потому как у нас данные не расходятся на несколько порядков. Да и считаются средние для 3-5 значений, так что все равно обычная средняя выходит похожа на геометрическую. А если считать геометрическую, то там геометрическое SD не показательно, надо считать CI95%, которые, насколько я помню там тоже немного иначе считаются. Короче наверно обычные mean+SD таки удобнее.

Я не мучаю, я сам не знаю (честно). Времени просто нет делать длительное исследование вопроса. А Вы все-таки попробуйте написать "interaction.plot" без скобок и нажать Enter...

Попробовал, это мне выдало исходный код функции Жалко конечно, но буду копать... Если найду - обещаю поделиться с общественностью тут.

Кстати, может у Вас есть доступ к следующим книжкам:
1. “R Graphics"
by Paul Murrell
http://www.stat.auckland.ac.nz/~paul/RGrap.../rgraphics.html

2. “Statistics: An Introduction using R”
Michael J. Crawley, March 2005
http://eu.wiley.com/WileyCDA/WileyTitle/pr...escription.html

3. Survival Analysis using S,
Tableman, Kim & Portnoy, Chapman & Hall, 2003
http://www.jstatsoft.org/v11/b05/v11b05.pdf

4. An R and S-PLUS Companion to Applied Regression
John Fox, Sage Publications, 2002
http://socserv.mcmaster.ca/jfox/Books/Companion/index.html

5 Introductory Statistics with R
by Peter Dalgaard, Springer, 2002.
http://staff.pubhealth.ku.dk/~pd/ISwR.html

6. Может Вы знаете. Есть две книжки:

“Practical Regression and Anova using R”, издание 2002 (http://www.stat.lsa.umich.edu/~faraway/book/)
“Linear Models with R”, издание 2005 (http://www.stat.lsa.umich.edu/~faraway/LMR/)
обе написаны Julian Faraway.

Но первая – free, а вторая – за деньги. Судя по содержанию, вторая – просто переработанное переиздание первой. Причем по количеству страниц почти не отличаются. Вы в курсе, а насколько они отличаются в плане полноты и доступности информации? “Practical Regression and Anova using R” еще полностью не прочитал, но первое впечатление – хорошее. Вопрос – стоит ли покупать более новое издание за деньги?

И еще. В начале этого поста Вы предлагали написать русский FAQ по R. Мне кажется это преждевременно. Потому как FAQ пишется когда кол-во пользователей набирает какую-то критическую массу и они начинают задавать однотипные повторяющееся вопросы. А русскоязычных пользователей R (профессиональных математиков и статистиков оставим за скобками, им FAQ не нужен) - пока еще исчезающе мало в том море Excel, STATISTICA и SPSS. Мне кажется что бы действительно пригодилось на данном этапе - простенькая книжка с самыми основами R, с описанием на примерах как провести в R простейший и наиболее часто встречающийся статистический анализ: t-test, ANOVA, непараметрическая статистика, описательные статистики, постройка графиков, анализ выживаемости, etc. Этим можно сразу здорово популяризовать R в русскоязычной среде. Потому как русководств таких уже достаточно много на английском, но на русском - нет.

Ну, я Вам сам накопал. См. приложение. Правда, ни space(), ни jitter() не подошли, пришлось писать самому.

Есть ксерокс предыдущей его книги про S-PLUS, весьма хорошая вещь. Он эколог, кстати, и большая часть его примеров -- биологические.

Есть на полке.

Книжки очень похожи, так что покупать за деньги не стоит.

Очень здравое мнение. Я, собственно, еще в конце позапрошлого года сделал такую попытку, и даже заключил договор с некоторым компьютерным издательством. Но пороху, увы, не хватило, написал только одну главу Я теперь всем подряд предлагаю соавторство, вот и Вам тоже Может, напишем вместе? Могу прямо тут вывесить план книжки, можно обсудить.
Есть и другой вариант: просто перевести книжку Dalgaard, в которой есть все то, чего Вы хотите. И, наконец, можно договориться с Peter (по переписке он производит впечатление вменяемого человека) и сделать гибрид: часть книжки как перевод, а часть написать самим. Это мне нравится больше, потому что русскоязычного читателя надо гораздо дольше вводить в курс дела: рассказывать про типы данных, про импорт/экспорт, а я, кроме того, очень хочу написать про многомерные методы.

Снимаю шляпу! Это именно то что нужно. Я немного модифицировал Ваш код,
добавил возможность отображения SD (часто полезно видеть, чтобы оценить разброс в группах) и исправил ошибку, которая появлялась когда количество rows и columns неодинаковое в табличке, по которой строим график. Помещаю тут, может пригодится кому (архив внизу). Вот такая красотень теперь выходит (иллюстрация внизу):


Спасибо за ответы по книжкам. А какую из них Вы бы советовали в качестве пособия по статистике и началам анализа в R для студентов естественнонаучников:

2. “Statistics: An Introduction using R”
Michael J. Crawley, March 2005

или

5 Introductory Statistics with R
by Peter Dalgaard, Springer, 2002.?

Необходимый минимум: описательные статистики, тесты для двух проб, пропорций, ANOVA, непараметрическая статистика, хорошее введение в анализ выживаемости. Больше склоняюсь к Crawley, так как он поновее, потолще и содержание мне больше понравилось. Как Вы считаете?


Предложение, конечно, заманчивое. Да только я не назвал бы себя специалистом в статистике, хотя она мне и нравится. Я скорее потребитель, да и знакомство с R у меня ограничивается двумя неделями Может как изучу его поплотнее, тогда можно будет вернуться к этой теме, потому как в принципе идея интересная. А как Вы считаете, она бы пользовалась спросом? Целевую аудиторию какую планируете?

Накопилось у меня еще пару вопросов к Вам:

1. Где в R (или может в Tinn-R) можно посмотреть табличку с цветами + их назван ия? Удобно было бы выбирать при постройке графиков. Про colors() знаю, но она только названия дает.

2. Аналогично, где можно посмотреть коды R для всяческого рода символов - греческие буквы, математические выражения, etc. Чтобы можно было потом в expression() подставлять.

3. Читая дискуссию http://tolstoy.newcastle.edu.au/R/help/05/08/11237.html наткнулся на такой пакет - coin (http://www.maths.lth.se/help/R/.R/library/.../doc/index.html). Весьма занятная вещь, Вы может работали с ним? Правильно ли я понимаю, что с его помощью можно производить непараметрический двухфакторый анализ, в случаях когда интеракции не предполагаются (не важны)? Если я правильно понял, то это получается вроде грубой замены two-way ANOVA without interactions когда данные не подпадают под нормальное распределение?

Спасибо за улучшенный код!
Что касается книжек, то я посоветовал бы Dalgaard -- весь минимум там есть, а чем книжка тоньше, тем студентам лучше Crawley, скорее, для аспирантов.
Встроенных таблиц с цветами и символами в R нет, надо искать в Сети, по книжкам, или делать самостоятельно. Впрочем, я еще посмотрю.
Пакетом coin не пользовался, так что ничего сказать не могу.
Свою книжку планировал для широкой аудитории пользователей статистических методов, для тех, кто не боится компьютера. Мучительно пытался решить, не начать ли с RCommander -- так можно работе в R обучить абсолютного чайника. Все же решил начать с командной строки

Plantago,

Я решил немножко оживить тему. Для затравки задам Вам вопрос - что Вы посоветуете для power/sample size calculation в R? Какие пакеты? В стандартном stats нашел только powr.t.test, power.anova (причем только one-way) и вроде для chi-test. Есть какой-нибудь может специализированный пакет?

А то для SAS есть UnifyPow (http://www.bio.ri.ccf.org/Power/) c очень богатым выбором поддерживающих методов, но в R он не работает

UnifyPow принципиально не обновлялся с 1998 года.
Для R есть asypow, pwr, ssanv.

Спасибо, из них трех pwr мне показался самым оптимальным для моих задач.

А еще другого плана вопрос, Вы не могли бы подсказать, как правильно давать ссылку на R в публикации при описании статистического анализа в материалах и методах. Где-то видел это уже, а где - забыл

И еще, если Вы в анализе использовали какие-то сторонние пакеты, а не только стандартные, то их также нужно упоминать? Если да, то как - имя, версия, ссылка в инете?

Я обычно цитирую R и не цитирую пакеты, хотя это, наверное, неправильно. Есть хорошая команда citation() , где все написано.

А работал ли кто-нибудь с данными microarray в пакете R?
Я сама только начинаю, то есть опыт практически нулевой, а тратить время зря не хочется...
Может быть кто-то даст мне советы, с чего начать, какие пакеты установить и тд и тп.
Наверное и Bioconductor тоже надо использовать будет, так что любая информация о этой софтине также будет полезна.

Bioconductor -- это набор пакетов к R. К сожалению, с microarray никогда не работал Но более широкие справки дать, наверное, могу. Так что пишите.

http://molbiol.ru/forums/index.php?showtop...960#entry423960
Dalgaard_P._Introductory_Statistics_With_R_

http://www.stat.berkeley.edu/~terry/zarray...re/smacode.html

хттп://лудщиг-сун2.унил.ч/~дарлене/Рмини/счед.хтмл
Добавка
есть еше такая штучка :-)
http://bioinf.ucd.ie/people/aedin/R/

Мысль действительно очень здравая. Небольшой толковой книги на русском по R (введения) действительно катастрофически не хватает. Периодически наталкиваюсь на попытки отдельных людей написать руководство, но все они переходили в стадию долгостроя (одному такое потянуть действительно тяжело).

А еще больше не хватает русскоговорящего сообщества (кому можно задать вопрос и оперативно получить ответ).

Начинание очень нужное и полезное, однако в виде форума нецелесообразно.
Форум очень быстро разрастется, к тому-же он линейный. Поэтому предлагаю создать сайт на wiki-движке.

На первое время собрать там ссылки на английские ресурсы и обрывки того, что есть на русском (в том числе и ссылки на обзорные статьи по R). Затем начать переводить небольшие статьи (главы) и параллельно работать над предложенной книгой.

На себя готов взять заботу о хостинге (регистрация имени, размещение, поддержка) и помощь с переводом.

Если есть люди, которые смогут помочь в настройке Wiki-движка, то я могу заняться этим уже с 1 января.



Предлагаю свою помощь по переводу.
Абсолютно согласен, что надо давать и введение в статистику (типы данных, импорт/экспорт и т.д.)

Кое-что есть у меня на сайте ( http://www.r-project.org/other-docs.html , внизу ), вчера выложил книжку по основам статистики "для маленьких" ( http://herba.msu.ru/shipunov/school/sch-ru.htm , часть про R сейчас пишется). Wiki люблю за идею, но положительного опыта работы нет.
Может, пойдете к нам с Полиной в книжку соавтором?
Про перевод написал Вам письмо.

Уважаемые, а вам не приходилось работать с модулем PLS для R? Я внимательно изучил мануал, вроде все сделал как нужно, а ничего не получается...
Что самое интересное, я не могу понять, в чем дело...
Не могли бы вы мне дать образец скрипта для работы с PLS?
Или подсказать, где можно найти информацию о работе с этим модулем...

Что _конкретно_ не получается? Какие данные, как обрабатывали? Приведите образец. Могу ответить и на другие вопросы, но сначала давайте разберемся с конкретными Вашими проблемами.

Есть данные о взаимодействие пептидов (9-мерных) с молекулой МНС. Количественной характеристикой степени взаимодействия является pIC50, предикторы--характеристики аминокислот, из которых состоит пептид (каждая аминокислота описывается вектором из шести свойств). Таким образом имеем таблицу:
seq, pIc50, p11, p12... p16... p96

таблица записана в виде файла *.csv
читаем таблицу:
hla<-read.csv("mydata.csv", header=TRUE, rawnames=1)
attach(hla)
формула:
pIC50 ~ p11+p12+...p16+...p96 (в формуле указаны имена столбцов)
mod<-mvr(formula(см. выше), ncomp, data=hla, method = "simpls",
validation = "LOO", model = TRUE)
После попытки записи анализа выдается ошибка, что R не может посчитать среднее значение в столбцах, потому что не все значения являются числами... После проверки и нескольких экспериментов я пришел к выводу, что программе не нравится организация данных: она упорно пытается проводить анализ с использованием названий строк! Может быть нужно преобразовать frame в матрицу?

Пример с конкретными данными смогу написать только завтра, поскольку результаты экспериментов с R остались на другой машине...