Все форумы > Биофизика и матметоды в биологии > версия для печати
Web-адрес темы: http://molbiol.ru/forums/index.php?s=8148751828567d698ac58105ecdcf55e&act=ST&f=43&t=102724

R Help

plantago 26.06.2006 23:06
user posted image

Уважаемый Посетитель,

помните, что все мы безусловно рады общению с вами.

1) У среды статистического анализа R обширная и одновременно точная документация. Дело в том, что R разрабатывали и документировали крайне компетентные и востребованные отраслью специалисты. Если что-то внесено в документацию, то оно в 99% случаев понадобится пользователю.

Поэтому всегда, когда в первый раз приходится воспользоваться какой-либо командой, наберите в интерпретаторе (иногда это отдельное окно среды разработки) сочетание ?имя_команды, и вы получите возможность прочитать справку, куда внесены все нюансы, которые "прочувствовали на себе" предыдущие поколения пользователей.

Документация R писалась добровольцами в течении более 15 лет и не содержит ни одной просто так написанной строки. Очень полезно просто прочитать страницы помощи всех команд R, входящих в пакет {base}, например, тут: http://finzi.psych.upenn.edu/R/library/bas...ml/00Index.html ( http://finzi.psych.upenn.edu/R/library/base/html/00Index.html ) Это немного.

2) Если вы никогда до этого не читали "Введение в R" https://cran.r-project.org/manuals.html ( https://cran.r-project.org/manuals.html ) , то обязательно сделайте это. Дело в том, что логика использования R, заложенная при его создании, подразумевает соблюдение нескольких простых правил. Если пытаться их игнорировать, то ничего, кроме разочарования и постоянного неудобства, испытать не удастся.

3) Если вы ищете какие-то возможности R по анализу ваших данных, полезно начать поиск с страницы https://cran.r-project.org/web/views/ ( https://cran.r-project.org/web/views/ ) , на которой собраны постоянно обновляемые обзоры растущих возможностей R в различных областях прикладного анализа данных.

4) Воспользуйтесь поиском по форуму с помощью гугла:

"ключи поиска" site:http://molbiol.ru/forums/

Если же (1), (2), (3) и (4) не помогли, вы можете написать свой вопрос сюда, в этот раздел форума.

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

И еще -- мы не решаем контрольные работы.

===

Исходная версия первого сообщения:

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



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

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

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



plantago 27.06.2006 23:15
Вместе искать веселее wink.gif
Вот про 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/library/base/html/double.html )
===
Посмотрите еще здесь -- http://finzi.psych.upenn.edu/R/tmp/Rhelp02...hive/71569.html ( http://finzi.psych.upenn.edu/R/tmp/Rhelp02a/archive/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 ( http://cran.arsmachinandi.it/bin/windows/contrib/r-release/data.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 ( http://cran.at.r-project.org/bin/windows/contrib/r-release/g.data_1.6.zip )
===
HTH



Amadeus 30.06.2006 22:40
Идея с 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?

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



Картинки:
Прикреплённое изображение

plantago 01.07.2006 01:56
1. plot(g) wink.gif
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 ( http://finzi.psych.upenn.edu/R/library/multtest/html/mt.rawp2adjp.html ) ). Посмотрите еще thread здесь -- http://finzi.psych.upenn.edu/R/tmp/Rhelp02...hive/65720.html ( http://finzi.psych.upenn.edu/R/tmp/Rhelp02a/archive/65720.html ) И еще вот чего нашел: http://faculty.washington.edu/~jstorey/qvalue ( 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."
Ввожу данные:
LINE = c('A','A' ...

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



Amadeus 01.07.2006 04:53
1. plot(g) wink.gif

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

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)]
    }

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

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

3. Нет, но народ советует использовать пакет multtest (например, http://finzi.psych.upenn.edu/R/library/mul....rawp2adjp.html ( http://finzi.psych.upenn.edu/R/library/multtest/html/mt.rawp2adjp.html )

Вроде как там есть Hochberg, но не уверен что именно тот что есть в SPSS. Тем более что вот в этой ссылке:
). Посмотрите еще thread здесь -- http://finzi.psych.upenn.edu/R/tmp/Rhelp02...hive/65720.html ( http://finzi.psych.upenn.edu/R/tmp/Rhelp02a/archive/65720.html )

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

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

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."

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

(plantago @ 30.06.2006 23:56)
Неужели вот так, "руками", и вводили данные?! Прочитав такое, народ Бог знает что о R может подумать...

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



plantago 01.07.2006 05:46
(Amadeus @ 30.06.2006 18:53)
Ссылка на исходное сообщение  Хм, у меня эта функция выводит четыре графика, и среди них ни одного P-P Plot, есть только Q-Q Plot

Так Вам шашечки или ехать? QQ Plot показывает то же самое, и даже эффективнее: http://finzi.psych.upenn.edu/R/tmp/Rhelp02...hive/57078.html ( http://finzi.psych.upenn.edu/R/tmp/Rhelp02a/archive/57078.html )
Положим я хочу проверить  homogeneity of variances для фактора DRUG и LINE. На levene.test(g,LINE) ругается frown.gif

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

В SPSS-ном виде.
Жалко что этот multtest такое скудное описание имеет

Вот Вам статья про него: http://www.bepress.com/cgi/viewcontent.cgi...text=ucbbiostat ( http://www.bepress.com/cgi/viewcontent.cgi?article=1164&context=ucbbiostat )
Вот именно для таких случаев и пригодился бы Hochberg.

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

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



Amadeus 01.07.2006 15:48
(plantago @ 01.07.2006 03:46)
Ссылка на исходное сообщение  Так Вам шашечки или ехать? QQ Plot показывает то же самое, и даже эффективнее: http://finzi.psych.upenn.edu/R/tmp/Rhelp02...hive/57078.html ( http://finzi.psych.upenn.edu/R/tmp/Rhelp02a/archive/57078.html )

Почитал я про этот 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? Как Вы, например, их делаете?

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

Попробовал - заработало, спасибо. Но что-то он почему то выдает результат абсолютно отличный от 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"?

Вот Вам статья про него: http://www.bepress.com/cgi/viewcontent.cgi...text=ucbbiostat ( http://www.bepress.com/cgi/viewcontent.cgi?article=1164&context=ucbbiostat )

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

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

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

Зато в R есть способы читать из буфера обмена. Под виндами это совсем просто:
> x <- read.delim("clipboard")
Под всеми платформами можно копировать по одной колонке:
> x <- scan()
Затем вставка, затем 2 раза ENTER.
Вот спасибо, это заработало, действительно удобно smile.gif



Картинки:
Прикреплённое изображение Прикреплённое изображение

plantago 03.07.2006 02:54
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 ( http://socserv.mcmaster.ca/jfox ) ). Но боюсь, что без особых резутьтатов -- про отношение к STATISTICA среди статистиков я уже писал, а самое главное -- Вы никогда не узнаете, что действительно использует STATISTICA, тогда как в случае с R это очень просто.
Кстати, если уж про тесты:
> ?bartlett.test
и все тесты, которые там в "See Also".
ANOVA пользуюсь нечасто, штатных средств R мне вполне хватает.



Amadeus 03.07.2006 18:51
(plantago @ 03.07.2006 00:54)
Ссылка на исходное сообщение  1. Про QQ Plot.
Я уже писал, что строить надо командой
> plot(g)
Если Вам нужны не все пять, нужно сделать так:
> plot(g, which=2)
Чтобы получить справку по этой команде, надо сделать так:
> ?plot.lm
потому что R -- по-настоящему объект-ориентированный язык.

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

2. Про R vs. STATISTICA
Можно задать про это вопрос  в R Help, а можно автору кода "Levene Test" -- John Fox ( http://socserv.mcmaster.ca/jfox ( http://socserv.mcmaster.ca/jfox ) ).

Разобрался сам. Дело в том, что 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 smile.gif

Почитал интернет по этому вопросу и оказалось, что второй вариант (реализованный в R, через medians) - точнее (например тут: http://tolstoy.newcastle.edu.au/R/help/03b/1901.html) ( 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... Но что-то мне не хочется демо ставить только ради этого smile.gif

А вообще в ходе поисков обнаружил что в пакете 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 мне вполне хватает.
Последнюю реплику не совсем понял. А разве ANOVA не является штатным средством R? Я думал это один из самых часто используемых параметрических методов? А чем сравниваете несколько групп, если такая нужда есть? Или сразу регрессию проводите? Но там ведь assumptions тоже надо сравнивать...


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



plantago 03.07.2006 22:41
(Amadeus @ 03.07.2006 08:51)
Ссылка на исходное сообщениеТ.е. если все точки (или по крайней мере подавляющее большинство) располагаются на линии или очень близко к ней, то данные можно считать нормальными?

Да
Или наде еще смотреть на угол наклона линии на этом графике

Нет
Разобрался сам.

Отлично! И спасибо за инфу о тестах.
Последнюю реплику не совсем понял.

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

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



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

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

(plantago @ 03.07.2006 20:41)
http://finzi.psych.upenn.edu/R/Rhelp02a/archive/15472.html ( http://finzi.psych.upenn.edu/R/Rhelp02a/archive/15472.html )
Прочитал, спасибо.
Только там форматировать предлагают каждый раз - это немного неудобно для меня. Я бы предпочел глобально это один раз переопределить для всех репортов. Решение такое нашел в доках:
options(scipen=20)
правда длинные они тогда иногда числа получаются, но все равно, для моих eyeballs оно так нагляднее и приятнее smile.gif

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

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 каждый новый график открывать в новом окошке, а не стирать старый?



plantago 04.07.2006 02:10
(Amadeus @ 03.07.2006 15:04)
Ссылка на исходное сообщение А какими тестами/графиками assumptions проверяете?

QQ Plot, Levene test wink.gif
1. Вы строете графики в R или GNUPlot? Какой из них мощнее?

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

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

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



Amadeus 04.07.2006 03:25
(plantago @ 04.07.2006 00:10)
Ссылка на исходное сообщение  QQ Plot, Levene test wink.gif

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


Gnuplot давно забросил. R, конечно же, мощнее. Вот, для забавы, посмотрите график:
user posted image
и код к нему http://www.stat.auckland.ac.nz/~paul/RGrap...examples-once.R ( http://www.stat.auckland.ac.nz/~paul/RGraphics/examples-once.R )

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

В R не очень много интерактивности в графиках (один из примеров ниже), но есть замечательный GGobi ( http://www.ggobi.org ( http://www.ggobi.org ) ), который сию интерактивность добавляет.

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

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

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

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

Перед графиком
> x11()
либо
> windows()

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



plantago 04.07.2006 03:52
(Amadeus @ 03.07.2006 17:25)
Ссылка на исходное сообщение А чем пользуетесь для постройки SVG графиков (если строите)?

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

Дык, это ж Paul Murrell, который с нуля сделал Trellis graphics в R!
А есть к нему какой-нибудь туториал для чайников? Такой который можно безболезненно прожевать?  smile.gif

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

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

?mtext



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



plantago 06.07.2006 02:19
В пакете sem есть bootstrap для structural equation models, что близко к FA. Вы можете выдать свой объект за sem-объект или сделать по аналогии, посмотрев код пакета.
Пакет boot предназначен для применения bootstrap к каким угодно задачам, поэтому Вы можете сами написать (думаю, не очень сложный) код под свою задачу. Вот здесь ( http://cran.at.r-project.org/doc/Rnews/Rnews_2002-3.pdf ( 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 ( http://tolstoy.newcastle.edu.au/R/help/05/02/11520.html )
Объяснено на примере PCA, но, думаю, с FA будет похоже.



Den-N 06.07.2006 19:48
Большое спасибо, попробую разобраться!



Amadeus 08.07.2006 17:09
(plantago @ 04.07.2006 01:52)
Естественно, особенно если Вы предпочитаете автоматизировать! Поменяются точки, и Ваша надпись закроет что-нибудь ценное...
?mtext

Спасибо.

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

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 confused.gif Как его заставить 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?

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



Картинки:
Прикреплённое изображение Прикреплённое изображение

plantago 13.07.2006 03:56
Извините, что отвечаю поздно и мало -- много работы...
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 ( http://finzi.psych.upenn.edu/R/library/R2HTML/DESCRIPTION ) ) Можно еще Rpad ( http://www.rpad.org/Rpad ( http://www.rpad.org/Rpad ) ), но он, по-моему, уж слишком наворочен.



plantago 13.07.2006 06:26
Кстати, не обратили внимание -- столбики на Вашем втором графике кажутся наклоненными друг к другу из-за зрительного эффекта wink.gif



Amadeus 16.07.2006 06:53
(plantago @ 13.07.2006 01:56)
Ссылка на исходное сообщение  Извините, что отвечаю поздно и мало -- много работы...


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

1) >  civ$DRUG <- factor(civ$DRUG, levels=c("X","D","M"))
> interaction.plot(civ$DRUG,civ$LINE,civ$DOSE_ORIG)

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

3) Не понял, объясните подробнее.

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

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

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

6) А зачем? Удобнее?

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

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

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



Картинки:
Прикреплённое изображение

plantago 18.07.2006 02:29
Вот как сделать столбики:
===
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))
===
Я уж Ваши данные в скрипт загнал, чтобы быстрее было smile.gif



Картинки:
Прикреплённое изображение

plantago 18.07.2006 06:14
Кажется, понял Ваш вопрос про разнесение точек!
Это то, что нужно?
===
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



Tolmacheva 18.07.2006 16:21
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.



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



Tolmacheva 19.07.2006 11:31
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' smile.gif Izvinite za bespokoystvo.



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



Amadeus 19.07.2006 17:03
(plantago @ 18.07.2006 00:29)
Ссылка на исходное сообщение  Вот как сделать столбики:
Я уж Ваши данные в скрипт загнал, чтобы быстрее было smile.gif


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()

===

Ляпота! cool.gif

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

Да, Вы правильно уловили что именно мне нужно. Именно функция 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. А эти наши данные хотя и отчитываются по логарифмической шкале, но после преобразования опять представляют собою нормальные концентрации. Как вы считаете?



Картинки (кликните, чтобы посмотреть полноразмерную):
4.7к -

Картинки:
Прикреплённое изображение

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



Amadeus 20.07.2006 16:54
(plantago @ 20.07.2006 06:18)
Ссылка на исходное сообщение  Ого, как у Вас здорово все выходит!
С учетом этого не буду предлагать готового решения. Попробуйте сделать так: выньте из R код interaction.plot() -- для этого достаточно ввечти эту команду без скобок, а потом засуньте внутрь туда куда-нибудь space().

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

Попробовал так как вы сказали:
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 для категориальных данных может? confused.gif

Второй вопрос для меня практически ясен: то, как работать с данными, зависит прежде всего не от того, КАК они были получены, а от того, ЧТО они собою представляют. С этой точки зрения правильно, конечно, вычислять обычные средние.

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



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



Amadeus 21.07.2006 20:24
(plantago @ 21.07.2006 09:48)
Ссылка на исходное сообщение  Я не мучаю, я сам не знаю (честно). Времени просто нет делать длительное исследование вопроса. А Вы все-таки попробуйте написать "interaction.plot" без скобок и нажать Enter...

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

Кстати, может у Вас есть доступ к следующим книжкам:
1. “R Graphics"
by Paul Murrell
http://www.stat.auckland.ac.nz/~paul/RGrap.../rgraphics.html ( 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 ( http://eu.wiley.com/WileyCDA/WileyTitle/productCd-0470022973,subjectCd-ST05,descCd-description.html )

3. Survival Analysis using S,
Tableman, Kim & Portnoy, Chapman & Hall, 2003
http://www.jstatsoft.org/v11/b05/v11b05.pdf ( 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 ( 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 ( 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 в русскоязычной среде. Потому как русководств таких уже достаточно много на английском, но на русском - нет.



plantago 22.07.2006 13:00
(Amadeus @ 21.07.2006 10:24)
Ссылка на исходное сообщение  Попробовал, это мне выдало исходный код функции smile.gif Жалко конечно, но буду копать...

Ну, я Вам сам накопал. См. приложение. Правда, ни space(), ни jitter() не подошли, пришлось писать самому.
2. “Statistics: An Introduction using R”
Michael J. Crawley, March 2005

Есть ксерокс предыдущей его книги про S-PLUS, весьма хорошая вещь. Он эколог, кстати, и большая часть его примеров -- биологические.
5 Introductory Statistics with R
by Peter Dalgaard, Springer, 2002.

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

Книжки очень похожи, так что покупать за деньги не стоит.
Мне кажется что бы действительно пригодилось на данном этапе - простенькая книжка с самыми основами R, с описанием на примерах как провести в R простейший и наиболее часто встречающийся статистический анализ: t-test, ANOVA, непараметрическая статистика, описательные статистики, постройка графиков, анализ выживаемости, etc. Этим можно сразу здорово популяризовать R в русскоязычной среде.

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



Файл/ы:

Скачать файл amadeus.r.zip
Размер:: 1.55к
кол-во скачиваний: 635




Amadeus 22.07.2006 18:36
(plantago @ 22.07.2006 11:00)
Ссылка на исходное сообщение  Ну, я Вам сам накопал. См. приложение. Правда, ни space(), ни jitter() не подошли, пришлось писать самому.


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


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

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


или

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


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

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

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

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

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

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

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



Картинки (кликните, чтобы посмотреть полноразмерную):
6.18к -

Файл/ы:

Скачать файл iplot2.zip
Размер:: 1.94к
кол-во скачиваний: 591




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



Amadeus 05.08.2006 14:41
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 он не работает frown.gif



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



Amadeus 13.08.2006 17:29
(plantago @ 09.08.2006 03:46)
Ссылка на исходное сообщение  UnifyPow принципиально не обновлялся с 1998 года.
Для R есть asypow, pwr, ssanv.

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

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

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



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



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



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



guest: Piter_ 17.11.2006 23:58
http://molbiol.ru/forums/index.php?showtop...960#entry423960 ( http://molbiol.ru/forums/index.php?showtopic=131481&st=0&p=423960#entry423960 )
Dalgaard_P._Introductory_Statistics_With_R_



Guest 18.11.2006 00:16
(Pryanik @ 03.11.2006 10:41)
Ссылка на исходное сообщение  А работал ли кто-нибудь с данными microarray в пакете R?
Я сама только начинаю, то есть опыт практически нулевой, а тратить время зря не хочется...
Может быть кто-то даст мне советы, с чего начать, какие пакеты установить и тд и тп.
Наверное и Bioconductor тоже надо использовать будет, так что любая информация о этой софтине также будет полезна.

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



Piter- 18.11.2006 00:17
хттп://лудщиг-сун2.унил.ч/~дарлене/Рмини/счед.хтмл ( http://ludwig-sun2.unil.ch/~darlene/Rmini/sched.html )
Добавка
есть еше такая штучка :-)
http://bioinf.ucd.ie/people/aedin/R/ ( http://bioinf.ucd.ie/people/aedin/R/ )



gav 15.12.2006 17:00
Очень здравое мнение. Я, собственно, еще в конце позапрошлого года сделал такую попытку, и даже заключил договор с некоторым компьютерным издательством. Но пороху, увы, не хватило, написал только одну главу frown.gif Я теперь всем подряд предлагаю соавторство, вот и Вам тоже wink.gif Может, напишем вместе? Могу прямо тут вывесить план книжки, можно обсудить.


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

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

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

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

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

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


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

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



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



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



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



Yason 05.02.2007 12:49
Есть данные о взаимодействие пептидов (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 остались на другой машине...



plantago 05.02.2007 13:01
Видна, по крайней мере, одна ошибка:
hla<-read.csv("mydata.csv", header=TRUE, rawnames=1)
Параметра "rawnames" у read.csv() нет. Есть row.names.
Вполне серьезно может быть, что причина в этом.
Но лучше всего, конечно, посмотреть на конкретные данные и конкретную реакцию системы.



Yason 06.02.2007 06:33
это я написал неправильно. там было указано row.names
я сегодня еще раз внимательно проверил все численные значения и обнаружил в одной строке (из 2500!) таблицы ошибку, после исправления все заработало!
jump.gif
но не надейтесь wink.gif, если у меня что-то не получится, снова буду мучить вас вопросами!



plantago 06.02.2007 07:28
You are welcome smile.gif



Guest 29.08.2007 20:54
А кто-нибудь когда-нибудь писал функции на R, которые бы задействовали код C/C++? Как это делается? Или где это можно посмотреть? Может у кого примеры найдутся?



plantago 29.08.2007 22:45
Вот здесь все подробно описано:
http://cran.r-project.org/doc/manuals/R-exts.html ( http://cran.r-project.org/doc/manuals/R-exts.html )



Pryanik 29.08.2007 23:09
Спасибо, этот мануал я уже нашла. smile.gif В мануале написано, что для компиляции C кода под R надо использовать функцию .C, но не ясно, где в аргументах прописыват имя файла, в котором находится компилируемая функция. Если пробовать загрузить файл, самый простенкий (с программои "Hello, World!", скажем), то выдает такую вот ерунду

> source("hello.c")
Error in parse(file, n = -1, NULL, "?", srcfile, encoding) :
hello.c: syntax error, unexpected SYMBOL, expecting '\n' or ';' at
2: #include<stdio.h>

Тект проги стандартный

#include<stdio.h>
int main(){
printf("Hello, People!\n");
return 0;
}



plantago 30.08.2007 03:31
0) Я ни разу этого не делал, так что рекомендации приблизительные.
1) Сначала надо скомилировать C в DLL (или SO, если используется *NIX), при помощи R CMD SHLIB в командной строке.
2) Потом надо загрузить получившийся файл функцией dyn.load("hello.dll")
3) Потом надо вызвать его при помощи .C(), .External() или тому подобных функций.



guest: Pryanik 30.08.2007 07:34
Please, подробнее, я чайник.
До этого только пользовалась, но никогда не програмировала ничего сложного. А тут приходится. smile.gif

У меня сейчас Mac, если это важно, уже аж второй ден. smile.gif



plantago 30.08.2007 07:55
Дело в том, что я вот попробовал как раз на Mac'е -- и ничего не вышло, слошные ошибки. А теперь нашел вот такой текст -- http://www.biostat.jhsph.edu/~rpeng/docs/interface.pdf ( http://www.biostat.jhsph.edu/~rpeng/docs/interface.pdf )
Там есть очень простой пример (и не один), попробуйте и напишите, получилось ли.



plantago 30.08.2007 23:57
Ура! У меня получился пример из последнего текста. Так что теперь могу и Вам что-то советовать smile.gif Имейте в виду: у Вас на компьютере (если это Mac) должна быть последняя полная версия R (2.5.1) и последняя версия Xcode (2.4.1), иначе не заработает (проверено).



Guest 01.09.2007 04:31
У меня тоже получились примеры из текста. А что есть Xcode (2.4.1)? Sorry. smile.gif

Я делала все по тексту. Все работает.



plantago 01.09.2007 05:14
Ну, раз они получились, то смотреть, какой у Вас Xcode, не обязательно wink.gif
Вообще говоря, Xcode -- это универсальное средство разработки (в основном на C/C++, но есть и Python), выпускаемое фирмой Apple для Mac (скачивается с их сайта бесплатно и/или устанавливается вместе с дистрибутивом операционки).



Guest 01.09.2007 05:18
Thanks :)



Guest 01.09.2007 05:20
я писала код в текстовом редакторе, сохраняла с расширением .c, а далее по тексту мануала.



Guest 11.09.2007 09:37
Plantago, а вы что-нибудь кроме примеров пытались делать, более или менее слогное?



plantago 12.09.2007 00:26
Нет.



Pryanik 13.09.2007 19:36
Ясно. Просто не понятно, например, как использовать всои функции, а не библиотечные, как узнать как остальные библиотеки (кроме тех, что в примерах) подключать и есче куча вопросов smile.gif

ага. функции можно просто так использовать. хоть отделным фаилом подключать, хоть все в одном описывать.
А как с библиотеками быть, так и не знаю...



plantago 14.09.2007 00:36
Очень трудно отвечать на такие общие вопросы. Что Вы называете "библиотеками"? Сишные библиотеки? Что Вы называете функциями? Функции R?



Pryanik 14.09.2007 02:42
Библиотеки -- сишные. Не знаю... tcl.h, sstream любые, кроме тех что в мануале есть. Как их подключить, чтобы при вызове С кода из R они тоге подключались.
Функции сишные. Например я реализую какую-нибудь функцию, а потом ее исполсую в другои, которую из R вызываю. Я это имела в виду.
Или, на пример, обьекты. Если у меня описан класс, а моя ф-ция, которую я потом из R вызываю, имеет аргумент такого типа. Как мне в .С этот аргумент обозначать?



plantago 14.09.2007 07:10
А просто прописать библиотеки в исходниках, скомпилировать и вызвать из R Вы не можете?



Guest 14.09.2007 07:30
я уге посмотрела исходники какого-то пакета, где есть ф-ции с .С()
В обсчем, зря вас потревогила. Все и так работает.
Просто для math.h специалное имя. Вот я и решила, что для остальных библиотек тоге.



Guest 20.11.2007 12:51
а есть в R что-нибудь похожее на то, что делает PAML?
я имею в виду детекцию сайтов или lineages, находящихся под положительной селекцией.



plantago 21.11.2007 11:21
Попробуйте поискать в Bioconductor или в пакетах, заточенных под разные генетические задачи (их немало на CRAN).



smirnovatatiana 25.11.2007 17:24
скажите пожалуйста, а GUI к нему не появился?? И где искать? А то я еще даже не поняла нужен ли мне именно R, а уже надо учиться им командовать smile.gif)



smirnovatatiana 26.11.2007 00:47
о, есть R Commander, ура! И еще RKWard. осталось поставить и победить установку, которая с первого раза, увы, не хочет... А кто-нибудь здесь их использует?? Я просто в легкой панике: статистика мне не нужна, нужна обработка данных более тривиальная, но тысячи столбцов с цифрами в gnumeric'e/oocalc'e - это похоже на ночной кошмар.
На самом деле мне вот что нужно: есть измерения свечения митохондрий в динамике (куча митохондрий, из разных клеток, с разных опытов), кажется в виде css, или что там AnalyzeParticle в ImageJ на выход дает... Для начала нужно найти те митохондрии, свечение которых давало скачкИ. Все банально, но очень уж их много, потому как из картинки вытаскиваются (почти)автоматически, это я вроде почти победила, ура...
Таблички тоже хочется все скопом, но просто опыта работы с более спец программами чем редактор таблиц - ноль, никогда не сталкивалась. "Статистика", явно, выход, но вся остальная жизнь в линуксе, переезжать неохота frown.gif((
спасибо большое.



plantago 26.11.2007 02:58
С RCommander я работаю в том смысле, что делал его русский перевод. Так что на конкретные вопросы готов отвечать.
Но раз у Вас такие объемы данных (кстати, какие?), то может быть, лучше вот это -- http://rattle.togaware.com/ ( http://rattle.togaware.com/ ) ? К R очень много GUI, вот тут есть список: http://www.sciviews.org/_rgui/ ( http://www.sciviews.org/_rgui/ )
Правильно ли я понимаю, что нахождение скачков -- это фактически поиск локальных максимумов (или минимумов)?



Guest 26.11.2007 11:45
большое спасибо!
а где можно Ваш перевод найти?
Объемы - ожидаю порядка 5-30тыс функций, каждая по 20-100 точкам, разбиты на группы (митохондрии в одной клетке) по 50-100примерно.
Скачки - это разница между соседними значениями, по величине бОльшая, чем порог. Это - отбор претендентов на рассмотрение.
Дальше - проверка всех функций той же группы на наличии каких-либо перепадов, в т.чл. низкоамплитудных, в то же время что и найденный всплеск. После этого хочу посмотреть на графики, а потом уже все остальное
Вообще я правильно поняла что R???



plantago 27.11.2007 07:42
Перевод ставится автоматически, если Вы используете R на русскоязычной локали.
Объемы серьезные, но R должен справиться.
Скачки - это разница между соседними значениями, по величине бОльшая, чем порог.

А, ну это легко. Предположим, у Вас такие числа: 1,5,17,45,97,183,192. Порог -- 25.
Тогда делаем в R:
CODE

x <- c(1,5,17,45,97,183,192)
which(c(0, diff(x)) > 25)
[1] 4 5 6

4,5,6 -- номера значений, превышающие порог.
Вообще я правильно поняла что R???

Вопрос не понял smile.gif



Guest 27.11.2007 13:09
УРА! кажется счастье есть... ох, как же я не люблю думать.. и из ImageJ он грузит таблички - сам и без проблем, прямо через clipboard! Правда ставить пришлось все это на винду... потом буду разбираться с установкой на (ALT)линукс. Почему-то не хочет. А у меня всё там(((
Осталось понять чего собственно надо-то)))))
Огромное спасибо! R рулит! потом еще вернусь, с глупыми вопросами.



Pryanik 07.01.2008 23:28
Люди, кто нибудь имеет опыт делания пакетов для R, с функциями использующими С++ код, который в свою очередь написан с использованием external libraries. Я уже нашла что для этого надо сделать configure скрипт, а чтобы сделать его наде сна4ала сделать configure.ac, а как это сделать я представляю слабо...
У меня в С++ используется gsl.
Могет быть кто-то знает пакеты, в которых в С коде изпользуется эта библиотека?



plantago 08.01.2008 00:08
Почему бы Вам не спросить это в R-Help?



Pryanik 08.01.2008 00:19
Я спросила уже, даже в [R-devel], жду. Просто решила еще и тут спросить, может кто скажет, сейчас еще пойду постер на доску обьявлений повешу smile.gif



Guest 17.01.2008 23:28
Простите за глупый вопрос, а кто-нить знает как в гистограмме покрасить один столбец в другой цвет. Например максимальный, или второй справа, или все что больше определенного значения.



Pryanik 18.01.2008 21:53
> c<-rnorm(20)
> hist(c,col=c("red", "blue","green","red", "white","white"))



Guest 21.01.2008 17:24
(Pryanik @ 18.01.2008 20:53)
Ссылка на исходное сообщение  > c<-rnorm(20)
> hist(c,col=c("red", "blue","green","red", "white","white"))

Спасибо, Полин.



Pryanik 20.04.2008 06:59
А есть ли пакеты для environmental microbiology/metagenomics?
Посмотрела бегло на сайтах R и BioC, но что-то ничего не обнаружила, можек кто пользовал, знает.



guest: Гость 28.05.2008 02:22
Доброго времени суток!
Вопрос по установке RCommander: что необходимо сделать, чтобы при повторной усановке не загружались повторно необходимые пакеты. Дело в том, что при указании установки из имеющихся zip-архивов, все необходимые (и уже скаченные) пакеты не устанавливаются.
R 2.7.0 для win32, R-Commander 1.3-14



plantago 28.05.2008 07:12
Я прочитал Ваше сообщение несколько раз, но так ничего и не понял. Попробуйте, пожалуйста, описать по пунктам, что Вы делаете и что Вы хотите получить.



guest: гость 31.05.2008 20:32
(plantago @ 28.05.2008 06:12)
Ссылка на исходное сообщение  Я прочитал Ваше сообщение несколько раз, но так ничего и не понял. Попробуйте, пожалуйста, описать по пунктам, что Вы делаете и что Вы хотите получить.

Прошу прошения за предыдущее сообщение.

Проблема в следующем:
Имеются zip-арихивы, которые автоматически были скачены при установке RCommander.
При установке R Commander на другой компьютер, указывается расположение пакетов в локальных zip-архивах.
В консоли получаем следущее:

> local({pkg <- select.list(sort(.packages(all.available = TRUE)))
+ if(nchar(pkg)) library(pkg, character.only=TRUE)})
Загрузка требуемого пакета: tcltk
Загружаю интерфейс Tcl/Tk... готово
Error in gzfile(file, "r") : не могу открыть соединение

Версия R Commander 1.3-13


Присоединяю пакет: 'Rcmdr'


The following object(s) are masked from package:tcltk :

ttkentry,
ttkframe,
ttkradiobutton,
ttkscrollbar


Вопрос: что необходимо для устанвки R Commander без подключения к интернету и повторного скачивания необходимых пакетов?
WinXP SP2, R 2.7.0, 1.3-13

Заранее благодарен!



plantago 01.06.2008 02:02
Я попробую повторить Ваши шаги на своих системах, и тогда отпишу.



plantago 05.06.2008 06:52
Попробовал. Понял. Если Вы хотите установить коммандер и все желаемые им пакеты офф-лайн, то их надо сохранить в папку (это Вы делали), а потом поставить через R (а не через R Commander!) командой "Установить из локальных zip-файлов". Можно попробовать установить и через R Commander (как делали Вы), но тогда в папке с ZIP'ами должен быть файл PACKAGES. Как его делать, должно быть написано где-то в документации к коммандеру.



noname00 20.06.2008 01:35
Прекрасная идея с FAQ по R. Прошло много времени, но далее конструкции Вопрос>Ответ нет продвижения. От всей души желаю большей популяризации R!



plantago 20.06.2008 04:38
FAQ -- это и есть "Вопрос - Ответ"...



Guest 20.06.2008 17:44
(plantago @ 20.06.2008 03:38)
Ссылка на исходное сообщение  FAQ -- это и есть "Вопрос - Ответ"...

Спасибо за расшифровку smile.gif Хотелось бы увидеть упорядоченную структуру по темам.



Rjhdby 17.07.2008 16:33
Есть некий массив вот такого вида
CODE

--X1--X2--X3--X4--X5
A  1   1   1   7   1
B 15   0   0   4   0
C  1   7   1  41   5
D  3   0   2   0   3

Допустим "вес" элемента определяется следующим образом
X?=A+B+C+D

Как мне отсортировать этот массив по "весам"?
CODE

--X4--X1--X5--X2--X3
A  7   1   1   1   1
B  4  15   0   0   0
C 41   1   5   7   1
D  0   3   3   0   2




plantago 17.07.2008 22:07
CODE

> a <- data.frame(X1=c(1,15,1,3), X2=c(1,0,7,0), X3=c(1,0,1,2), X4=c(7,4,41,0), X5=c(1,0,5,3))
> row.names(a) <- c("A","B","C","D")
> a[,rev(order(colSums(a)))]
 X4 X1 X5 X2 X3
A  7  1  1  1  1
B  4 15  0  0  0
C 41  1  5  7  1
D  0  3  3  0  2




Rjhdby 18.07.2008 10:03
Дайте я вас расцелую! smile.gif



guest: LuK 31.07.2008 09:27
Всем привет. Собрал немного литературы по языку R, надеюсь поможет в изучении smile.gif

Все книги лежат одним архивом формата 7z на файлообменнике по адресу: http://ifolder.ru/7501417 ( http://ifolder.ru/7501417 )

Список книг:

Baayen R. Analyzing Linguistic Data - A Practical Introduction to Statistics using R.pdf
Chambers J. Software for Data Analysis - Programming with R.pdf
Crawley M. The R Book.pdf
Dalgaard P. Introductory statistics with R.pdf
Everitt B. et al. A Handbook of Statistical Analyses Using R.pdf
Kause A. et al. The Basics of S and S-Plus.pdf
Maindonald J. et al. Data Analysis and Graphics Using R.djvu
Marques J. Applied Statistics Using SPSS, STATISTICA, MATLAB and R.pdf
Murrell P. R Graphics.pdf
Sarkar D. Lattice - Multivariate Data Visualization with R.pdf
Seefeld K. et al. Statistics Using R with Biological Examples.pdf
Spector P. Data Manipulation with R.pdf
Venables W. et al. Modern applied statistics with S.pdf
Venables W. et al. S Programming.pdf
Venables W., Smith D. An Introduction to R.pdf



zollin 22.09.2008 19:26
Здравствуйте,
кто-нибудь знает как в R можно проводить множественные сравнения непараметриченскими методами, т.е. в каких пакетах имеются необходимые функции?

UPD: Уже нашел пакет 'npmc' ( http://finzi.psych.upenn.edu/R/library/npmc/html/00Index.html ). Только мануал там слишком лаконичный. Где можно подробнее узнать о использованных там тестах?



plantago 22.09.2008 20:05
Спросить в R Help, написать автору.



mpyat 10.10.2008 10:54
Коллеги!

Никто не работал в R с пакетами для параллельных вычислений? Вчера обнаружил прелюбопытнейшую (на первый взгляд) статью - http://www.biomedcentral.com/1471-2105/9/390 ( http://www.biomedcentral.com/1471-2105/9/390 ) (вышла 22 сентября).

выжимка абстракта:

Results

We have designed and implemented an R add-on package, R/parallel, that extends R by adding user-friendly parallel computing capabilities. With R/parallel any bioinformatician can now easily automate the parallel execution of loops and benefit from the multicore processor power of today's desktop computers. Using a single and simple function, R/parallel can be integrated directly with other existing R packages. With no need to change the implemented algorithms, the processing time can be approximately reduced N-fold, N being the number of available processor cores.
Conclusion

R/parallel saves bioinformaticians time in their daily tasks of analyzing experimental data. It achieves this objective on two fronts: first, by reducing development time of parallel programs by avoiding reimplementation of existing methods and second, by reducing processing time by speeding up computations on current desktop computers. Future work is focused on extending the envelope of R/parallel by interconnecting and aggregating the power of several computers, both existing office computers and computing clusters.

Конечно, ускориться в два раза почти бесплатно на своей рабочей машине - это очень здорово! Скачал их пакет, пришлось новый R устанавливать (жить под 2.6 он не захотел). Кстати, в CRAN пакета тоже нет, только у них на сайте.

Согласно документации накропал простенькую тестовую функцию:
myfunc <- function( ncycles, nsort)
{
Result=numeric(0);

if( "rparallel" %in% names( getLoadedDLLs()) ) {
runParallel(resultVar="Result", resultOp="rbind" )
#resultVar - which variables within the loop will store the calculation results after each iteration
# resultOp - how these variables have to be 'operated' or 'reduced'.

#As long as each iteration is independent from the others, each iteration will generate
#an independent set of result variables. All these variable sets have to be
#reduced to a single one with the values that it will be obtained in the case
#the loop was run sequentially.
#Examples of reduce operations (i.e. values of resultOp) are: max,'+' or rbind.
}
else
{
# 2. Start of loop
for(i in 1:ncycles) { #Make some calculations
tmparr = rnorm(nsort)
tmpsort = sort(tmparr)
tempResult = tmpsort[1];
Result=rbind(Result,tempResult);
}
}
return( Result )
}

Так вот, параллельная версия действительно загружает оба ядра, но работает на ~10% медленнее (!) чем обычная версия которая грузит только одно ядро. Cтатья "highly aссessed", там говорят что у них все замечательно ускоряется

Что я делаю не так? Надеюсь, тема вызовет интерес smile.gif



PS2004R 10.10.2008 13:16
(mpyat @ 10.10.2008 09:54)
Ссылка на исходное сообщение  Коллеги!

Никто не работал в R с пакетами для параллельных вычислений? Вчера обнаружил прелюбопытнейшую (на первый взгляд) статью - http://www.biomedcentral.com/1471-2105/9/390 ( http://www.biomedcentral.com/1471-2105/9/390 ) (вышла 22 сентября).

выжимка абстракта:

Results

We have designed and implemented an R add-on package, R/parallel, that extends R by adding user-friendly parallel computing capabilities. With R/parallel any bioinformatician can now easily automate the parallel execution of loops and benefit from the multicore processor power of today's desktop computers. Using a single and simple function, R/parallel can be integrated directly with other existing R packages. With no need to change the implemented algorithms, the processing time can be approximately reduced N-fold, N being the number of available processor cores.
Conclusion

R/parallel saves bioinformaticians time in their daily tasks of analyzing experimental data. It achieves this objective on two fronts: first, by reducing development time of parallel programs by avoiding reimplementation of existing methods and second, by reducing processing time by speeding up computations on current desktop computers. Future work is focused on extending the envelope of R/parallel by interconnecting and aggregating the power of several computers, both existing office computers and computing clusters.

Конечно, ускориться в два раза почти бесплатно на своей рабочей машине - это очень здорово! Скачал их пакет, пришлось новый R устанавливать (жить под 2.6 он не захотел). Кстати, в CRAN пакета тоже нет, только у них на сайте.

Согласно документации накропал простенькую тестовую функцию:
myfunc <- function( ncycles, nsort)
  {
    Result=numeric(0);

    if( "rparallel" %in% names( getLoadedDLLs()) ) {
      runParallel(resultVar="Result", resultOp="rbind" )
      #resultVar - which variables within the loop will store the calculation results after each iteration
      # resultOp - how these variables have to be 'operated' or 'reduced'.
     
      #As long as each iteration is independent from the others, each iteration will generate
      #an independent set of result variables. All these variable sets have to be
      #reduced to a single one with the values that it will be obtained in the case
      #the loop was run sequentially.
      #Examples of reduce operations (i.e. values of resultOp) are: max,'+' or rbind.
    }
    else
    {
      # 2. Start of loop
      for(i in 1:ncycles) {    #Make some calculations
        tmparr = rnorm(nsort)
        tmpsort = sort(tmparr)
        tempResult = tmpsort[1];
        Result=rbind(Result,tempResult);
      }
    }
    return( Result )
  }
 
Так вот, параллельная версия действительно загружает оба ядра, но работает на ~10% медленнее (!) чем обычная версия которая грузит только одно ядро. Cтатья "highly aссessed", там говорят что у них все замечательно ускоряется

Что я делаю не так? Надеюсь, тема вызовет интерес smile.gif



тема обсуждалась в R-help, из того что пробовал сам: pnmath собирается и работает на моем дебиан линуксе.

"Juan Pablo Romero Méndez" <jpablo.romero@gmail.com> writes:

> > Hello,
> >
> > The problem I'm working now requires to operate on big matrices.
> >
> > I've noticed that there are some packages that allows to run some
> > commands in parallel. I've tried snow and NetWorkSpaces, without much
> > success (they are far more slower that the normal functions)

Do you mean like this?

> > library(Rmpi)
> > mpi.spawn.Rslaves(nsl=2) # dual core on my laptop
> > m <- matrix(0, 10000, 1000)
> > system.time(x1 <- apply(m, 2, sum), gcFirst=TRUE)
user system elapsed
0.644 0.148 1.017
> > system.time(x2 <- mpi.parApply(m, 2, sum), gcFirst=TRUE)
user system elapsed
5.188 2.844 10.693

? (This is with Rmpi, a third alternative you did not mention;
'elapsed' time seems to be relevant here.)

The basic problem is that the overhead of dividing the matrix up and
communicating between processes outweighs the already-efficient
computation being performed.

One solution is to organize your code into 'coarse' grains, so the FUN
in apply does (considerably) more work.

A second approach is to develop a better algorithm / use an
appropriate R paradigm, e.g.,

> > system.time(x3 <- colSums(m), gcFirst=TRUE)
user system elapsed
0.060 0.000 0.088

(or even faster, x4 <- rep(0, ncol(m)) wink.gif

A third approach, if your calculations make heavy use of linear
algebra, is to build R with a vectorized BLAS library; see the R
Installation and Administration guide.

A fourth possibility is to use Tierney's 'pnmath' library mentioned in
this thread

https://stat.ethz.ch/pipermail/r-help/2007-...ber/148756.html ( https://stat.ethz.ch/pipermail/r-help/2007-December/148756.html )

The README file needs to be consulted for the not-exactly-trivial (on
my system) task of installing the package. Specific functions are
parallelized, provided the length of the calculation makes it seem
worth-while.

> > system.time(exp(m), gcFirst=TRUE)
user system elapsed
0.108 0.000 0.106
> > library(pnmath)
> > system.time(exp(m), gcFirst=TRUE)
user system elapsed
0.096 0.004 0.052

(elapsed time about 2x faster). Both BLAS and pnmath make much better
use of resources, since they do not require multiple R instances.

None of these approaches would make a colSums faster -- the work is
just too small for the overhead.

Martin

> > My problem is very simple, it doesn't require any communication
> > between parallel tasks; only that it divides simetricaly the task
> > between the available cores. Also, I don't want to run the code in a
> > cluster, just my multicore machine (4 cores).
> >
> > What solution would you propose, given your experience?
> >
> > Regards,
> >
> > Juan Pablo
> >
> > ______________________________________________
> > R-help@r-project.org mailing list
> > https://stat.ethz.ch/mailman/listinfo/r-help ( https://stat.ethz.ch/mailman/listinfo/r-help )
> > PLEASE do read the posting guide http://www.R-project.org/posting-guide.html ( http://www.R-project.org/posting-guide.html )
> > and provide commented, minimal, self-contained, reproducible code.

-- Martin Morgan Computational Biology / Fred Hutchinson Cancer Research Center 1100 Fairview Ave. N. PO Box 19024 Seattle, WA 98109 Location: Arnold Building M2 B169 Phone: (206) 667-2793 ______________________________________________ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help ( https://stat.ethz.ch/mailman/listinfo/r-help ) PLEASE do read the posting guide http://www.R-project.org/posting-guide.html ( http://www.R-project.org/posting-guide.html ) and provide commented, minimal, self-contained, reproducible code.

pnmath currently uses up to 8 threads (i.e. 1, 2, 4, or 8).
getNumPnmathThreads() should tell you the maximum number used on your
system, which should be 8 if the number of processors is being
identified correctly. With the size of m this calculation should be
using 8 threads, but the exp calculation is fairly fast, so the
overhead is noticable. On a Linux box with 4 dual-core AMD processors
I get

> m <- matrix(0, 10000, 1000)
> mean(replicate(10, system.time(exp(m), gcFirst=TRUE))["elapsed",])
[1] 0.3859
> library(pnmath)
> mean(replicate(10, system.time(exp(m), gcFirst=TRUE))["elapsed",])
[1] 0.0775

A similar example using qbeta, a slower function, gives

> p <- matrix(0.5,1000,1000)
> setNumPnmathThreads(1)
[1] 1
> mean(replicate(10, system.time(qbeta(p,2,3), gcFirst=TRUE))["elapsed",])
[1] 7.334
> setNumPnmathThreads(8)
[1] 8
> mean(replicate(10, system.time(qbeta(p,2,3), gcFirst=TRUE))["elapsed",])
[1] 0.9576


On an 8-core Intel/OS X box the improvement for exp is much less, but
is similar for qbeta.

luke


On Thu, 10 Jul 2008, Martin Morgan wrote:

> "Juan Pablo Romero Méndez" <jpablo.romero@gmail.com> writes:
>
>> Just out of curiosity, what system do you have?
>>
>> These are the results in my machine:
>>
>>> system.time(exp(m), gcFirst=TRUE)
>> user system elapsed
>> 0.52 0.04 0.56
>>> library(pnmath)
>>> system.time(exp(m), gcFirst=TRUE)
>> user system elapsed
>> 0.660 0.016 0.175
>>
>
> from cat /proc/cpuinfo, the original results were from a 32 bit
> dual-core system
>
> model name : Intel® Core™2 CPU T7600 @ 2.33GHz
>
> Here's a second set of results on a 64-bit system with 16 core (4 core
> on 4 physical processors, I think)
>
>> mean(replicate(10, system.time(exp(m), gcFirst=TRUE))["elapsed",])
> [1] 0.165
>> mean(replicate(10, system.time(exp(m), gcFirst=TRUE))["elapsed",])
> [1] 0.0397
>
> model name : Intel® Xeon® CPU X7350 @ 2.93GHz
>
> One thing is that for me in single-thread mode the faster processor
> actually evaluates slower. This could be because of 64-bit issues,
> other hardware design aspects, the way I've compiled R on the two
> platforms, or other system activities on the larger machine.
>
> A second thing is that it appears that the larger machine only
> accelerates 4-fold, rather than a naive 16-fold; I think this is from
> decisions in the pnmath code about the number of processors to use,
> although I'm not sure.
>
> A final thing is that running intensive tests on my laptop generates
> enough extra heat to increase the fan speed and laptop temperature. I
> sort of wonder whether consumer laptops / desktops are engineered for
> sustained use of their multiple core (although I guess the gaming
> community makes heavy use of multiple cores).
>
> Martin
>
>
>
>> Juan Pablo
>>
>>
>>>
>>>> system.time(exp(m), gcFirst=TRUE)
>>> user system elapsed
>>> 0.108 0.000 0.106
>>>> library(pnmath)
>>>> system.time(exp(m), gcFirst=TRUE)
>>> user system elapsed
>>> 0.096 0.004 0.052
>>>
>>> (elapsed time about 2x faster). Both BLAS and pnmath make much better
>>> use of resources, since they do not require multiple R instances.
>>>
>>
>> ______________________________________________
>> R-help@r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-help ( https://stat.ethz.ch/mailman/listinfo/r-help )
>> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html ( http://www.R-project.org/posting-guide.html )
>> and provide commented, minimal, self-contained, reproducible code.
>
>

--
Luke Tierney
Chair, Statistics and Actuarial Science
Ralph E. Wareham Professor of Mathematical Sciences
University of Iowa Phone: 319-335-3386
Department of Statistics and Fax: 319-335-3017
Actuarial Science
241 Schaeffer Hall email: luke@stat.uiowa.edu
Iowa City, IA 52242 WWW: http://www.stat.uiowa.edu ( http://www.stat.uiowa.edu )



______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help ( https://stat.ethz.ch/mailman/listinfo/r-help )
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html ( http://www.R-project.org/posting-guide.html )
and provide commented, minimal, self-contained, reproducible code.



guest: LuK 11.03.2009 00:01
Всем привет! Я создал во Вконтакте группу, посвящённую R. Предлагаю пользователям соцсети Вконтакте туда подтягиваться =)

Группа называется "Язык программирования R для статистических вычислений (The R Project)" и находится по адресу http://vkontakte.ru/club8142131 ( http://vkontakte.ru/club8142131 )

С уважением,
LuK.



guest: Алёна 30.04.2009 22:52
Вижу, тут собрались очень опытные люди! Подскажите, пожалуйста, реализованы ли в R методы оценки неизвестных параметров распределений? И если да, то как?
Я хочу оценить параметры для гамма-распределения по имеющейся статистике методом максимального правдоподобия.
Спасибо!



plantago 30.04.2009 23:30
http://finzi.psych.upenn.edu/R/library/STA...l/gammaMLE.html ( http://finzi.psych.upenn.edu/R/library/STAR/html/gammaMLE.html )
Подходит?



tensor 09.08.2009 17:59
Уважаемые коллеги, подскажите есть ли возможность или какой-либо модуль для расчета U Манна-Уитни. На иностранных сайтах все отправляют на работу с модулем вилкоксон.тест, но ведь тест вилкоксона насколько я помню предназначен для работы с взаимосвязанными группами, в то время как М-У для независимых. При этом никакой возможности вывести значение U нету. Мы работаем с достаточно малыми выборками распредление которых отлично от нормального поэтому использование столь жесткого критерия как Вилкоксона местами неприемлемо.
Подскажите пожалуйста как произвести расчет U Манна-Уитни в R.

И второй вопрос - насколько я понял при импорте таблиц из других форматов в источнике не должно быть пустых ячеек (выпадающих значений). Что вставить на их место для корректного импорта в R?
Копирование каждой переменной (варианты) по отдельности - неприемлемо, потому что порой количество их выходит за 300 - например в офтальмологических исследованиях, иммунном фенотипировании, онкологии. Есть ли возможность автоматической замены значений пустых ячеек на NA при импорте данных?



-daf- 11.08.2009 14:58
Для взаимосвязанных групп (т.е x[i] соответствует y[i]) применяется 'Wilcoxon signed-rank test', а для независимых 'Mann–Whitney U test' или как его еще называют 'Wilcoxon rank-sum test'. Вам я так понял нужно второе. Оба этих теста можно сделать с помощью функции wilcox.test из стандартного пакета 'stats'. Для 'Wilcoxon rank-sum test' используйте wilcox.test(x,y), но в качестве статистики он выдаст W, а не U. К U легко перейти так как
U=W + m*(m+1)/2, где m - размер x.
Второй вопрос: а каким способом вы импортируете данные? функция read.table с этим легко справляется



voliadis 11.08.2009 20:56
Немного об U-статистике. Как только ее не пробовали вычислять (так ( http://www.nabble.com/mann-whitney-U-test---U-value-td24442713.html#a24443086 ), так ( http://www.nabble.com/Mann-Whitney-U-td12151974.html ) и так ( http://www.nabble.com/-R--The-W-statistic-in-wilcox.exact-td6664544.html )).

Возмем пример из http://www.nabble.com/Mann-Whitney-U-td12151974.html ( http://www.nabble.com/Mann-Whitney-U-td12151974.html )

group1<-c(1.34,1.47,1.48,1.49,1.62,1.67,1.7,1.7,1.7,1.73,1.81,1.84,1.9,1.96,2,2,2.19,2.29,2.29,2.41,2.41,2.46,2.5,2.6,2.8,2.8,3.07,3.3)
group2<-c(0.98,1.18,1.25,1.33,1.38,1.4,1.49,1.57,1.72,1.75,1.8,1.82,1.86,1.9,1.97,2.04,2.14,2.18,2.49,2.5,2.55,2.57,2.64,2.73,2.77,2.9,2.94)

Если сделать так: wilcox.test(group1, group2), то W = 405.5. Однако, SPSS и Statistica дают результат U = 350.5. Формула, указанная выше, здесь не сработает.

Но если в R написать wilcox.test(group2, group1), то получим W = 350.5. Мы получили, что W=U.

В любом случае p-value = 0.6494 ( R ) и 0,643 (SPSS, Statistica).

Вывод: Какую бы величину U или W не давали нам программы, величина p одинакова. И это хорошо umnik.gif



tensor 16.08.2009 14:36
На одном из форумов (правда не помню на каком) я кстати видел объяснение тому феномену, что различается значение р. В статистике большиство методов сходных с дисперсионным анализом идут классическим расчетом через среднее арифметическое, в R же используются в ряде случаев более новые методы доработанные для использования медианы.



-daf- 18.08.2009 04:04
По поводу того, что формула не верна - это не так. Просто существует путаница в определении статистик W и U. Обычно W определяют как количество пар (x[i],y[j]) таких, что x[i]>y[j], а если существуют такие (i,j), что x[i]=y[j], то еще добавляют 1/2 от количества пар таких что x[i]=y[j]. То есть
W <- function(x,y) { sum(sapply(x,function(xi) sum(y<xi)+sum(y==xi)/2)) }
Именно эта статистика и вычисляется в 'wilcox.test'. U-статистика обычно определяется как сумма ранков элементов x в объедененном векторе т.е.(в c(x,y)). Ранк обычно определяется также как в функции 'rank', т.е. для одинаковых значение вычисляется средний ранк.
U <- function(x,y){sum(rank(c(x,y))[seq_along(x)])}
Путем не хитрых манипуляций можно понять что эти две статистики связаны соотношением
U(x,y) = W(x,y) + length(x)*(length(x)+1)/2
Вот например:
> W <- function(x,y) { sum(sapply(x,function(xi) sum(y<xi)+sum(y==xi)/2)) }
> U <- function(x,y){sum(rank(c(x,y))[seq_along(x)])}
> group1<-c(1.34,1.47,1.48,1.49,1.62,1.67,1.7,1.7,1.7,1.73,1.81,1.84,1.9,1.96,2,2,2.19,2.29,2.29,2.41,2.41,2.46,2.5,2.6,2.8,2.8,3.07,3.3)
> group2<-c(0.98,1.18,1.25,1.33,1.38,1.4,1.49,1.57,1.72,1.75,1.8,1.82,1.86,1.9,1.97,2.04,2.14,2.18,2.49,2.5,2.55,2.57,2.64,2.73,2.77,2.9,2.94)
> wilcox.test(group1,group2)$statistic
W
405.5
> W(group1,group2)
[1] 405.5
> U(group1,group2)
[1] 811.5
> W(group1,group2)+length(group1)*(length(group1)+1)/2
[1] 811.5

кстати есть еще пакет 'coin':
> library(coin)
> data <- c(group1, group2)
> groups <- factor(c(rep('group1', length(group1)),rep('group2', length(group2))))
> statistic(wilcox_test(data~groups),type='linear')
group1 811.5

он как раз вычисляет U-статистику.
SPSS и Statistica по-видимому выдают другие статистики, я с ними не работал,
но похоже что они выдают W(y,x) (все равно что в определении знак > заменить на <).
Значения p конечно будут везде похожими, но точное значение p вычислить сложно,
зато статистика вычисляется точно и если не запутаться в определениях то пользоваться ей гораздо проще.
Вот например я получил несколько значений p для нашего примера:
> wilcox.test(group1, group2, exact=F)$p.value
[1] 0.6493547
> wilcox.exact(group1, group2, exact=T)$p.value
[1] 0.6458658
> pvalue(wilcox_test(data~groups, distribution='exact'))
[1] 0.6488684
> pvalue(wilcox_test(data~groups, distribution='approximate'))
[1] 0.626
> pvalue(wilcox_test(data~groups, distribution='asymptotic'))
[1] 0.6433081



dlinnosheee 12.09.2009 21:46
Вопрос профессионалам R:
Есть ли в R специальные механизмы для создания онлайн-приложений?
Типа чтобы пользователь на сайте ввел данные, они обработались на сервере и на выходе результаты (всякие графики) опять на сайте. Если такая фича есть, то насколько просто ей пользоваться (по сравнению с использованием для вычислений C++, Fortran или JAVA).



plantago 12.09.2009 22:11
http://cran.r-project.org/doc/FAQ/R-FAQ.ht...-Web-Interfaces ( http://cran.r-project.org/doc/FAQ/R-FAQ.html#R-Web-Interfaces )



tensor 13.09.2009 06:47
Такая возможность непременно есть. Так как R поддерживает импорт данных из БД MySQL, PSQL, ODBC. А так же существуют модули, позволяющие включать скрипты R в языки программирования. Если пришлете в личку мейл - отправлю вам несколько книжек, где описано как это делается.



dlinnosheee 13.09.2009 13:57
я немного вник в R, умею уже делать в нем элементарные вещи. По ходу я продолжаю отсматривать другие скриптовые языки. Возникли вопросы:

1) Что проще в освоении и использовании, R или SAGE на базе Python?
http://en.wikipedia.org/wiki/Software_for_...Experimentation ( http://en.wikipedia.org/wiki/Software_for_Algebra_and_Geometry_Experimentation )

2) Есть ли в R или в пакетах к нему встроенный оператор аналитического дифференцирования (как в Matematica, чтобы результат был не цифра, а математическое выражение).

3) Как в R хранится код созданный пользователем (в каких файлах, куда он девается восле выполнения команд во встроенном базовом GUI)?



PS2004R 14.09.2009 07:32
(dlinnosheee @ 13.09.2009 13:57)
Ссылка на исходное сообщение  я немного вник в R, умею уже делать в нем элементарные вещи. По ходу я продолжаю отсматривать другие скриптовые языки. Возникли вопросы:

1) Что проще в освоении и использовании, R или SAGE на базе Python?
http://en.wikipedia.org/wiki/Software_for_...Experimentation ( http://en.wikipedia.org/wiki/Software_for_Algebra_and_Geometry_Experimentation )

2) Есть ли в R или в пакетах к нему встроенный оператор аналитического дифференцирования (как в Matematica, чтобы результат был не цифра, а математическое выражение).

3) Как в R хранится код созданный пользователем (в каких файлах, куда он девается восле выполнения команд во встроенном базовом GUI)?

1) Как то очень флеймообразующе звучит... Для ответа на вопрос надо указывать что конкретно (с приведением кода) надо решить.
2) http://www.google.com/search?q=maxima&doma...G=Google+Search ( http://www.google.com/search?q=maxima&domains=r-project.org&sitesearch=r-project.org&btnG=Google+Search )
http://www.google.com/search?hl=ru&safe=of...h=r-project.org ( http://www.google.com/search?hl=ru&safe=off&domains=r-project.org&sitesearch=r-project.org&q=symbolic+maxima&btnG=%D0%9F%D0%BE%D0%B8%D1%81%D0%BA&sitesearch=r-project.org )
3) Все объявленные функции хранятся между сеансами в дампе памяти. Очень удобно и классически для сред анализа.
Второй подход --- писать их в текстовом файле в виде программы. Облегчает запуск на других платформах, решая вопросы совместимости.



PS2004R 14.09.2009 10:28
(dlinnosheee @ 13.09.2009 13:57)
Ссылка на исходное сообщение  я немного вник в R, умею уже делать в нем элементарные вещи. По ходу я продолжаю отсматривать другие скриптовые языки. Возникли вопросы:

1) Что проще в освоении и использовании, R или SAGE на базе Python?
http://en.wikipedia.org/wiki/Software_for_...Experimentation ( http://en.wikipedia.org/wiki/Software_for_Algebra_and_Geometry_Experimentation )

Кстати, если очень привычка к питону, то можно работать с R из питона.

python-rpy - Python interface to the GNU R language and environment
python-rpy-doc - Python interface to the GNU R language (documentation package)



dlinnosheee 15.09.2009 10:01
(PS2004R @ 14.09.2009 08:32)
Ссылка на исходное сообщение  
2)
http://www.google.com/search?q=maxima&doma...G=Google+Search ( http://www.google.com/search?q=maxima&domains=r-project.org&sitesearch=r-project.org&btnG=Google+Search )

http://www.google.com/search?hl=ru&safe=of...h=r-project.org ( http://www.google.com/search?hl=ru&safe=off&domains=r-project.org&sitesearch=r-project.org&q=symbolic+maxima&btnG=%D0%9F%D0%BE%D0%B8%D1%81%D0%BA&sitesearch=r-project.org )


Так maxima это не R. Это скорее ближе к SAGE. И как ее к R ее подключить - это надо разбираться еще...



PS2004R 15.09.2009 11:46
(dlinnosheee @ 15.09.2009 10:01)
Ссылка на исходное сообщение  Так maxima это не R. Это скорее ближе к SAGE. И как ее к R ее подключить - это надо разбираться еще...

Первые линки в поиске соответственно что может делать R,
http://wiki.r-project.org/rwiki/doku.php?i...slations:maxima ( http://wiki.r-project.org/rwiki/doku.php?id=getting-started:translations:maxima )
http://wiki.r-project.org/rwiki/doku.php?i...s:mathematica2r ( http://wiki.r-project.org/rwiki/doku.php?id=getting-started:translations:mathematica2r )
Так что надо посчитать то в R? Напоминаю, этот форум называется "R Help". Про возможность работать с R из питона написано. Sage обсуждать здесь...
Загадка: зеленое, сидит под мостом, на т называется.



dlinnosheee 15.09.2009 22:10
(PS2004R @ 15.09.2009 12:46)
Ссылка на исходное сообщение 
Так что надо посчитать то в R? Напоминаю, этот форум называется "R Help". Про возможность работать с R из питона написано. Sage обсуждать здесь...


Моя задача формулируется не "посчитать", а создать онлайн-приложение. Математика - динамическое программирование, расчеты длительные, объем памяти большой. Считаться думаю будет на сервере, а клиент вводит только свои данные. На данный момент часть ядра которая завязана на вычисления работает на фортране. Скриптовых языков типа R или Python я не знаю, но готов выучить и переписать на нем (это для тренировки тоже). Решил остановиться на R потому что (1) есть вроде возмодность создавать веб-интерфейсы и (2) популярный язык вроде становится, в жизни пригодяится. С практической точки зрения, меня сейчас интересует как ввод-вывод онлайн организовать.



PS2004R 16.09.2009 07:28
(dlinnosheee @ 15.09.2009 22:10)
Ссылка на исходное сообщение  Моя задача формулируется не "посчитать", а создать онлайн-приложение. Математика - динамическое программирование, расчеты длительные, объем памяти большой. Считаться думаю будет на сервере, а клиент вводит только свои данные. На данный момент часть ядра которая завязана на вычисления работает на фортране. Скриптовых языков типа R или Python я не знаю, но готов выучить и переписать на нем (это для тренировки тоже). Решил остановиться на R потому что (1) есть вроде возмодность создавать веб-интерфейсы и (2) популярный язык вроде становится, в жизни пригодяится. С практической точки зрения, меня сейчас интересует как ввод-вывод онлайн организовать.


CRAN Task View: Optimization and Mathematical Programming

Ну и техника написания интерфейса к существующему фортран коду.



PS2004R 29.09.2009 11:52
(dlinnosheee @ 13.09.2009 13:57)
Ссылка на исходное сообщение  
2) Есть ли в R или в пакетах к нему встроенный оператор аналитического дифференцирования (как в Matematica, чтобы результат был не цифра, а математическое выражение).



http://cran.r-project.org/web/packages/Ryacas/index.html ( http://cran.r-project.org/web/packages/Ryacas/index.html )

вот попалось в списке рассылке



marinakom 08.11.2009 14:34
Коллеги!
Для построения множественной линейной регрессии пошаговым способом (включением, исключением или stepwise) рекомендуют функцию stepAIC( ) из пакета MASS.
Вопросы:
1) Что есть AIC? как оно расшифровывается?
2) У меня сложилось ощущение, что метод поиска лучших предикторов в R чуток иной, нежели F-включения или исключения, реализованный в других статистических пакетах, да или нет?
3) Получается, что данная функция только помогает выбрать предикторы, а чтоб проверить модель на их основе, надо её просто построить и до ума довести руками, так? Не задавала бы этого вопроса, если б незначимые предикторы не оставались в подобной модели (при построении модели методом исключения -- backward).
4) Может, ещё какой путь есть для пошаговой регрессии?
5) Как стандартизованные (beta) коэффициенты регрессии найти (кроме как ручками стандартизовать каждую переменную и перестраивать модель)?



PS2004R 08.11.2009 16:39
(marinakom @ 08.11.2009 13:34)
Ссылка на исходное сообщение  Коллеги!
Для построения множественной линейной регрессии пошаговым способом (включением, исключением или stepwise) рекомендуют  функцию stepAIC( ) из пакета MASS.
  Вопросы:
1) Что есть AIC? как оно расшифровывается?


http://en.wikipedia.org/wiki/Akaike_information_criterion ( http://en.wikipedia.org/wiki/Akaike_information_criterion )



ИгорьЧерниенко 08.02.2010 09:58
Рад приветствовать!
Такая проблема: в виндоуз при попытке подключить RPostgreSQL требует какую-то длльку, при скачивании оной и установки вручную требует уже другую длльку. Не помню названий, если нужно, вышлю подробности. Причем в версии 2.8 такой проблемы нет, но наблюдается во всех более поздних. В чем может быть дело?



plantago 08.02.2010 10:13
Напишите, какие DLL, посмотрим.



Pryanik 09.02.2010 03:42
Несколько вопросов:
1. Как в R сделать 2 графика с одним заголовком?
2. Как в R вставить один график в другой?



Guest 09.02.2010 04:05
Хотел подключить пакет чтобы написать, что за DLL, а он заработал. Дело было, видимо, в системе. Но все равно, вопрос пока оставляю открытым :о)



ИгорьЧерниенко 09.02.2010 04:06
Прошу прощения, забыл подписаться



PS2004R 09.02.2010 15:37
(Guest @ 09.02.2010 04:05)
Ссылка на исходное сообщение  Хотел подключить пакет чтобы написать, что за DLL, а он заработал. Дело было, видимо, в системе. Но все равно, вопрос пока оставляю открытым :о)


значит аналог ldconfig линуксового, для винды это перегрузка системы smile.gif



PS2004R 09.02.2010 15:43
(Pryanik @ 09.02.2010 03:42)
Ссылка на исходное сообщение  Несколько вопросов:
1.    Как в R сделать 2 графика с одним заголовком?
2.    Как в R вставить один график в другой?


?layout для нескольких графиков, не совсем понял про заголовки... Обычно заголовок для сложной иллюстрации генерируется уже в текстовом процессоре. Например в LaTeX.

На либрусеке есть книжки по R, в том числе по графике.



plantago 09.02.2010 16:33
Либрусек вроде бы недоступен уже почти неделю, но книжки по R лежат на либгене: http://gen.lib.rus.ec/ ( http://gen.lib.rus.ec/ )
===
Можно еще не обновлять дисплей (опция add=TRUE).



PS2004R 10.02.2010 10:43
http://lib.rus.ec/b/156264 ( http://lib.rus.ec/b/156264 ) это pdf с книгой Murrell P. R Graphics. Если не отвечает, скажите куда её выслать.



bubnilkin 16.03.2010 10:34
Скажите, пожалуйста, есть ли возможность представить в R столбиковую диаграмму с интервалами Габриэля?

http://udel.edu/~mcdonald/statanovaunplanned.html ( http://udel.edu/~mcdonald/statanovaunplanned.html )

Gabriel, K.R. 1978. A simple method of multiple comparison of means. J. Amer. Stat. Assoc. 73: 724-729.

А как упомянуть в статье на то что статистика делалась в R ?



plantago 16.03.2010 11:01
Если сами интервалы известны, то нарисовать такую диаграмму -- не проблема (см. примеры выше по треду). Если же их надо рассчитать, то похоже, что такого в R нет. Вместо этого рекомендуется использовать более продвинутые методы множественных сравнений, например, из пакета multcomp
===
Упомянуть -- что-то типа
All statistical calculations and graphs were made in the R environment and language (R Development Core Team, 2007).
...
R Development Core Team. (2007) R: a language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. URL: http://www.R-project.org ( http://www.R-project.org ).



tensor 16.03.2010 11:01
А как упомянуть в статье на то что статистика делалась в R ?
========
Я обычно получаю данные с использованием citation() и citation(название модуля) и привожу ссылки на разработчика в виде (Автор, год).



bubnilkin 16.03.2010 14:33
> plantago, tensor:

Спасибо smile.gif



guest: Иван 21.03.2010 22:10
Товарищи, здраствуйте!

Хочу попросить помощи по такому вопросу: как в R можно в виде таблицы вывести результаты вычислений? Примерно как это делается функцией summary().
Только начал знакомиться с R, буду рад любой помощи!



plantago 22.03.2010 01:30
http://finzi.psych.upenn.edu/R/library/R2HTML/html/HTML.html ( http://finzi.psych.upenn.edu/R/library/R2HTML/html/HTML.html )
http://finzi.psych.upenn.edu/R/library/Hmisc/html/latex.html ( http://finzi.psych.upenn.edu/R/library/Hmisc/html/latex.html )
http://finzi.psych.upenn.edu/R/library/odf...l/odfTable.html ( http://finzi.psych.upenn.edu/R/library/odfWeave/html/odfTable.html )
http://cran.at.r-project.org/web/packages/r2lh/index.html ( http://cran.at.r-project.org/web/packages/r2lh/index.html )
http://cran.at.r-project.org/web/packages/...ools/index.html ( http://cran.at.r-project.org/web/packages/reporttools/index.html )



bubnilkin 22.03.2010 12:56
добрый день:)

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


Возможно вопрос глупый, т.к. с цитометрией и с R "на Вы"



PS2004R 22.03.2010 15:17
(bubnilkin @ 22.03.2010 12:56)
Ссылка на исходное сообщение  добрый день:)

Суть: создать (полу)автоматический анализ данных экспериментов по проточной цитометрии. Сами пакеты, к удивлению, давно уже оказывается есть. Но может уже есть готовые интерфейсы (если это не противоречит философии R)? Не хотелось бы изобретать велосипед. Т.е.: указал путь к папке (а там сотни файликов...) нажал кнопочку и пошло экспортироваться и т.п.
Возможно вопрос глупый, т.к. с цитометрией и с R "на Вы"


Есть способ написать TCL/TK интерфейс, если я правильно понял, и внутри вызывать функции пакетов. Это если я правильно понял. smile.gif

Лучше напишите какие пакеты Вы планируете использовать для цитометрии.



bubnilkin 22.03.2010 15:56
(PS2004R @ 22.03.2010 15:17)
Ссылка на исходное сообщение  Есть способ написать TCL/TK интерфейс, если я правильно понял, и внутри вызывать функции пакетов. Это если я правильно понял. smile.gif

Лучше напишите какие пакеты Вы планируете использовать для цитометрии.


Я думал что-нить попроще из этого списка выбрать:

http://bioconductor.org/docs/workflows/flowoverview/ ( http://bioconductor.org/docs/workflows/flowoverview/ )

но пока с функционалом не разобрался

"внутри вызывать функции пакетов" хм... это = "отмечать галочками тот/иной тест"? если да, то примерно так smile.gif



PS2004R 22.03.2010 16:02
(bubnilkin @ 22.03.2010 15:56)
Ссылка на исходное сообщение  Я думал что-нить попроще из этого списка выбрать:

http://bioconductor.org/docs/workflows/flowoverview/ ( http://bioconductor.org/docs/workflows/flowoverview/ )

но пока с функционалом не разобрался

"внутри вызывать функции пакетов" хм... это = "отмечать галочками тот/иной тест"? если да, то примерно так smile.gif


http://bioconductor.org/packages/release/b...To-flowCore.pdf ( http://bioconductor.org/packages/release/bioc/vignettes/flowCore/inst/doc/HowTo-flowCore.pdf )

вот это Вы читали уже?

PS собственно стандартный графический интерфейс у него уже есть, там любому приросшему к мышке просто раздолье smile.gif http://bioconductor.org/docs/workflows/flo...y/tutorial.mpeg ( http://bioconductor.org/docs/workflows/flowcytometry/tutorial.mpeg ) , если сложится быстрый набор команд с небольшим выбором, то можно легко вставить вставку графического интерфеса в нужный момент.



guest: bubnilkin 22.03.2010 16:18
(PS2004R @ 22.03.2010 16:02)
Ссылка на исходное сообщение  http://bioconductor.org/packages/release/b...To-flowCore.pdf ( http://bioconductor.org/packages/release/bioc/vignettes/flowCore/inst/doc/HowTo-flowCore.pdf )

вот это Вы читали уже?


спасибо что ткнули (я оказывается ужо скачал, но пока не прочитал....) сижу с основами -- изучаю такие понятия как VECTOR, MATRIX, DATAFRAME. Дорогу осилит идущий



Guest 22.03.2010 16:21
smile.gif [url=http://bioconductor.org/docs/workflows/flowcytometry/tutorial.mpeg-->
(PS2004R @ 22.03.2010
PS собственно стандартный графический интерфейс у него уже есть, там любому приросшему к мышке просто раздолье smile.gif [url=http://bioconductor.org/docs/workflows/flowcytometry/tutorial.mpeg)
http://bioconductor.org/docs/workflows/flo...y/tutorial.mpeg[/url] ,  если сложится быстрый набор команд с небольшим выбором, то можно легко вставить вставку графического интерфеса в нужный момент.


фильмец просмотрел.... но (не в обиду создавшему), какойто он совсем уж полуавтоматический....

а сам GUI можно создавать средствами R? confused.gif



PS2004R 22.03.2010 16:38
(Guest @ 22.03.2010 16:21)

фильмец просмотрел.... но (не в обиду создавшему), какойто он совсем уж полуавтоматический....

а сам GUI можно создавать средствами R? confused.gif


Да, легко... Например rpanel пакет, есть и gtk, и сразу платформонезависимо можно. Хоть сразу на web и на десктоп распространяй результат. Статья в Линуксформате есть в этом году, по моему второй номер.



ИгорьЧерниенко 01.04.2010 07:14
Рад приветствовать!
Существует ли пакет для разделения смеси нормальных распределений? Порылся на кране - так и не разобрался.



plantago 01.04.2010 07:20
fastICA?



mpyat 01.04.2010 11:09
mclust?



-daf- 01.04.2010 14:07
(ИгорьЧерниенко @ 01.04.2010 08:14)
Ссылка на исходное сообщение  Рад приветствовать!
Существует ли пакет для разделения смеси нормальных распределений? Порылся на кране - так и не разобрался.

mixtools



ИгорьЧерниенко 03.04.2010 15:58
Большое спасибо! Как раз то, что надо.



Alexander bk 04.04.2010 21:30
Для тех, кто переходит с Виндовс на открытое ПО и сталкивается с проблемой поиска аналогов MS Office, Statistica и MathCad неплохо было бы какое-нибудь начальное руководство, чтобы человек смог установить и увидеть все прелести этой программы.
Я понимаю, что есть фанаты консоли, но мне бы например хотелось бы графический интерфейс. Если такое руководство есть, то поделитесь ссылкой пожалуйста. Если нет, то подскажите, от куда начинать.

Спасибо.



PS2004R 04.04.2010 23:48
(Alexander bk @ 04.04.2010 21:30)
Ссылка на исходное сообщение неплохо было бы какое-нибудь начальное руководство, чтобы человек смог установить и увидеть все прелести этой программы.



По R просто масса литературы. Очень правда странно что Вы начинаете вопрос с каких то условий.

Графический интерфейс есть Rcmdr, но его можно и самому писать.

Обычно используют связку Sweave - LaTeX , она наиболее продуктивна в работе. Очень помогает при этом работа в связке emacs-ess-noweb-auctex-reftex, что позволяет работать с реально сложными проектами не теряя нити.

Сайт программы http://cran.r-project.org/ ( http://cran.r-project.org/ )

Найти свою область можно в разделе http://cran.r-project.org/web/views/ ( http://cran.r-project.org/web/views/ ) .

Введение в сам пакет http://cran.r-project.org/doc/manuals/R-intro.html ( http://cran.r-project.org/doc/manuals/R-intro.html ) перевод на русский http://m7876.wiki.zoho.com/Introduction-to-R.html ( http://m7876.wiki.zoho.com/Introduction-to-R.html )

Масса книг есть на либрусеке.



plantago 05.04.2010 04:33
Мы пишем книгу по R на русском. Сделана на 70%, но идет медленно.



Alexander bk 05.04.2010 22:38
Обычно используют связку Sweave - LaTeX , она наиболее продуктивна в работе. Очень помогает при этом работа в связке emacs-ess-noweb-auctex-reftex, что позволяет работать с реально сложными проектами не теряя нити.


Извините, но я боюсь, что эти слова могут быть не совсем понятны некоторым начинающим пользователям. Да, вы можете отправить искать в Гугл, заранее спасибо. Хотим мы этого или нет, но единственное преимущество Виндовс и приложений под него в их простоте и доступности пользователю и это преимущество решающее при выборе.

Допустим, что я как обычный пользователь, только что перешедший на Линукс и потративший три недели просто на то, чтобы настроить звук в своей ОС. Но когда я начинаю применять это для своей работы, то сталкиваюсь с еще большей кучей трудностей. В итоге, ввиду отсутствия времени, я махаю на все это рукой и возвращаюсь к МС Эксель, Статистике и Маткаду. Это как пример.

Я пишу это не потому, что я такой упертый. Просто уже достало, что везде нас обязывают, если это презентация, то МС ПауэрПоинт, если графики, то МС Эксель, если статобработка, то Статистика. Доходит до анекдота, когда приносишь отчет, сделанный в другой среде и его отклоняют просто потому, что это не Майкрософт. Хотим мы или нет, но свободное ПО надо распространять и надо делать его как можно более доступным. Иначе кроме как нас самих оно никому не будет нужно.

Я к тому, что давайте создавать руководства и ветки на форуме, но хорошо бы, чтобы это было как минимум две отдельных темы: одна для новичков и отдельно для тех, кто уже знаком и долгие годы занимается этим.

2Plantago, спасибо большое, надеюсь, мы все скоро увидим эту книгу, вы вообще великий человек.



plantago 05.04.2010 23:53
Вот есть такой "простецкий" текст. Пригодится?



Файл/ы:

rabotaem_v_r.pdf
Размер:: 150.51к
кол-во скачиваний: 6




PS2004R 06.04.2010 07:30
(Alexander bk @ 05.04.2010 22:38)
Ссылка на исходное сообщение  Извините, но я боюсь, что эти слова могут быть не совсем понятны некоторым начинающим пользователям. Да, вы можете отправить искать в Гугл, заранее спасибо. Хотим мы этого или нет, но единственное преимущество Виндовс и приложений под него в их простоте и доступности пользователю и это преимущество решающее при выборе.



Я ничего не понимаю frown.gif
При чём тут виндовс и линукс? При чём коммерческий софт и свободный? Всё что я написал работает под Виндовс. Вы даже можете купить коммерческий emacs и R если это необходимо.

Откуда какие то предположения про гугл? Отправил я Вас прочитать официальный сайт R (конкретный раздел про области применения конкретных средств R), и русскоязычный перевод руководства для начинающих. Вы хоть ссылки открывали?

Теперь попробую погадать я.

Вы что ? Одновременно ведёте дискуссию ещё с кем то о пользе Виндовс и ответили не туда? smile.gif
Вы потролить пришли?



Pryanik 06.04.2010 08:25
Товарищи, как сделать в [R barplot] покрашеный разными цветами в разных комбинациях?
Была уверена, что вот это

barplot(rbind(c(1,2,1,2),c(1,1,1,1)),col=rbind(c("red","red","red","green"),c("green","green","green","red")))


даст 4 столбика. У 1-3 низ будет красный а верх зеленый, а у 4 наоборот. Но нет, выдает все одинаковые столбы, низ красный, верх зеленый. Как это исправить?

Заранее спасибо!



Файл/ы:

1.pdf
Размер:: 2.35к
кол-во скачиваний: 3




plantago 06.04.2010 08:52
Получите smile.gif
barplot(rbind(c(1,2,1,0),c(1,1,1,0)),col=c("red","green"))
barplot(rbind(c(0,0,0,2),c(0,0,0,1)),col=c("green","red"), add=T, axes=F)
===
Проблема в том, что barplot() раскрашивает _каждый_ столбик одинаково, так, как определено в векторе col. У Вас просто читались первые два значения (по столбцам), остальные игнорировались.



Pryanik 06.04.2010 09:58
хм... как-то сильно сложно для R.

А как это работает? Почему первый плот вызывает новое графическое окно, а второй нет?



PS2004R 06.04.2010 11:49
(Pryanik @ 06.04.2010 09:58)
Ссылка на исходное сообщение  хм... как-то сильно сложно для R.

А как это работает? Почему первый плот вызывает новое графическое окно, а второй нет?


имеет место быть параметр add = TRUE из ?barplot видим, что он "logical specifying if bars should be added to an already existing plot; defaults to ‘FALSE’."



Pryanik 06.04.2010 19:41
упс... читать разучилась )))



PS2004R 07.04.2010 14:25
Обнаружил на ютубе кучу роликов "про R"

Как мотивация изучения R понравился вот этот, иногда говорю практически также smile.gif

http://www.youtube.com/watch?v=W2GZFeYGU3s&feature=related ( http://www.youtube.com/watch?v=W2GZFeYGU3s&feature=related )



PS2004R 06.05.2010 11:57
Одолел склероз, помогите найти frown.gif

В R несколько систем графики: базовая, latice, grid, iplot ...

Есть еще одна система типа latice.

Предложена относительно недавно, реализована в R в виде языка описания сложных графиков. Язык на довольно солидной теоретической основе. По нему автор реализации R несколько презентаций выложил с состоявшихся воркшопов. Красивые сложные графики. Например карта южной и центральной америк + pcp годовых температур + выделенные на pcp типы кривых температур делят территорию карты на зоны, при этом команда в одну строку помещается. Автор реализации выпустил бумажную книгу.

Я эти презентации увы только просмотрел и не сохранил. Теперь не могу вспомнить ни одного ключевого слова и соответственно найти в поисковиках ничего разумного.

И хотел же сюда запостить, но подумал что не совсем будет по теме frown.gif

Зато теперь по теме получается smile.gif



plantago 06.05.2010 13:32
http://had.co.nz/ggplot2/ ( http://had.co.nz/ggplot2/ ) ?



Guest 08.05.2010 09:30
Скажите, пожалуйста, есть ли версия R для Win98?



PS2004R 08.05.2010 12:56
(Guest @ 08.05.2010 09:30)
Ссылка на исходное сообщение  Скажите, пожалуйста, есть ли версия R для Win98?


в FAQ написано

"R 2.6.2 was the last version for Windows 95. 98, ME and NT4."

ps http://cran.r-project.org/bin/windows/base/old/2.6.2 ( http://cran.r-project.org/bin/windows/base/old/2.6.2 )



dlinnosheee 09.05.2010 20:40
Скажите пожалуйста, как в R рещается такая задача:

построение графика, который меняется с течением времени (есть несколько временнЫх точек, и для каждой точки свой график, хочется аостроить киношку как они друг друга последовательно сменяют.

Желательно самое просто решение. Может быть существует готовый скрипт который это делает?



plantago 09.05.2010 23:22
Постройте несколько графиков, отконвертируйте в gif, сделайте animated gif.



dlinnosheee 10.05.2010 00:40
(plantago @ 10.05.2010 00:22)
Ссылка на исходное сообщение  Постройте несколько графиков, отконвертируйте в gif, сделайте animated gif.

А если у меня сотни графиков?
И потом, гиф не катит, мне надо с ними поиграться, параметры поменять, все такое



plantago 10.05.2010 03:01
Тогда можно http://finzi.psych.upenn.edu/R/library/ani...on-package.html ( http://finzi.psych.upenn.edu/R/library/animation/html/animation-package.html )



PS2004R 10.05.2010 08:28
(dlinnosheee @ 10.05.2010 00:40)
Ссылка на исходное сообщение  А если у меня сотни графиков?
И потом, гиф не катит, мне надо с ними поиграться, параметры поменять, все такое


Напишите функцию выдающую анимированный гиф и вызов просмотрщика для него (или просто последовательность графиков с паузами).

Если хочется графического интерфеса, вызывайте функцию из rpanel, задав параметры в виде всяческих ползунков.



Guest 11.05.2010 14:37
(PS2004R @ 08.05.2010 12:56)
Ссылка на исходное сообщение  в FAQ написано

"R 2.6.2 was the last version for Windows 95. 98, ME and NT4."

ps http://cran.r-project.org/bin/windows/base/old/2.6.2 ( http://cran.r-project.org/bin/windows/base/old/2.6.2 )

Спасибо!



dlinnosheee 12.05.2010 00:44
с последним цветочком промахнулся smile.gif



Игорь Черниенко 02.06.2010 06:52
Рад приветствовать!
А нет ли возможности в CircStats (или чем-то другом) сориентировать круговые и розы диаграммы по азимуту?



plantago 02.06.2010 08:22
В документации нет. В коде тоже. Вообще, rose.diag() функция несложная, можно попробовать модифицировать код. Там только аккуратно надо с lines().
===
Upd. В пакете circular есть своя rose.diag() и там есть параметр rotate wink.gif



Guest 02.06.2010 12:32
В circular все даже проще - при создании вектора соответствующего типа задаешь template="geographics" и будет счастье. Спасибо за наводку!



Lea1 17.06.2010 17:29
Добрый день.

Мне нужно генерировать 2895 initial values for Winbug, который дожны выглядеть как 0,0,0,0,0... решила воспользоваться R. Мои действия:

a<-matrix(0,nrow=2895,ncol=1)
write.table(a, file="d:/.../test.txt", sep=",", row.names = FALSE, col.names = TRUE, quote=FALSE)


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



PS2004R 17.06.2010 17:37
(Lea1 @ 17.06.2010 17:29)
Ссылка на исходное сообщение  Добрый день.

Мне нужно генерировать 2895 initial values for Winbug, который дожны выглядеть как 0,0,0,0,0... решила воспользоваться R. Мои действия:

a<-matrix(0,nrow=2895,ncol=1)
write.table(a, file="d:/.../test.txt", sep=",", row.names = FALSE, col.names = TRUE, quote=FALSE)
Но получила
0
0
0...
Очень нужны запятые, вручную простовлять утомительно. Не понимаю, почему не получились. Помогите. плжалуйста,  добавить запятые в текстовый файл


write.table(t(rep(0,2000)), file="/tmp/test.txt", sep=",", row.names = FALSE, col.names = FALSE, quote=FALSE)

если идти Вашим путем smile.gif

upd или если я вдруг не так понял

write.table(cbind(rep(0,2000),rep('',2000)), file="/tmp/test.txt", sep=",", row.names = FALSE, col.names = FALSE, quote=FALSE)



Игорь Черниенко 18.06.2010 09:31
Здравствуйте.
Можно ли в R сравнивать коэффициенты линейной регрессии? Для 2-х уравнений, а лучше - более чем 2-х.



plantago 18.06.2010 10:23
Не вполне понял Вас, но если Вы возьмете lm$coefficients, то увидите все коэффициенты данной линейной модели, а дальше можно запомнить их в переменную и сравнивать.



PS2004R 18.06.2010 14:31
(Игорь Черниенко @ 18.06.2010 09:31)
Ссылка на исходное сообщение  Здравствуйте.
Можно ли в R сравнивать коэффициенты линейной регрессии? Для 2-х уравнений, а лучше - более чем 2-х.


Две модели можно сравнивать между собой.

m1 <- lm( .... )
m2 <- lm( .... )

anova(m1,m2)



Cucumaria 18.06.2010 22:14
А еще есть специальный тест, чтобы сравнить углы наклона прямых в регрессиях, только я не помню, как он называется. И в Р он есть, надо вспомнить. Для сравнения двух коэффициентов мне попадался тест, а вот чтобы больше...



Cucumaria 18.06.2010 22:18
Эта функция rdif.nul в пакете psychometric.



PS2004R 19.06.2010 15:39
(Cucumaria @ 18.06.2010 22:18)
Ссылка на исходное сообщение  Эта функция rdif.nul в пакете psychometric.


всё таки он сравнивает два коэффициента корреляции

rdif.nul(r1, r2, n1, n2)

Arguments
r1 Correlation 1
r2 Correlation 2
n1 Sample size for r1
n2 Sample size for r2
Details

First converts r to z' for each correlation. Then constructs a z test for the difference z <- (z1 - z2)/sqrt(1/(n1-3)+1/(n2-3))



dlinnosheee 18.07.2010 19:36
а не пробегала ли ту случайно книга в pdf по R и биокондуктор?
желательно что-нибудь из этого:
http://www.bioconductor.org/pub ( http://www.bioconductor.org/pub )



voliadis 20.07.2010 09:29
Собрал описания 100 книг по использованию R. Книги на английском языке. Полных текстов нет. Есть только ссылки на каждую книгу в Google Books и amazon

http://voliadis.ru/bibliography/r/list ( http://voliadis.ru/bibliography/r/list )



plantago 20.07.2010 10:07
http://www.r-project.org/doc/bib/R-books.html ( http://www.r-project.org/doc/bib/R-books.html )



bubnilkin 23.07.2010 15:44
Добрый день!

У меня есть значения SE, CI и M можно ли построить ящик с усами в R только на этих данных (т.е. БЕЗ вариационных рядов)?



PS2004R 23.07.2010 17:46
(bubnilkin @ 23.07.2010 15:44)
Ссылка на исходное сообщение  Добрый день!

У меня есть значения SE, CI и M можно ли построить ящик с усами в R только на этих данных (т.е. БЕЗ вариационных рядов)?


?bxp



bubnilkin 02.08.2010 12:11
а как отобразить список файлов с определённым разрешением (.fcs)?

функция

list.files(path = "...", pattern = "[*.fcs]")

выдаёт все файлы в директории path почему-то...



plantago 02.08.2010 12:39
list.files(path = "...", pattern = "*\\.fcs")
или
list.files(path = "...", pattern = glob2rx("*.fcs"))



bubnilkin 02.08.2010 13:32
благодарю smile.gif с первой строкой с 4-раза получилось smile.gif



voliadis 06.08.2010 20:15
Здравствуйте,

В описании функции aov() сказано следующее:

Fit an analysis of variance model by a call to ‘lm’ for each stratum.

О чем идет речь я понимаю, но выразить по-русски не могу. Подскажите, пожалуйста, правильный перевод этой фразы. teapot.gif



tensor 06.08.2010 21:11
(voliadis @ 06.08.2010 20:15)
Ссылка на исходное сообщение  Здравствуйте,

В описании функции aov() сказано следующее:

Fit an analysis of variance model by a call to ‘lm’ for each stratum.

О чем идет речь я понимаю, но выразить по-русски не могу. Подскажите, пожалуйста, правильный перевод этой фразы. teapot.gif

Оценка {сходства|соответствия} [распределений] с использованием {анализа дисперсии|дисперсионного анализа} {вызовом|обращением к} [функции] 'lm' для каждой {из} групп{ы}.
В квадратных скобках то, что подразумевается. В фигурных - варианты использования.
PS Спасибо Вам за те книжки, что вы мне когда-то переслали! До сих пор успешно использую.



dlinnosheee 07.08.2010 13:42
Video tutorial. Statistics with R

part 1: vector arithmetics tutorial: http://www.youtube.com/watch?v=mL27TAJGlW ( http://www.youtube.com/watch?v=mL27TAJGlW )
part 2: help search commands tutorial: http://www.youtube.com/watch?v=d6nSr9L9bVU ( http://www.youtube.com/watch?v=d6nSr9L9bVU )
part 3: plot and history tutorial: http://www.youtube.com/watch?v=NfH5peM1RtI ( http://www.youtube.com/watch?v=NfH5peM1RtI )
part 4: R CRAN web tutorial: http://www.youtube.com/watch?v=AipnE4s8sKk ( http://www.youtube.com/watch?v=AipnE4s8sKk )
part 5: external files tutorial: http://www.youtube.com/watch?v=qb3LpWxvubQ ( http://www.youtube.com/watch?v=qb3LpWxvubQ )
part 6: matrix operations tutorial: http://www.youtube.com/watch?v=GFsgBqjAMvQ ( http://www.youtube.com/watch?v=GFsgBqjAMvQ )
part 7: large dataset tutorial: http://www.youtube.com/watch?v=B5h9MHqpGOQ ( http://www.youtube.com/watch?v=B5h9MHqpGOQ )



dlinnosheee 07.08.2010 14:44
У меня вопрос: использует ли кто-то тут R commander?



PS2004R 07.08.2010 21:10
(dlinnosheee @ 07.08.2010 14:44)
Ссылка на исходное сообщение  У меня вопрос: использует ли кто-то тут R commander?


Я понимаю что отвечаю не на вопрос (:

Наверное всё таки Sweave() более эффективная техника. Если использовать связку emacs+ess+austex+ reftex+noweb и какую нибудь из систем контроля версий, то работа вообще становится праздником (:

Да и сам подход, когда отбрасывается первое достоинство R --- язык, в котором любая манипуляция с данными занимает как правило строчку кода максимум в 100 символов (если без форматирования <:), и имеется встроенная автокомплектация всего что только можно...

Уж лучше иметь шаблоны типового анализа для Sweave().



plantago 08.08.2010 00:05
(dlinnosheee @ 07.08.2010 07:44)
Ссылка на исходное сообщение  У меня вопрос: использует ли кто-то тут R commander?

Я его на русский перевожу smile.gif



dlinnosheee 08.08.2010 11:50
(plantago @ 08.08.2010 01:05)
Ссылка на исходное сообщение  Я его на русский перевожу smile.gif

У меня такой глюк с ним: пытаюсь сделать Load File, нацеливаю его на текстовый файл (большой текстовый файл который распакован из зип файла, и другими вьюверами этот текстовый файл спокойно читается), а в ответ R commander пишет "bad restore file magic number (file may be corrupted) -- no data loaded". И все на этом, загрузить этот файл я не могу. Не встречались с таким?



PS2004R 10.08.2010 11:34
(dlinnosheee @ 08.08.2010 11:50)
Ссылка на исходное сообщение  У меня такой глюк с ним: пытаюсь сделать Load File, нацеливаю его на текстовый файл (большой текстовый файл который распакован из зип файла, и другими вьюверами этот текстовый файл спокойно читается), а в ответ R commander пишет "bad restore file magic number (file may be corrupted) -- no data loaded". И все на этом, загрузить этот файл я не могу. Не встречались с таким?


судя по всему оно использует load() , а не что то типа read.table()...


?load
Reload datasets written with the function ‘save’



plantago 10.08.2010 19:56
Надо использовать не load, а import.



flegmatic 17.08.2010 12:03
Как в R можно решить следующую задачу?

Имеются две переменные: $esan (групирующий фактор с 8-ю уровнями) и $reus.

Пример:
"esan" "reus"
a 9.3
a 9.05
a 7.78
b 7.11
b 7.14
c 8.12
c 7.5
d 7.84
e 7.8
f 7.52
g 8.84
g 6.98
h 6.1
h 6.89

Нужно сравнить попарно все 8 уровней/подгруп с помощью wilcox.test по переменной "reus" и получить результат в виде таблицы с p-значениями примерно такого вида:

a b c
b 0.773 - -
c 0.112 0.084 -
d 0.908 0.848 0.095

Что-то похожее нашёл здесь: https://stat.ethz.ch/pipermail/r-help/2007-July/136627.html ( https://stat.ethz.ch/pipermail/r-help/2007-July/136627.html ) (если убрать переменную "Grp"), но адаптировать не сумел.



PS2004R 17.08.2010 15:24
(flegmatic @ 17.08.2010 12:03)
Ссылка на исходное сообщение  Как в R можно решить следующую задачу?

Имеются две переменные: $esan (групирующий фактор с 8-ю уровнями) и $reus.

Пример:
"esan" "reus"
a 9.3
a 9.05
a 7.78
b 7.11
b 7.14
c 8.12
c 7.5
d 7.84
e 7.8
f 7.52
g 8.84
g 6.98
h 6.1
h 6.89

Нужно сравнить попарно все 8 уровней/подгруп с помощью wilcox.test по переменной "reus" и получить результат в виде таблицы с p-значениями примерно такого вида:

    a    b    c   
b 0.773 -    -   
c 0.112 0.084 -   
d 0.908 0.848 0.095

Что-то похожее нашёл здесь: https://stat.ethz.ch/pipermail/r-help/2007-July/136627.html ( https://stat.ethz.ch/pipermail/r-help/2007-July/136627.html ) (если убрать переменную "Grp"), но адаптировать не сумел.



> myDat <- read.table("data.txt", head= TRUE)
> myDat
esan reus
1 a 9.30
2 a 9.05
3 a 7.78
4 b 7.11
5 b 7.14
6 c 8.12
7 c 7.50
8 d 7.84
9 e 7.80
10 f 7.52
11 g 8.84
12 g 6.98
13 h 6.10
14 h 6.89

> pairwise.wilcox.test(as.numeric(as.character(myDat$reus)),
+ myDat$esan,p.adjust.method = "none")$p.value

a b c d e f g
b 0.2 NA NA NA NA NA NA
c 0.4 0.3333333 NA NA NA NA NA
d 1.0 0.6666667 1.0000000 NA NA NA NA
e 1.0 0.6666667 1.0000000 1.0000000 NA NA NA
f 0.5 0.6666667 1.0000000 1.0000000 1.0000000 NA NA
g 0.4 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 NA
h 0.2 0.3333333 0.3333333 0.6666667 0.6666667 0.6666667 0.3333333
>



PS2004R 17.08.2010 18:17
Но методологически так делать не правильно. Тест должен только подтверждать, или опровергать, уже выдвинутую гипотезу, а не служить причиной её выдвижения.



flegmatic 17.08.2010 19:13
(PS2004R @ 17.08.2010 15:24)
Ссылка на исходное сообщение
> pairwise.wilcox.test(as.numeric(as.character(myDat$reus)),
+ myDat$esan,p.adjust.method = "none")$p.value


Спасибо! Но, если воспроизводить точь-в-точь Ваше решение (PS2004R), то всё получается, а когда пытаюсь применить данное решение на своих конкретных данных, то получаю "ERROR: not enough (finite) 'x' observations" (это происходит в R Commander). Мой рабочий файл (datа.csv) отличается от data.txt только количеством переменных и обьёмом выборки, поэтому не вижу в чём проблема.

myDat <- read.table("datа.csv", head=TRUE, sep=",", na.strings="NA", dec=".",
+ strip.white=TRUE)
myDat
pairwise.wilcox.test(as.numeric(as.character(myDat$reus)),
+ myDat$esan,p.adjust.method = "none")$p.value



PS2004R 17.08.2010 19:29
(flegmatic @ 17.08.2010 19:13)
Мой рабочий файл (datа.csv) отличается от data.txt только количеством переменных и обьёмом выборки, поэтому не вижу в чём проблема.



as.numeric() должна решать эту проблему...

Если данные не секретные выложите их, проблема похоже в них и их преобразовании.



flegmatic 17.08.2010 20:27
Данные не секретные.

Не могу сообразить как as.numeric() может решить проблему frown.gif.



Файл/ы:

Скачать файл data.txt
Размер:: 15.05к
кол-во скачиваний: 233




PS2004R 17.08.2010 21:07
[quote=flegmatic,17.08.2010 20:27]Ссылка на исходное сообщение  Данные не секретные.

Не могу сообразить как as.numeric() может решить проблему frown.gif.
[/]


Прилось почистить от NA

> data <- read.csv("/home/petrov/uni/proba.r/data1.txt")

> data.na.omit <- na.omit(cbind(data$reus,data$esan))

> pairwise.wilcox.test(data.na.omit[,1],
+ data.na.omit[,2],
+ p.adjust.method = "none")$p.value

1 2 3 4 7
2 3.551600e-04 NA NA NA NA
3 5.368394e-08 0.004353211 NA NA NA
4 1.467191e-10 0.004062074 0.56463777 NA NA
7 5.171187e-03 0.361756594 0.34506796 0.53028646 NA
8 4.695635e-02 0.938276125 0.04076322 0.06934225 0.4114888

Предупреждения
1: In wilcox.test.default(xi, xj, paired = paired, ...) :
не могу подсчитать точное p-значение при наличии повторяющихся наблюдений
2: In wilcox.test.default(xi, xj, paired = paired, ...) :
не могу подсчитать точное p-значение при наличии повторяющихся наблюдений
3: In wilcox.test.default(xi, xj, paired = paired, ...) :
не могу подсчитать точное p-значение при наличии повторяющихся наблюдений



flegmatic 17.08.2010 22:17
Спасибо большое, PS2004R!

Всё получается!

И последнее. Не хочется злоупотреблять Вашим вниманием, но хотелось бы решить и другую проблему. Как применить Ваш код (вместе с чисткой от NA) ко всем переменным, т.е. обработать одним скриптом всю таблицу?



flegmatic 18.08.2010 10:02
Для PS2004R и др.

Вот что получилось на данный момент (замечание: код перенят и адаптирован, но СПАМ-фильтр почему-то не разрешает вставку ссылки. Страницу можно найти при помощи Google по названию: "Beautiful Correlation Tables in R"):
myDat <- read.csv("/home/xxxx/Documents/data.txt")
myDat.na.omit <- na.omit(cbind(myDat$reus,myDat$esan))
P <- pairwise.wilcox.test(as.numeric(as.character(myDat.na.omit[,1])),
+ myDat.na.omit[,2],p.adjust.method = "none")$p.value
## truncate the matrix that holds the p-values to three decimal
P <- format(round(cbind(rep(1.111, ncol(P)), P), 3))[,-1]
## define notions for significance levels; spacing is important.
mystars <- ifelse(P < .001, "***", ifelse(P < .01, "** ", ifelse(P < .05, "*  ", "   ")))
## build a new matrix that includes the p-values with their apropriate stars
Pnew <- matrix(paste(P, mystars, sep=""), ncol=ncol(P))
Pnew <- as.data.frame(Pnew)
rownames(Pnew) <- rownames(P)
colnames(Pnew) <- paste(colnames(P), "", sep="")
Pnew

Это выдаваемый результат:
          1        2        3        4        7
2 0.000***    NA       NA       NA       NA  
3 0.000*** 0.004**     NA       NA       NA  
4 0.000*** 0.004**  0.565       NA       NA  
7 0.005**  0.362    0.345    0.530       NA  
8 0.047*   0.938    0.041*   0.069    0.411

Хотелось бы решить следующее (умения пока не хватает):
1. Почистить таблицу от NA.
2. Превратить данный код в функцию.
3. И, как писал выше, сделать возможным одновременное применение данной функции к n-количеству переменных.



Файл/ы:

Скачать файл data.txt
Размер:: 15.05к
кол-во скачиваний: 162




PS2004R 18.08.2010 10:34
(flegmatic @ 18.08.2010 10:02)
Хотелось бы решить следующее (умения пока не хватает):
1. Почистить таблицу от NA.
2. Превратить данный код в функцию.
3. И, как писал выше, сделать возможным одновременное применение данной функции к n-количеству переменных.


Чистить придётся попарно, иначе при удачном стечении обстоятельсв от таблицы просто ничего не останется.

функция:

CODE

pw <- function (f,x) {data.na.omit <- na.omit(cbind(f,x))
                              pairwise.wilcox.test(data.na.omit[,1],
                                                          data.na.omit[,2],
                                                          p.adjust.method = "none")$p.value}
pw(data$reus,data$esan)


просто в цикл её поместить, поскольку Вам не надо объединять результат в одну структуру.

for (i in c(...например номера колонок переменных....)) pw(группирующая, data(,i))



flegmatic 18.08.2010 12:02
(PS2004R @ 18.08.2010 10:34)
Ссылка на исходное сообщение
CODE

pw <- function (f,x) {data.na.omit <- na.omit(cbind(f,x))
                              pairwise.wilcox.test(data.na.omit[,1],
                                                          data.na.omit[,2],
                                                          p.adjust.method = "none")$p.value}
pw(data$reus,data$esan)


Функция работает, спасибо, но цикл не получается. Делал три варианта, все нерабочие:
CODE

for(i in c(length(names(myDat)))) pw(myDat$esan, myDat[,i])
for(i in seq(3, length(names(myDat)), 1)) pw(myDat$esan, myDat[,i])
for(i in 3:length(names(myDat))) pw(myDat$esan, myDat[,i])




PS2004R 18.08.2010 12:55
CODE


> names(myDat)
[1] "esan" "reus"
> length(names(myDat))
[1] 2
> 1:length(names(myDat))
[1] 1 2

> pw <- function (f,x) {data.na.omit <- na.omit(cbind(f,x))
+                       print(pairwise.wilcox.test(data.na.omit[,1],
+                                                  data.na.omit[,2],
+                                                  p.adjust.method = "none")$p.value)}
> for (i in 2:length(names(myDat))) pw(data$reus, data[,i])
            1           2          3          4         7
2 3.551600e-04          NA         NA         NA        NA
3 5.368394e-08 0.004353211         NA         NA        NA
4 1.467191e-10 0.004062074 0.56463777         NA        NA
7 5.171187e-03 0.361756594 0.34506796 0.53028646        NA
8 4.695635e-02 0.938276125 0.04076322 0.06934225 0.4114888
Предупреждения
1: In wilcox.test.default(xi, xj, paired = paired, ...) :
 cannot compute exact p-value with ties
2: In wilcox.test.default(xi, xj, paired = paired, ...) :
 cannot compute exact p-value with ties
3: In wilcox.test.default(xi, xj, paired = paired, ...) :
 cannot compute exact p-value with ties
>



Крайне рекомендую "An Introduction to R". Я его просто распечатал в своё время, переплел и последовательно прочитал (оно до сих пор у меня есть smile.gif. Там очень хорошая база для последующего использования возможностей R.

Время на изучение этого руководства никогда в последующем не будет Вами оценено как потраченное зря.



flegmatic 18.08.2010 13:50
(PS2004R @ 18.08.2010 12:55)
Ссылка на исходное сообщениеКрайне рекомендую "An Introduction to R". Я его просто распечатал в своё время, переплел и последовательно прочитал (оно до сих пор у меня есть smile.gif. Там очень хорошая база  для последующего использования возможностей R.

Время на изучение этого руководства никогда в последующем не будет Вами оценено как потраченное зря.

Спасибо! Я с Вами полностью согласен! Некоторое время назад я начинал его изучать, но так и не довёл это дело до конца. Обещаю (сам себе) это обязательно сделать!



PS2004R 19.08.2010 00:25
http://rnd.cnews.ru/army/news/top/index_sc...10/08/16/405340 ( http://rnd.cnews.ru/army/news/top/index_science.shtml?2010/08/16/405340 )

R в новостях smile.gif



tensor 19.08.2010 10:12
Простите, что лезу с небольшой критикой, но тут надо будет вносить поправку Бонферрони, потому что возникнет эффект множественных сравнений. И после ее внесения уровни статистической значимости круто изменятся на гораздо менее оптимистичные. Так что осторожнее, а то можно ненароком себя и других обмануть.



flegmatic 19.08.2010 10:16
1) tensor, спасибо за помощь!

2) Следующий код (спасибо PS2004R!):
CODE
pw <- function(f,x) {data.na.omit <- na.omit(cbind(f,x))
P <- pairwise.wilcox.test(data.na.omit[,1], data.na.omit[,2],
p.adjust.method = "none", paired=FALSE)$p.value
## truncate the matrix that holds the p-values to three decimal
P <- format(round(cbind(rep(1.111, ncol(P)), P), 3))[,-1]
## define notions for significance levels; spacing is important.
stars <- ifelse(P < .001, "***", ifelse(P < .01, "** ", ifelse(P < .05, "*  ", "   ")))
## build a new matrix that includes the p-values with their apropriate stars
Pnew <- matrix(paste(P, stars, sep=""), ncol=ncol(P))
Pnew <- as.data.frame(Pnew)
rownames(Pnew) <- rownames(P)
colnames(Pnew) <- paste(colnames(P), "", sep="")
print(Pnew)}

for(i in 3:length(names(myDat))) {
cat("\n")
print(names(myDat[i]))
pw(myDat[,i], myDat$esantionul)}

выдаёт следующий результат (кол-во знаков после запятой сокращено до 3-х, звёздочками обозначены три уровня статистической значимости):
CODE

[1] "reus"
        1        2        3        4        7
2 0.000***    NA       NA       NA       NA  
3 0.000*** 0.004**     NA       NA       NA  
4 0.000*** 0.004**  0.565       NA       NA  
7 0.005**  0.362    0.345    0.530       NA  
8 0.047*   0.938    0.041*   0.069    0.411  

[1] "larg"
        1        2        3        4        5        6        7
2 0.542       NA       NA       NA       NA       NA       NA  
3 0.222    0.100       NA       NA       NA       NA       NA  
4 0.511    0.220    0.460       NA       NA       NA       NA  
5 0.109    0.344    0.012*   0.030*      NA       NA       NA  
6 0.027*   0.133    0.003**  0.006**  0.508       NA       NA  
7 0.713    0.469    0.687    1.000    0.123    0.033*      NA  
8 0.931    0.869    0.408    0.629    0.259    0.081    0.759

Видно, что строки и колонки (групирующая переменная) обозначены цифрами, причём в первой таблице отсутствуют строки/колонки номер 5 и 6. Каждое цифровое обозначение указывает на конкретную подгрупу наблюдений. Так как цифры слишком уж безликие имена, хотелось бы заменить их на более легко воспринимаемые/понимаемые (в первую очередь самим исследователем smile.gif). Например: 1 заменить на I, 2 на II, 3 на III, 4 на IV, 5 на 1av, 6 на 2av, 7 на gr1, 8 на gr2. Выполнить замену в самой групирующей переменной не получилось, R выдаёт что она должна быть числовой. Значит, замену нужно выполнить при помощи кода? Но как это сделать?



PS2004R 19.08.2010 10:58
CODE

p.adjust(p, method = p.adjust.methods, n = length(p))
   
    p.adjust.methods
    # c("holm", "hochberg", "hommel", "bonferroni", "BH", "BY",
    #   "fdr", "none")


CODE


data$esan[data$esan==8] <- "acht"
> data$esan
 [1] "1"    "1"    "1"    "1"    "1"    "1"    "1"    "1"    "1"    "1"  
[11] "1"    "1"    "1"    "1"    "1"    "1"    "1"    "1"    "1"    "1"  
[21] "1"    "1"    "1"    "1"    "1"    "1"    "1"    "1"    "1"    "1"  
[31] "1"    "1"    "1"    "1"    "1"    "1"    "1"    "1"    "1"    "1"  
[41] "1"    "1"    "1"    "1"    "1"    "1"    "1"    "1"    "1"    "1"  
[51] "1"    "1"    "1"    "1"    "1"    "1"    "1"    "1"    "1"    "1"  
[61] "1"    "1"    "1"    "1"    "1"    "1"    "1"    "1"    "1"    "1"  
[71] "1"    "5"    "5"    "5"    "5"    "5"    "5"    "5"    "5"    "5"  
[81] "5"    "5"    "5"    "5"    "5"    "5"    "5"    "5"    "2"    "2"  
[91] "2"    "2"    "2"    "2"    "2"    "2"    "2"    "2"    "2"    "2"  
[101] "2"    "2"    "2"    "2"    "2"    "2"    "2"    "2"    "2"    "2"  
[111] "2"    "2"    "2"    "2"    "2"    "2"    "2"    "2"    "2"    "2"  
[121] "2"    "2"    "2"    "2"    "2"    "2"    "2"    "2"    "2"    "2"  
[131] "2"    "2"    "2"    "2"    "2"    "2"    "2"    "2"    "2"    "2"  
[141] "2"    "2"    "2"    "6"    "6"    "6"    "6"    "6"    "6"    "6"  
[151] "6"    "6"    "6"    "6"    "6"    "6"    "3"    "3"    "3"    "3"  
[161] "3"    "3"    "3"    "3"    "3"    "3"    "3"    "3"    "3"    "3"  
[171] "3"    "3"    "3"    "3"    "3"    "3"    "3"    "3"    "3"    "3"  
[181] "3"    "3"    "3"    "3"    "3"    "3"    "3"    "4"    "4"    "4"  
[191] "4"    "4"    "4"    "4"    "4"    "4"    "4"    "4"    "4"    "4"  
[201] "4"    "4"    "4"    "4"    "4"    "4"    "4"    "4"    "4"    "4"  
[211] "4"    "4"    "4"    "4"    "4"    "4"    "4"    "4"    "4"    "4"  
[221] "4"    "4"    "4"    "4"    "4"    "4"    "4"    "4"    "4"    "4"  
[231] "4"    "4"    "4"    "4"    "4"    "4"    "4"    "4"    "4"    "4"  
[241] "4"    "4"    "4"    "4"    "4"    "4"    "4"    "4"    "4"    "4"  
[251] "4"    "4"    "7"    "7"    "7"    "7"    "7"    "7"    "7"    "7"  
[261] "7"    "7"    "7"    "acht" "acht" "acht" "acht" "acht" "acht" "acht"
[271] "acht" "acht" "acht" "acht"



ну и преобразовать в фактор потом

CODE


attach(airquality)
    Month <- factor(Month, labels = month.abb[5:9])
    ## These give warnings because of ties :
    pairwise.wilcox.test(Ozone, Month)
    pairwise.wilcox.test(Ozone, Month, p.adj = "bonf")
    detach()


CODE

labels: _either_ an optional vector of labels for the levels (in the
         same order as ‘levels’ after removing those in ‘exclude’),
         _or_ a character string of length 1.




flegmatic 19.08.2010 11:27
(tensor @ 19.08.2010 10:12)
Ссылка на исходное сообщение  Простите, что лезу с небольшой критикой, но тут надо будет вносить поправку Бонферрони, потому что возникнет эффект множественных сравнений. И после ее внесения уровни статистической значимости круто изменятся на гораздо менее оптимистичные. Так что осторожнее, а то можно ненароком себя и других обмануть.

Кошмар Бонферрони.ppt (4Mb) ( http://www.vigg.ru/get_file.php?id=196 )

P.S. Извините за отклонение от главной темы, но материал показался интересным!
Для просмотревших: как применяются альтернативы для поправки Бонферрони (и для других поправок, наверное), указанные в данной презентации, средствами R?

PS2004R,
CODE

p.adjust(p, method = p.adjust.methods, n = length(p))
 
   p.adjust.methods
   # c("holm", "hochberg", "hommel", "bonferroni", "BH", "BY",
   #   "fdr", "none")

"fdr" в Вашем сообщении и "FDR-контроль" в презентации это одно и то же?



flegmatic 19.08.2010 12:04
PS2004R, делаю таким образом:
CODE
myDat$esan[myDat$esan==8] <- "acht"
myDat$esan <- factor(myDat$esan)
>myDat$esan
[1] 1    1    1    1    1    1    1    1    1    1    1    1    1    1    1    1    1    1    1    1    1    1    1    1    1    1    1    1    1    1    1  
[32] 1    1    1    1    1    1    1    1    1    1    1    1    1    1    1    1    1    1    1    1    1    1    1    1    1    1    1    1    1    1    1  
[63] 1    1    1    1    1    1    1    1    1    5    5    5    5    5    5    5    5    5    5    5    5    5    5    5    5    5    2    2    2    2    2  
[94] 2    2    2    2    2    2    2    2    2    2    2    2    2    2    2    2    2    2    2    2    2    2    2    2    2    2    2    2    2    2    2  
[125] 2    2    2    2    2    2    2    2    2    2    2    2    2    2    2    2    2    2    2    6    6    6    6    6    6    6    6    6    6    6    6  
[156] 6    3    3    3    3    3    3    3    3    3    3    3    3    3    3    3    3    3    3    3    3    3    3    3    3    3    3    3    3    3    3  
[187] 3    4    4    4    4    4    4    4    4    4    4    4    4    4    4    4    4    4    4    4    4    4    4    4    4    4    4    4    4    4    4  
[218] 4    4    4    4    4    4    4    4    4    4    4    4    4    4    4    4    4    4    4    4    4    4    4    4    4    4    4    4    4    4    4  
[249] 4    4    4    4    7    7    7    7    7    7    7    7    7    7    7    acht acht acht acht acht acht acht acht acht acht acht
Levels: 1 2 3 4 5 6 7 acht
Похоже что всё в порядке, но после обработки функцией pw получаем ту же картину, т.е. вместо "acht" имеем цифру "8":
CODE
[1] "reus"
        1        2        3        4        7
2 0.002**     NA       NA       NA       NA  
3 0.000*** 0.013*      NA       NA       NA  
4 0.000*** 0.013*   0.605       NA       NA  
7 0.013*   0.493    0.493    0.605       NA  
8 0.088    0.938    0.087    0.116    0.514

Я что-то делаю не так?



PS2004R 19.08.2010 12:52
Вот видите какие "трудности" порождает методологически неправильный путь? smile.gif

Придется преобразовать матрицу в датафрейм (он разрешает разнотипные столбцы-переменные), и прибить в функции всё гвоздями.

CODE

pw <- function (f,x) {data.na.omit <- as.data.frame(na.omit(cbind(f,x)))
                                 print(pairwise.wilcox.test(as.numeric(data.na.omit[,1]),
                                                                       as.factor(data.na.omit[,2]),
                                                                       p.adjust.method = "none")$p.value)}



наверное так лучше

CODE

pw <- function (f,x) {data.na.omit <- na.omit(data.frame(f,x))
                                 print(pairwise.wilcox.test(as.numeric(data.na.omit[,1]),
                                                                        as.factor(data.na.omit[,2]),
                                                                        p.adjust.method = "none")$p.value)}





flegmatic 19.08.2010 13:08
Решил проблему так:
CODE

pw <- function(f,x) {data.na.omit <- na.omit(cbind(f,x))
P <- pairwise.wilcox.test(data.na.omit[,1], data.na.omit[,2],
p.adjust.method = "none", paired=FALSE)$p.value
## truncate the matrix that holds the p-values to three decimal
P <- format(round(cbind(rep(1.111, ncol(P)), P), 3))[,-1]
## define notions for significance levels; spacing is important.
stars <- ifelse(P < .001, "***", ifelse(P < .01, "** ", ifelse(P < .05, "*  ", "   ")))
## build a new matrix that includes the p-values with their apropriate stars
Pnew <- matrix(paste(P, stars, sep=""), ncol=ncol(P))
Pnew <- as.data.frame(Pnew)
rownames(Pnew) <- rownames(P)
colnames(Pnew) <- paste(colnames(P), "", sep="")
[color=green]colnames(Pnew)[colnames(Pnew)==1] <- "I"
colnames(Pnew)[colnames(Pnew)==2] <- "II"
colnames(Pnew)[colnames(Pnew)==3] <- "III"
colnames(Pnew)[colnames(Pnew)==4] <- "IV"
colnames(Pnew)[colnames(Pnew)==5] <- "1av"
colnames(Pnew)[colnames(Pnew)==6] <- "2av"
colnames(Pnew)[colnames(Pnew)==7] <- "gr1"
colnames(Pnew)[colnames(Pnew)==8] <- "gr2"
rownames(Pnew)[rownames(Pnew)==1] <- "I"
rownames(Pnew)[rownames(Pnew)==2] <- "II"
rownames(Pnew)[rownames(Pnew)==3] <- "III"
rownames(Pnew)[rownames(Pnew)==4] <- "IV"
rownames(Pnew)[rownames(Pnew)==5] <- "1av"
rownames(Pnew)[rownames(Pnew)==6] <- "2av"
rownames(Pnew)[rownames(Pnew)==7] <- "gr1"
rownames(Pnew)[rownames(Pnew)==8] <- "gr2"[/color]
print(Pnew)}
Решение не совсем нравится, т.к. теряется универсальность функции, но на данный момент оно меня устраивает.



flegmatic 19.08.2010 13:12
(PS2004R @ 19.08.2010 12:52)
Ссылка на исходное сообщение  Вот видите какие "трудности" порождает методологически неправильный путь? smile.gif

Придется преобразовать матрицу в датафрейм (он разрешает разнотипные столбцы-переменные), и прибить в функции всё гвоздями.


PS2004R, Ваше последнее сообщение ещё изучаю, но спасибо.



PS2004R 19.08.2010 13:39
Огромное спасибо за красивую ссылку. Метод "коврового бомбометания" тестами ужасен по своей сути frown.gif.

Фактически ситуация когда способ проверки результата применяется как метод выявления (разведки).

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



Я предпочитаю применять бутсреп (в конце Автор презентации об этом и говорит smile.gif. Даже простейший метод складного ножа с исключением одной, двух вариант в малых выборках способен выявить сомнительные случаи.


результат ?p.adjust

CODE


The adjustment methods include the Bonferroni correction
    (‘"bonferroni"’) in which the p-values are multiplied by the
    number of comparisons.  Less conservative corrections are also
    included by Holm (1979) (‘"holm"’), Hochberg (1988)
    (‘"hochberg"’), Hommel (1988) (‘"hommel"’), Benjamini & Hochberg
    (1995) (‘"BH"’ or its alias ‘"fdr"’), and Benjamini & Yekutieli
    (2001) (‘"BY"’), respectively.  A pass-through option (‘"none"’)
    is also included.  The set of methods are contained in the
    ‘p.adjust.methods’ vector for the benefit of methods that need to
    have the method as an option and pass it on to ‘p.adjust’.

    The first four methods are designed to give strong control of the
    family wise error rate.  There seems no reason to use the
    unmodified Bonferroni correction because it is dominated by Holm's
    method, which is also valid under arbitrary assumptions.

    Hochberg's and Hommel's methods are valid when the hypothesis
    tests are independent or when they are non-negatively associated
    (Sarkar, 1998; Sarkar and Chang, 1997).  Hommel's method is more
    powerful than Hochberg's, but the difference is usually small and
    the Hochberg p-values are faster to compute.

    The ‘"BH"’ (aka ‘"fdr"’) and ‘"BY"’ method of Benjamini, Hochberg,
    and Yekutieli control the false discovery rate, the expected
    proportion of false discoveries amongst the rejected hypotheses.
    The false discovery rate is a less stringent condition than the
    family wise error rate, so these methods are more powerful than
    the others.

    Note that you can set ‘n’ larger than ‘length(p)’ which means the
    unobserved p-values are assumed to be greater than all the
    observed p for ‘"bonferroni"’ and ‘"holm"’ methods and equal to 1
    for the other methods.





PS2004R 19.08.2010 13:43
(flegmatic @ 19.08.2010 13:12)
Ссылка на исходное сообщение  PS2004R, Ваше последнее сообщение ещё изучаю, но спасибо.


У меня оно выводит вот так, как и в самом первом примере.

CODE


> pw <- function (f,x) {data.na.omit <- na.omit(data.frame(f,x))
+                                  print(pairwise.wilcox.test(as.numeric(data.na.omit[,1]),
+                                                                        as.factor(data.na.omit[,2]),
+                                                                        p.adjust.method = "none")$p.value)}
> for (i in 2:length(names(myDat))) pw(data$reus, data[,i])
               1           2          3          4         7
2    3.551600e-04          NA         NA         NA        NA
3    5.368394e-08 0.004353211         NA         NA        NA
4    1.467191e-10 0.004062074 0.56463777         NA        NA
7    5.171187e-03 0.361756594 0.34506796 0.53028646        NA
acht 4.695635e-02 0.938276125 0.04076322 0.06934225 0.4114888





flegmatic 19.08.2010 13:52
PS2004R, большое спасибо, Вы мне очень помогли!!!
Вот что получилось:
CODE
> myDat <- read.csv("data.csv")
> myDat$esan <- factor(myDat$esan, labels=c('I','II','III','IV','1av','2av','gr1','gr2'))
> pw <- function (f,x) {data.na.omit <- na.omit(data.frame(f,x))
+ print(pairwise.wilcox.test(as.numeric(data.na.omit[,1]), data.na.omit[,2],
+ p.adjust.method = "none", paired=FALSE)$p.value)}
> for(i in 4:length(names(myDat))-1) {
+ cat("\n")
+ print(names(myDat[i]))
+ pw(myDat[,i], myDat$esan)}


[1] "reus"
              I          II        III         IV       gr1
II  3.551600e-04          NA         NA         NA        NA
III 5.368394e-08 0.004353211         NA         NA        NA
IV  1.467191e-10 0.004062074 0.56463777         NA        NA
gr1 5.171187e-03 0.361756594 0.34506796 0.53028646        NA
gr2 4.695635e-02 0.938276125 0.04076322 0.06934225 0.4114888

[1] "larg"
            I        II         III          IV       1av        2av       gr1
II  0.54203583        NA          NA          NA        NA         NA        NA
III 0.22233753 0.1002058          NA          NA        NA         NA        NA
IV  0.51062882 0.2202389 0.460456826          NA        NA         NA        NA
1av 0.10907085 0.3437917 0.011805059 0.029649001        NA         NA        NA
2av 0.02737454 0.1333209 0.002904960 0.005729777 0.5076156         NA        NA
gr1 0.71342935 0.4690854 0.687418537 1.000000000 0.1234447 0.03311754        NA
gr2 0.93137787 0.8694781 0.407991367 0.629225653 0.2594992 0.08075015 0.7588632




flegmatic 19.08.2010 14:19
(PS2004R @ 19.08.2010 13:39)
Ссылка на исходное сообщениеЯ предпочитаю применять бутсреп...

Можете прокоментировать?



PS2004R 19.08.2010 14:52
(flegmatic @ 19.08.2010 14:19)
Ссылка на исходное сообщение  Можете прокоментировать?


Как у Эфрона. Точечная оценка заменяется на densityplot() результата подачи на вход процедуры вычисления оценки массива перевыборок с возвращением исходных данных.



flegmatic 19.08.2010 15:38
(PS2004R @ 19.08.2010 14:52)
Ссылка на исходное сообщениеКак у Эфрона. Точечная оценка заменяется на densityplot() результата подачи на вход процедуры вычисления оценки массива перевыборок с возвращением исходных данных.

Переводу не поддаётся!!! smile.gif



PS2004R 19.08.2010 16:02
(flegmatic @ 19.08.2010 15:38)
Ссылка на исходное сообщение  Переводу не поддаётся!!!  smile.gif


http://free-books.dontexist.com/search?req...D&nametype=orig ( http://free-books.dontexist.com/search?req=%D1%8D%D1%84%D1%80%D0%BE%D0%BD&nametype=orig )

вот книжка на русском



flegmatic 20.08.2010 10:52
Возник следующий вопрос:
CODE
> print(names(myDat))
[1] "Конкретные.тесты."
[2] "Конкретные..тесты."
[3] "Конкретные...тесты."
(исправлено) Как заменить точки между словами: если одна точка, то на один пробел, а если две и более точек, то на " - " (пробел+дефис+пробел), а последнюю, она всегда одна, просто убрать?

P.S. Извиняюсь, больше исправлений не будет.



flegmatic 20.08.2010 12:12
Ещё один интересный (для меня smile.gif) вопрос:
Как указать R что уровни/подгрупы групирующей переменной должны обрабатываться и соответственно выводиться в определённой очерёдности, например: сперва "bacht" а потом "acht"? Как я понимаю, если не указывать очередь, то R сортирует уровни/подгрупы групирующей переменной по алфавиту или по величине, в зависимости от способа обозначения, так?



PS2004R 20.08.2010 12:28
(flegmatic @ 20.08.2010 12:12)
Ссылка на исходное сообщение  Ещё один интересный (для меня smile.gif) вопрос:
Как указать R что уровни/подгрупы групирующей переменной должны обрабатываться и соответственно выводиться в определённой очерёдности, например: сперва "bacht" а потом "acht"? Как я понимаю, если не указывать очередь, то R сортирует уровни/подгрупы групирующей переменной по алфавиту или по величине, в зависимости от способа обозначения, так?


А как же обещание себе прочитать руководство? smile.gif

?ordered



flegmatic 20.08.2010 12:32
(PS2004R @ 20.08.2010 12:28)
Ссылка на исходное сообщение  А как же обещание себе прочитать руководство? smile.gif

?ordered

Сдержу, но на это надо время.
Спасибо!



PS2004R 20.08.2010 12:34
(flegmatic @ 20.08.2010 10:52)
Ссылка на исходное сообщение  Возник следующий вопрос:
CODE
> print(names(myDat))
[1] "Конкретные.тесты."
[2] "Конкретные..тесты."
[3] "Конкретные...тесты."
(исправлено) Как заменить точки между словами: если одна точка, то на один пробел, а если две и более точек, то на " - " (пробел+дефис+пробел), а последнюю, она всегда одна, просто убрать?

P.S. Извиняюсь, больше исправлений не будет.


можно воспользоваться регулярными выражениями ?regex

надо разработать шаблон для функций g?sub() smile.gif



flegmatic 20.08.2010 13:51
(flegmatic @ 20.08.2010 10:52)
Ссылка на исходное сообщение  Возник следующий вопрос:
CODE
> print(names(myDat))
[1] "Конкретные.тесты."
[2] "Конкретные..тесты."
[3] "Конкретные...тесты."
Как заменить точки между словами: если одна точка, то на один пробел, если две и более точек, то на " - " (пробел+дефис+пробел), а последнюю, она всегда одна, просто убрать?

Решение (спасибо PS2004R):
CODE
names(myDat) <- gsub('\\.\\.\\.', ' - ', names(myDat))
names(myDat) <- gsub('\\.\\.', ' - ', names(myDat))
names(myDat) <- gsub('\\.', ' ', names(myDat))
names(myDat) <- sub(' +$', '', names(myDat))
names(myDat)


А с этим пока не справился frown.gif :
(PS2004R @ 20.08.2010 12:28)
Ссылка на исходное сообщение?ordered
Нужно найти решение для/на выше обсуждаемом примере (с pw функцией).



PS2004R 20.08.2010 14:36
(flegmatic @ 20.08.2010 13:51)
А с этим пока не справился  frown.gif :Нужно найти решение для/на выше обсуждаемом примере (с pw функцией).


"лучше день потерять, потом за пять минут долететь" (С) smile.gif

CODE

> sort(unique(data$esan))
[1] "1"    "2"    "3"    "4"    "5"    "6"    "7"    "acht"
> f <- ordered(data$esan, c("acht", "1", "2", "3", "4", "5", "6", "7" ))
> pw <- function (f,x) {data.na.omit <- na.omit(data.frame(f,x))
+                       print(pairwise.wilcox.test(as.numeric(data.na.omit[,1]),
+                                                  data.na.omit[,2],
+                                                  p.adjust.method = "fdr")$p.value)}
> pw(data$reus, f)
       acht            1          2         3        4
1 0.08804316           NA         NA        NA       NA
2 0.93827613 1.775800e-03         NA        NA       NA
3 0.08734976 4.026295e-07 0.01292797        NA       NA
4 0.11557041 2.200787e-09 0.01292797 0.6049690       NA
7 0.51436099 1.292797e-02 0.49330445 0.4933044 0.604969




flegmatic 20.08.2010 22:39
(PS2004R @ 20.08.2010 14:36)
Ссылка на исходное сообщение  "лучше день потерять, потом за пять минут долететь" (С) smile.gif
Стараюсь так и делать, просто дней свободных мало smile.gif. Спасибо!



PS2004R 27.08.2010 12:00
R признали программой года в конкурсе программ с открытым кодом.
http://www.infoworld.com/d/open-source/bos...=5#slideshowTop ( http://www.infoworld.com/d/open-source/bossie-awards-2010-the-best-open-source-application-development-software-140¤t=10&last=5#slideshowTop )



bubnilkin 31.08.2010 02:09
добрый день!

хочу проценты преобразовать с помощью Freeman-Tukey transformation [free(x, n).]

установил пакет 'binhf', но почему-то пишет: "Ошибка: не могу найти функцию "free"

и зачем n необходимо?


в справке (http://127.0.0.1:14218/library/stats/html/Binomial.html) для rbinom написано:

rbinom(n, size, prob)
n number of observations. {как я понимаю это n должно равняться объёму моей выборки}
size number of trials (zero or more). {зачем это необходимо?}
prob probability of success on each trial. {зачем это необходимо?}

заранее спасибо



PS2004R 31.08.2010 09:52
(bubnilkin @ 31.08.2010 02:09)
Ссылка на исходное сообщение  добрый день!

хочу проценты преобразовать с помощью Freeman-Tukey transformation [free(x, n).]

установил пакет 'binhf', но почему-то пишет: "Ошибка: не могу найти функцию "free"

и зачем n необходимо?
в справке (http://127.0.0.1:14218/library/stats/html/Binomial.html) для rbinom написано:

rbinom(n, size, prob)
n  number of observations. {как я понимаю это n должно равняться объёму моей выборки}
size  number of trials (zero or more). {зачем это необходимо?}
prob  probability of success on each trial.  {зачем это необходимо?}

заранее спасибо


Делаете ли Вы сначала library(binhf) ? Иначе free() не будет загружено.

rbinom() генерирует выборку заданного размера n.

Из ?rbinom следует:

The binomial distribution with ‘size’ = n and ‘prob’ = p has
density

p(x) = choose(n,x) p^x (1-p)^(n-x)

for x = 0, ..., n.


PS Если вдруг я не так понял вопрос:

size это размер серии независимых экспериментов с вероятностью успеха в каждом prob. А первый параметр в rbinom() это сколько случайных реализаций таких серий Вы заказываете. Каждая реализация будет представлена сверткой --- кол-во успешных независимых экспериментов данной серии.

PPS

If x is the numerator and n the denominator of a proportion, the Freeman-Tukey transformation (FT) is given by:

FT = arcsin( sqrt( x / n+1)) + arcsin( sqrt( x+1 / n+1))

, в случае процентов size <- 100 как я понимаю.



Guest 01.09.2010 10:04
спасибо за ответ smile.gif

т.е., как я понял, нужно всего лишь:
x (input data vector) - мои данные и n=100(%)

но вот пишу:

a <- c(read.delim("clipboard", header = FALSE))
b <- c(a)
d <- free(d, 100)

а выдаёт:

Ошибка в x/c : нечисловой аргумент для бинарного оператора

на я же перевёл a в b-вектор

???



PS2004R 01.09.2010 10:32
(Guest @ 01.09.2010 10:04)
Ссылка на исходное сообщение  спасибо за ответ smile.gif

т.е., как я понял, нужно всего лишь:
x (input data vector) - мои данные и n=100(%)

но вот пишу:

a <- c(read.delim("clipboard", header = FALSE))
b <- c(a)
d <- free(d, 100)

а выдаёт:

Ошибка в x/c : нечисловой аргумент для бинарного оператора

на я же перевёл a в b-вектор

???




This is a generic function which combines its arguments.

А для приведения используются функции типа as.к_чему_приводим()

Крайне рекомендую прочитать "An Introduction to R"



guest: bubnilkin 02.09.2010 12:16
вроде разобрался с тех. моментами. совет возьму на вооружение. спасибо

сделал вот так:

library("binhf")
a <- readClipboard()
b <- as.numeric(a)
d <- free(b, 100) ## т.к. вот как оказалось: https://stat.ethz.ch/pipermail/r-help/2006-...ber/114802.html ( https://stat.ethz.ch/pipermail/r-help/2006-October/114802.html )



PS2004R 02.09.2010 17:23
(guest: bubnilkin @ 02.09.2010 12:16)
Ссылка на исходное сообщение  вроде разобрался с тех. моментами. совет возьму на вооружение. спасибо

сделал вот так:

library("binhf")
a <- readClipboard()
b <- as.numeric(a)
d <- free(b, 100) ## т.к. вот как оказалось: https://stat.ethz.ch/pipermail/r-help/2006-...ber/114802.html ( https://stat.ethz.ch/pipermail/r-help/2006-October/114802.html )


лучше не загромождать промежуточными переменными, они только загромождает память
CODE

free(as.numeric(readClipboard()),
       100)




flegmatic 07.09.2010 18:53
Кто подскажет, как устанавливается пакет RcmdrPlugin.TextMining ( http://cran.r-project.org/web/packages/RcmdrPlugin.TextMining/index.html )? Сам не смог найти ответа.



PS2004R 07.09.2010 21:05
(flegmatic @ 07.09.2010 18:53)
Ссылка на исходное сообщение  Кто подскажет, как устанавливается пакет RcmdrPlugin.TextMining ( http://cran.r-project.org/web/packages/RcmdrPlugin.TextMining/index.html )? Сам не смог найти ответа.


install.packages(), в системе должен быть tcl/tk



daryd 10.09.2010 10:27
Здравствуйте!
Скажите пожалуйста, сталкивался ли кто-нибудь со следующей проблемой:
есть две последовательности (доза - ответ), они описываются экспоненциальной зависимостью вида y=a*exp(x*b). Нужно по этой зависимости предсказать значение y и его доверительный интервал при заданном значении x. Я знаю, что такая функция есть в пакете drc (dose-responce curves), но там за основу расчета взята логистическая функция. Другие зависимости этот пакет не использует.
Если такой расчет у при заданном х возможен, подскажите, как это сделать. Буду признательна за все направляющие пинки.



PS2004R 10.09.2010 11:00
(daryd @ 10.09.2010 10:27)
Ссылка на исходное сообщение  Здравствуйте!
Скажите пожалуйста, сталкивался ли кто-нибудь со следующей проблемой:
есть две последовательности (доза - ответ), они описываются экспоненциальной зависимостью вида y=a*exp(x*b). Нужно по этой зависимости предсказать значение y и его доверительный интервал при заданном значении x. Я знаю, что такая функция есть в пакете drc (dose-responce curves), но там за основу расчета взята логистическая функция. Другие зависимости этот пакет не использует.
Если такой расчет у при заданном х возможен, подскажите, как это сделать. Буду признательна за все направляющие пинки.


А остальные пакеты уже просмотрели?

bmd: Benchmark dose analysis for dose-response data

Benchmark dose analysis for continuous and quantal dose-response data.

DoseFinding: Planning and Analyzing Dose Finding experiments

The DoseFinding package provides functions for the design and analysis of dose-finding experiments (for example pharmaceutical Phase II clinical trials). It provides functions for: multiple contrast tests, fitting non-linear dose-response models, calculating optimal designs and an implementation of the MCPMod methodology. Currently only normally distributed homoscedastic endpoints are supported.

drfit: Dose-response data evaluation

drfit provides basic and easy-to-use functions for fitting dose-response curves to dose-response data, calculating some (eco)toxicological parameters and plotting the results. Functions that are fitted are the cumulative density function of the lognormal distribution (probit fit), of the logistic distribution (logit fit), of the weibull distribution (weibull fit) and a linear-logistic model ("linlogit" fit), derived from the latter, which is used to describe data showing stimulation at low doses (hormesis). In addition, functions checking, plotting and retrieving dose-response data retrieved from a database accessed via RODBC are included.

EffectiveDose: Estimation of the Effective Dose including Bootstrap confidence intervals

estimates the effective dose level for quantal bioassay data by nonparametric techniques and gives a bootstrap confidence interval



Darid 10.09.2010 11:59
Спасибо.
Отправляюсь смотреть.



ИгорьЧерниенко 27.09.2010 05:23
Здравствуйте!
Кто-нибудь может подсказать: имеется ли возможность наложить на карту, загруженную из шейп-файла, или в виде набора точек отрисованного линией, схему распределения построенную командой contourplot (levelplot) в lattice?
Заранее спасибо.



Guest 27.09.2010 05:27
Здравствуйте!
Кто-нибудь может подсказать: имеется ли возможность наложить на карту, загруженную из шейп-файла, или в виде набора точек отрисованного линией, схему распределения построенную командой contourplot (levelplot) в lattice?
Заранее спасибо.



Guest 27.09.2010 09:08
(ИгорьЧерниенко @ 27.09.2010 05:23)
Ссылка на исходное сообщение  Здравствуйте!
Кто-нибудь может подсказать: имеется ли возможность наложить на карту, загруженную из шейп-файла, или в виде набора точек отрисованного линией, схему распределения построенную командой contourplot (levelplot) в lattice?
Заранее спасибо.



например вот так:

CODE


> library(maps)
> coplot(lat ~ long | depth, data = quakes, number=4,
        panel=function(x, y, ...) {
          usr <- par("usr")
          rect(usr[1], usr[3], usr[2], usr[4], col="white")
          map("world2", regions=c("New Zealand", "Fiji"),
              add=TRUE, lwd=0.1, fill=TRUE, col="grey")
          text(180, -13, "Fiji", adj=1, cex=0.7)
          text(170, -35, "NZ", cex=0.7)
          points(x, y, pch=".")
        })


Paul Murrell, R Graphics

кроме points() можно что угодно вставлять что добавляет... например contour(, add = TRUE)



Darid 27.09.2010 14:49
Здравствуйте!
Уважаемые коллеги, подскажите, где можно найти R2 или другой критерий для оценки качества аппроксимации моделей в пакете drc.
Модели лог-логистическая и экспоненциальная из этого пакета.
Команда anova переадресует к команде modelFit, а та в свою очередь выдает ответ следующего содержания:
> mod<-drm(y~x,fct=LL.3())
> anova(mod)
Ошибка в anova.drc(mod) : Use the function modelFit()
> modelFit(mod)
No test available

ModelDf Log lik Df Chisq value p value

DRC model

Заранее спасибо.



Guest 27.09.2010 15:14
(Darid @ 27.09.2010 14:49)
Ссылка на исходное сообщение  

> mod<-drm(y~x,fct=LL.3())


а что summary(mod) показывает?



PS2004R 27.09.2010 15:27
(Guest @ 27.09.2010 15:14)
Ссылка на исходное сообщение  а что summary(mod) показывает?



сравнивать модели можно ещё с помощью mselect, похоже это и есть штатное средство выбора наилучшей модели.

Model selection by comparison of different models using the maximum log likelihood value, Akaike’s
information criterion (AIC), the estimated residual variance and the p-value from a lack-of-fit test
as criteria.


### Example with continuous data
## Fitting initial four-parameter log-logistic model
ryegrass.m1 <- drm(rootl ~ conc, data = ryegrass, fct = LL.4())
## Model selection
mselect(ryegrass.m1, list(LL.3(), LL.5(), W1.3(), W1.4(), W2.4(), baro5()))


а если anova(m1, m2) или anova.drc(m1,m2,) спросить?



Darid 28.09.2010 09:22
Спасибо!

summry(mod) показывает значения коэффициетнов в уравнении регрессии и их достоверность.

Функция mselect действительно выдает значения некоторых коэффициентов.
Предполагается, что наилучшей будет модель с наименьшим значением критерия AIC (так "считает" Википедия).
А что показывают logLik и Res var?
Правда ли то, что на основе этого теста следует выбрать экспоненциальную модель с двумя переменными (EXD.2)?

> mselect(mod1,list(LL.3(),EXD.3(),EXD.2()))
logLik AIC Lack of fit Res var
LL.3 -18.71964 45.43929 NA 21.54638
LL.3 -18.71964 45.43929 NA 21.54638
EXD.3 -19.41176 46.82351 NA 26.25758
EXD.2 -19.56717 45.13433 NA 21.95981



PS2004R 28.09.2010 10:04
(Darid @ 28.09.2010 09:24)
Ссылка на исходное сообщение  Спасибо!

summry(mod) показывает значения коэффициетнов в уравнении регрессии и их достоверность.

Функция mselect действительно выдает значения некоторых коэффициентов.
Предполагается, что наилучшей будет модель с наименьшим значением критерия AIC (так "считает" Википедия).
А что показывают logLik и Res var?
Правда ли то, что на основе этого теста следует выбрать экспоненциальную модель с двумя переменными (EXD.2)?

> mselect(mod1,list(LL.3(),EXD.3(),EXD.2()))
         logLik      AIC Lack of fit  Res var
LL.3  -18.71964 45.43929          NA 21.54638
LL.3  -18.71964 45.43929          NA 21.54638
EXD.3 -19.41176 46.82351          NA 26.25758
EXD.2 -19.56717 45.13433          NA 21.95981


Лучшая модель минимальной сложности дорогого стоит smile.gif Чем больше параметров, тем лучше подгонка модели, но тем меньше её прогностическая сила. Предсказание для точки не участвовавшей в подгонке модели в идеале должно иметь минимальную ошибку.

Чем больше logLik тем лучше подогнана модель, но чем больше в ней параметров тем хуже. AIC вводит единый критерий выбора модели.

Для выбора, по моему, надо обязательно сначала посмотреть:

plot(model.exd2)

plot(model.exd3, add = TRUE)

plot(model.ll3, add = TRUE)

Также можно посмотреть статистику при последовательном исключении одной точки из выборки. Смотреть надо как распределится AIK, вдруг он от одной единственной точки зависит?

Для исключённых точек можно смотреть
PR(model, вектор_неучавствовавших_в_текущем_расчете)

по результатам на пример: 1) можно исключить подозрительные точки, или 2) выбрать модель которая удовлетворительнее дает прогноз для исключённых точек.

Неплохо проделать то, что рекомендует руководство, для определения достоверности отличий моделей

modelFit(model.exd3)

modelFit(model.exd2)

anova(model.exd2,model.exd3)

... как то так smile.gif



Darid 28.09.2010 10:24
Большое спасибо.
Сейчас попробую.



Guest 28.09.2010 11:04
...кроме points() можно что угодно вставлять что добавляет... например contour(, add = TRUE)
[/quote]
Спасибо огромное, но это не совсем мой случай.
Имеются интерполированные в ячейки регулярной сетки при помощи gstat значения определенного параметра. Необходимо отобразить их в виде схемы распределения с привязкой к схеме побережья. Функция contour отрисовывать мою сетку категорически не хочет, зато она прекрасно рисуется contourplot и levelplot из lattice.
Карта побережья загружена из шейпа при помощи maptools, может быть представлена в виде полигонов или преобразована в точки - узлы полигонов, которые можно нарисовать при помощи plot(oj,type='l') или lines. Хотелось бы эти карты совместить. При использовании функции lines схема побережья на карту накладывается, но очертания побережья и схема распределения не имеют между собой ничего общего.

Пользуясь случаем, еще вопрос, для тех кто использует gstat.
Имеется код:
создаем объект:
cmm.gs<-gstat(id=names(cmm)[4],formula=cmm[,4]~1,locations=lon+lat,data=cmm)
(cmm - таблица данных со значениями некоторых параметров на станциях)
далее, добавляем данные в объект:
for(i in c(5:12)) {
cmm.gs<-gstat(cmm.gs,id=names(cmm)[i],formula=cmm[,i]~1,locations=lon+lat,data=cmm)
}
имена переменных в объект записываются как положено, а вот данные - только из одного столбца. Может, у меня в цикле что-то не так?
Заранее спасибо.



PS2004R 28.09.2010 11:35
[quote=Guest,28.09.2010 11:04]Ссылка на исходное сообщение  ...кроме points() можно что угодно вставлять что добавляет... например contour(, add = TRUE)
[/quote]
Спасибо огромное, но это не совсем мой случай.
Имеются интерполированные в ячейки регулярной сетки при помощи gstat значения определенного параметра. Необходимо отобразить их в виде схемы распределения с привязкой к схеме побережья. Функция contour отрисовывать мою сетку категорически не хочет, зато она прекрасно рисуется contourplot и levelplot из lattice.
Карта побережья загружена из шейпа при помощи maptools, может быть представлена в виде полигонов или преобразована в точки - узлы полигонов, которые можно нарисовать при помощи plot(oj,type='l') или lines. Хотелось бы эти карты совместить. При использовании функции lines схема побережья на карту накладывается, но очертания побережья и схема распределения не имеют между собой ничего общего.

[/quote]

что тогда мешает в levelplot вставить map?

PS тогда еще image можно попробовать, а что конкретное говорит contour?



PS2004R 28.09.2010 12:04
(Guest @ 28.09.2010 11:04)
Пользуясь случаем, еще вопрос, для тех кто использует gstat.
Имеется код:
создаем объект:
cmm.gs<-gstat(id=names(cmm)[4],formula=cmm[,4]~1,locations=lon+lat,data=cmm)
(cmm - таблица данных со значениями некоторых параметров на станциях)
далее, добавляем данные в объект:
for(i in c(5:12)) {
cmm.gs<-gstat(cmm.gs,id=names(cmm)[i],formula=cmm[,i]~1,locations=lon+lat,data=cmm)
}
имена переменных в объект записываются как положено, а вот данные - только из одного столбца. Может, у меня в цикле что-то не так?
Заранее спасибо.


я qstat не использую, хотя меня он и интересует smile.gif

1) достаточно 5:12 написать

2) почему Вы в cmm.gs пытаетесь поместить несколько объектов?

его значение это "Value an object of class gstat, which inherits from list". Это объект на основе списка. Логично будет получить несколько объектов класса qstat размещая их в общем списке.

CODE

rezult.list <- lapply(5:12, function (i) gstat(cmm.gs,id=names(cmm)[i], formula=cmm[,i]~1,locations=lon+lat,data=cmm) )



и обращаться соответственно rezult.list[[j]]$(что то из возвращаемых qstat)

3) можно и сразу извлекать результат из объекта возвращаемого qstat, и уже его объединять в нужную Вам структуру данных...

например qstat()$имя_параметра


data --- list; each element is a list with the formula, locations, data, nvars, beta, etc., for a variable

qstat()$data[[1]] это первая переменная по номеру, а qstat()$data[[1]]$beta конкретный параметр первой переменной. Проверить увы не могу...


4) Лучше добавить в FAQ необходимость в вопросе приводить полностью функциональный пример. Так значительно интереснее отвечать будет smile.gif



Guest 28.09.2010 15:48
[/quote]

что тогда мешает в levelplot вставить map?
как это? add=T?
PS тогда еще image можно попробовать, а что конкретное говорит contour?
[/quote]
контур требует, чтобы х и у были отсортированы по возрастанию. Я так понял, что и х и у должны равномерно возрастать, в моем случае это невозможно. Имаж, в сущности, тот же контур:о)



Guest 28.09.2010 15:54
1) достаточно 5:12 написать
спасибо :о)
2) почему Вы в cmm.gs пытаетесь поместить несколько объектов?

В gstat такой синтаксис - при добавлении поля в объект, если там уже есть другие поля, он указывается. Несколько полей нужны для расчета кокригинга :о)

3) можно и сразу извлекать результат из объекта возвращаемого qstat, и уже его объединять в нужную Вам структуру данных...

например qstat()$имя_параметра


data --- list; each element is a list with the formula, locations, data, nvars, beta, etc., for a variable

qstat()$data[[1]] это первая переменная по номеру, а qstat()$data[[1]]$beta конкретный параметр первой переменной. Проверить увы не могу...
О, с этим поразбираюсь
4) Лучше добавить в FAQ необходимость в вопросе приводить полностью функциональный пример. Так значительно интереснее отвечать будет smile.gif
[/quote]
Постараюсь учесть :0)|||
Спасибо за ответы!



PS2004R 28.09.2010 21:17
(Guest @ 28.09.2010 15:48)
Ссылка на исходное сообщение  


что тогда мешает в levelplot вставить map?
как это? add=T?
PS тогда еще image можно попробовать, а что конкретное говорит contour?

контур требует, чтобы х и у были отсортированы по возрастанию. Я так понял, что и х и у должны равномерно возрастать, в моем случае это невозможно. Имаж, в сущности, тот же контур:о)


Вот текст в пример засунутый

CODE

levelplot(rnorm(10) ~ 1:10 + sort(runif(10)),
            panel = function(x, y, z, ...) {panel.levelplot(x,y,z, ...);
                                                      ltext(5, 0.5, "Fiji")})


по идее должен и map залезать в панель

CODE

         .......        
         map("world2", regions=c("New Zealand", "Fiji"),
             add=TRUE, lwd=0.1, fill=TRUE, col="grey")
         .......



А по контуру лучше на конкретном примере смотреть всё же...

PS есть "родной" способ выводить карты в пакете latticeExtra



PS2004R 28.09.2010 22:27
(Guest @ 28.09.2010 15:54)
Ссылка на исходное сообщение  
2) почему Вы в cmm.gs пытаетесь поместить несколько объектов?

В gstat такой синтаксис - при добавлении поля в объект, если там уже есть другие поля, он указывается. Несколько полей нужны для расчета кокригинга :о)



# multivariable; thanks to M. Rufino:
meuse.g <- gstat(id = "zn", formula = log(zinc) ~ 1, data = meuse)
meuse.g <- gstat(meuse.g, "cu", log(copper) ~ 1, meuse)
meuse.g <- gstat(meuse.g, model = vgm(1, "Sph", 900, 1), fill.all = TRUE)
x <- variogram(meuse.g, cutoff = 1000)

а оно не требует именно такой порядок?

в цикле писать

for(i in c(5:12)) {
cmm.gs<-gstat(cmm.gs, names(cmm)[i], cmm[,i]~1, locations=lon+lat, cmm)
}

или "руками" цикл проходит корректно?



Guest 30.09.2010 02:10
В том то все и дело, что руками цикл корректно проходит



PS2004R 30.09.2010 07:53
(Guest @ 30.09.2010 02:10)
Ссылка на исходное сообщение  В том то все и дело, что руками цикл корректно проходит


это с окружением проблема...

Сначала сделать call(), потом его исполнить eval().

Без рабочего примера ничего попробовать не могу frown.gif.

PS можно почитать еще ?force , я for() вообще стараюсь не использовать... Дайте полный пример это будет интересно и позновательно не только мне [я считаю] smile.gif

PPS как то так надо попробовать [слава фезаму]

for(i in c(5:12)) {
cmm.gs<-gstat(force(cmm.gs), names(cmm)[i], cmm[,i]~1, locations=lon+lat, cmm)
}

или

for(i in c(5:12)) {
cmm.gs.tmp<-gstat(cmm.gs, names(cmm)[i], cmm[,i]~1, locations=lon+lat, cmm)
cmm.gs<-cmm.gs.tmp
}

или

for(i in c(5:12)) {
cmm.gs.tmp<-gstat(cmm.gs, names(cmm)[i], cmm[,i]~1, locations=lon+lat, cmm)
cmm.gs<-print(cmm.gs.tmp)
}

или ....



Darid 04.10.2010 12:50
Скажите пожалуйста!
Может ли критерий Акаике принимать отрицательные значения?

Спасибо.



PS2004R 04.10.2010 13:07
(Darid @ 04.10.2010 12:50)
Ссылка на исходное сообщение  Скажите пожалуйста!
Может ли критерий Акаике принимать отрицательные значения?

Спасибо.


Да.


It is a common misconception that the log-likelihood must be negative. If the likelihood is derived from a probability density it can quite reasonably exceed 1 which means that log-likelihood is positive, hence the deviance and the AIC are negative.



http://tolstoy.newcastle.edu.au/R/help/05/04/2614.html ( http://tolstoy.newcastle.edu.au/R/help/05/04/2614.html )



Lea1 07.10.2010 11:11
Добрый день,

У меня проблема с загрузкой sph file.
Мои шаги:

включаю пакеты maptools и gpclib,

v <-readShapePoly("D:/sph /G3G02.shp")
Ошибка в getinfo.shape(filen) : Error opening SHP file

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

Не могу понять, где ошибка. ((



PS2004R 07.10.2010 11:28
(Lea1 @ 07.10.2010 11:11)
Ссылка на исходное сообщение 

Не могу понять, где ошибка. ((


file name that you supply should not include the extension

?readShapePoly
CODE

fn: shapefile layer name, when writing omitting the extensions
         *.shp, *.shx and *.dbf, which are added in the function




Lea1 07.10.2010 11:41
Спасибо за ответ.

Попробовала
> v <-readShapePoly("D:/sph /G3G02")
Ошибка в getinfo.shape(filen) : Error opening SHP file

До этого работала с другим шейп файлом(3 месяца назад все было ок, а сейчас

> v <-readShapePoly("D:/sph/G3G01")
Field name: >;51 changed to: X...51
Ошибка в read.shape(filen = fn, verbose = verbose, repair = repair) :
File size and implied file size differ, consider trying repair=TRUE

Что это могло бы значить?))



Lea1 07.10.2010 11:58
Нашла ошибку, использовала для распаковки архиватор winrar.
Сейчас вопользовалась другим архиватором 7- Zip, шейп файл загружается!

Спасибо за внимание к моей проблеме и терпение к новичкам в R.



Игорь Черниенко 08.10.2010 03:13
(Guest @ 28.09.2010 11:04)
Ссылка на исходное сообщение  ...кроме points() можно что угодно вставлять что добавляет... например contour(, add = TRUE)

Спасибо огромное, но это не совсем мой случай.
Имеются интерполированные в ячейки регулярной сетки при помощи gstat значения определенного параметра.

С этим разобрался - просто преобразую результат в SpatialGridDataFrame и рисую функцией image, к нему lines добавляется нормально
С циклом разберусь попозже. Огромное спасибо за советы!



Lea1 08.10.2010 18:53
Добрый вечер!

Подскажите, пожалуйста, что значит команда

idd[2895]<-2369
idd[2139]<-1877

и
v[2895,]<-v[2369,]
v[2139]<-v[1877]
Заместить или, что то другое.

v<-v[-2643,] - удалить?



PS2004R 08.10.2010 21:19
(Lea1 @ 08.10.2010 18:53)
Ссылка на исходное сообщение  Добрый вечер!

Подскажите, пожалуйста, что значит команда

idd[2895]<-2369
idd[2139]<-1877

и
v[2895,]<-v[2369,]
v[2139]<-v[1877]
Заместить или, что то другое.

v<-v[-2643,] - удалить?



x[i] iтый элемент вектора x

x[i] <- z поместить в iтый элемент вектора значение z

x[i,] iтая строка таблицы x (в том числе многомерной таблицы)

соответственно x[i,]<- x[z,] поместить z строку вместо i строки

лучше всего прочитать "введение в R", вот ссылка на русском http://m7876.wiki.zoho.com/Introduction-to...912000000002015 ( http://m7876.wiki.zoho.com/Introduction-to-R.html?pid=78912000000002015 )

v[-i] за исключением iтого элемента вектора



Lea1 08.10.2010 22:03
Спасибо!

так как я изучала и работаю со статистической программой STATA, я пытаюсь найти аналог команды R в STATA., когда это возможно, конечно.

К сожалению, для меня)), STATA не обеспечивает все необходимые манипуляции с данными, приходится работать с R.



PS2004R 09.10.2010 18:31
(Lea1 @ 08.10.2010 22:03)
Ссылка на исходное сообщение  приходится работать с R.


Прочитайте главу про анализ тормозного пути машин, R среда во многом появилась из APL систем. Очень краткое и полное представление на русском языке что такое техника прямой работы с R.

http://flibusta.net/b/156597 ( http://flibusta.net/b/156597 )


Вторая техника, это "литературное программирование". В случае R это погружение кода на R в среду документа на LaTeX. Для обеспечения этого служит Sweave, ну и редактор с подходящими режимами. На выходе процесса получается оформленный отчёт или статья, или презентация. Причём при достаточной сноровки в LaTeX из одного и того же документа smile.gif



Darid 12.10.2010 14:12
Здраствуйте!
Уважаемые коллеги, подскажите, есть ли в R пакет, позволяющий рассчитать распределение энергии фотонов и электронов в объекте с помощью метода Монте-Карло (MCNP)?



PS2004R 12.10.2010 16:01
(Darid @ 12.10.2010 14:12)
Ссылка на исходное сообщение  Здраствуйте!
Уважаемые коллеги, подскажите, есть ли в R пакет, позволяющий рассчитать распределение энергии фотонов и электронов в объекте с помощью метода Монте-Карло  (MCNP)?


CODE

Программы MCNP/MCNPX также могут использоваться для расчетов наработки различных ядерных материалов и попадают под экспортные ограничения Министерства энергетики США, а доступ к этим программам (платный в общем случае) предоставляется по выполнению ряда требований.

Программа может предоставляться бесплатно сотрудникам Министерства энергетики США (и ряду отдельных категорий исследователей).



террористо?



Darid 13.10.2010 14:38
Скорее длячервядозурассчетисто.
Ух ты! Как все оказывается круто.
Эту методу использовали для расчета коэффициентов конверсии дозы для оценки дозовых нагрузок на растения и животных. Наше животное не совсем подходит под их референтные виды и мы хотели сами прпробовать. Жаль.



Guest 14.10.2010 15:36
http://www-rsicc.ornl.gov/rsiccnew/CodesAv...leElsewhere.htm ( http://www-rsicc.ornl.gov/rsiccnew/CodesAvailableElsewhere.htm )



Игорь Черниенко 25.10.2010 14:24
Здравствуйте!
Подскажите, можно ли из нескольких объектов-результатов выполнения cor.test вытащить коэффициент корреляции, доверительные интервалы и уровень значимости и свести их в таблицу вида r | low95 |up95 | p ?
__|______|____|___
Заранее благодарен



PS2004R 25.10.2010 15:45
(Игорь Черниенко @ 25.10.2010 14:24)
Ссылка на исходное сообщение  Здравствуйте!
Подскажите, можно ли из нескольких объектов-результатов выполнения cor.test вытащить коэффициент корреляции, доверительные интервалы и уровень значимости и свести их в таблицу вида r | low95 |up95 | p  ?
                                      __|______|____|___
Заранее благодарен


CODE


> x <- c(44.4, 45.9, 41.9, 53.3, 44.7, 44.1, 50.7, 45.2, 60.1)
> y <- c( 2.6,  3.1,  2.5,  5.0,  3.6,  4.0,  5.2,  2.8,  3.8)
> ct.ob <- cor.test(x, y)
> ct.ob$
ct.ob$statistic    ct.ob$p.value      ct.ob$null.value   ct.ob$method       ct.ob$conf.int    
ct.ob$parameter    ct.ob$estimate     ct.ob$alternative  ct.ob$data.name    
> ct.ob$conf.int
[1] -0.1497426  0.8955795
attr(,"conf.level")
[1] 0.95
> ct.ob$statistic
      t
1.841083
> ct.ob$estimate
     cor
0.5711816
> ct.ob$conf.int[[1]]
[1] -0.1497426
> ct.ob$conf.int[[2]]
[1] 0.8955795

> c(ct.ob$estimate[[1]],ct.ob$conf.int[[1]],ct.ob$conf.int[[2]])
[1]  0.5711816 -0.1497426  0.8955795



Как то так. Приведите пример как пытаетесь считать, в виде чего хранятся данные, тогда можно помочь подробнее...



Guest 25.10.2010 22:30
Большое спасибо! Это именно то что нужно, дальше вставлю этот код в цикл. Я еще не до конца освоился с синтаксисом. :о)



PS2004R 26.10.2010 13:45
(Guest @ 25.10.2010 22:30)
Ссылка на исходное сообщение  Большое спасибо! Это именно то что нужно, дальше вставлю этот код в цикл. Я еще не до конца освоился с синтаксисом. :о)


"какой мэханизм? "(С) то есть "цикл" smile.gif

Цикл в R "приют юдоли и печали" smile.gif

Скорее надо вставить в sapply()



Игорь Черниенко 26.10.2010 20:28
"Или так" (С) :о)|||



Игорь Черниенко 29.10.2010 06:54
Здравствуйте!
Имеется набор данных dth:
group ye n est
g1 2007 11.0 11.0
g1 2008 10.0 9.0
g1 2009 6.1 7.3
g1 2010 6.5 6.0
g2 2007 6.0 6.0
g2 2008 4.1 3.6
g2 2009 1.5 2.2
g2 2010 1.51.33878096089058
можно ли в lattice построить вот такой график
yplot(est + n ~ ye | group, pch=16,
scales=list(x=list(relation='same'),
y=list(relation='free')),
data=dth)
но чтобы один ряд был точечным а второй - линией?
заранее спасибо



Игорь Черниенко 29.10.2010 10:05
Нашел. Вот как надо, если кто не знает:
yplot(est + n ~ ye | group, pch=16,
scales=list(x=list(relation='same'),
y=list(relation='free')), type=c('l','p'), distribute.type=T,
data=dth)
Нашел тут книжку:
http://www.accuratefiles.com/fileinfo/gs1a50462h1i0# ( http://www.accuratefiles.com/fileinfo/gs1a50462h1i0# )
Только не надо нажимать на большую кнопку "Скачать" , там есть маленькая кнопка "download" :0)|||



PS2004R 29.10.2010 10:08
(Игорь Черниенко @ 29.10.2010 06:54)
Ссылка на исходное сообщение  Здравствуйте!
Имеется набор данных dth:
group ye n est
g1 2007 11.0 11.0
g1 2008 10.0 9.0
g1 2009 6.1 7.3
g1 2010 6.5 6.0
g2 2007 6.0 6.0
g2 2008 4.1 3.6
g2 2009 1.5 2.2
g2 2010 1.51.33878096089058
можно ли в lattice построить вот такой график
yplot(est + n ~ ye |  group, pch=16,
scales=list(x=list(relation='same'),
y=list(relation='free')),
data=dth)
но чтобы один ряд был точечным а второй - линией?
заранее спасибо



CODE


> yplot(est + n ~ ye | group, pch=16,
+ scales=list(x=list(relation='same'),
+ y=list(relation='free')),
+ data=dth)
Ошибка: не могу найти функцию "yplot"





PS2004R 29.10.2010 10:24
(Игорь Черниенко @ 29.10.2010 10:05)
Ссылка на исходное сообщение  Нашел. Вот как надо, если кто не знает:
yplot(est + n ~ ye | group, pch=16,
scales=list(x=list(relation='same'),
y=list(relation='free')), type=c('l','p'), distribute.type=T,
data=dth)
Нашел тут книжку:
http://www.accuratefiles.com/fileinfo/gs1a50462h1i0# ( http://www.accuratefiles.com/fileinfo/gs1a50462h1i0# )
Только не надо нажимать на большую кнопку "Скачать" , там есть маленькая кнопка "download" :0)|||




File removed...

что такое yplot?



Игорь Черниенко 29.10.2010 11:07
1. вкралась досадная опечатка: это xyplot

2. А у меня ссылка работает. Там книга Lattice Multivariate Data Visualization with R, страниц на 300



Guest 29.10.2010 11:23
и еще: в последней строке данных, которые я привел вместо g2 2010 1.51.33878096089058 должно быть g2 2010 1.5 1.34



PS2004R 29.10.2010 11:50
(Игорь Черниенко @ 29.10.2010 11:07)
Ссылка на исходное сообщение  
2. А у меня ссылка работает. Там книга Lattice Multivariate Data Visualization with R, страниц на 300


да, отключенный ява скрипт сыграл свою роль... спасибо за ссылку!



PS2004R 29.10.2010 12:06
(Guest @ 29.10.2010 11:23)
Ссылка на исходное сообщение  и еще: в последней строке данных, которые я привел вместо g2 2010 1.51.33878096089058 должно быть g2 2010 1.5 1.34


это я сразу исправил smile.gif

код лучше форматировать как s выражения

CODE


xyplot(est + n ~ ye | group,
         pch=16,
         scales=list(x=list(relation='same'),
                          y=list(relation='free')),
         type=c('l','p'),
         distribute.type=T,
         data=dth)





Darid 12.11.2010 09:49
Здравствуйте!
Скажите пожалуйста, можно ли в R не просто анализировать готовые карты, а "рисовать" свои?
Например, есть образцы почв. Для каждого есть координаты точки отбора и определено содержание некоторых элементов. Необходимо нарисовать карты площадного распределения элементов. (Я знаю, что это умеет делать ArcView, а в R получилось только раскрасить точки в разые цвета в зависимости от концентрации.) Можно ли получить картинку площадного распределения? Если, да - какой пакет для этого используют?
Спасибо.



PS2004R 12.11.2010 10:25
(Darid @ 12.11.2010 09:49)
Ссылка на исходное сообщение  Здравствуйте!
Скажите пожалуйста, можно ли в R не просто анализировать готовые карты, а "рисовать" свои?
Например, есть образцы почв. Для каждого есть координаты  точки отбора и определено содержание некоторых элементов. Необходимо нарисовать карты площадного распределения элементов. (Я знаю, что это умеет делать ArcView, а в R получилось только раскрасить точки в разые цвета в зависимости от концентрации.) Можно ли получить картинку площадного распределения? Если, да - какой пакет для этого используют?
Спасибо.


R рисует изолинии даже в базовой комплектации, это один из стандартных типов графиков.

?contour
?contourplot

и понадобится 2d сглаживание, оно есть в пакетах...

PS вот так contour(density2d(x,y))

а, еще z координата... ну тогда ?panel.voronoi , оно сегментирует автоматически и раскрашивает согласно z

CODE

   These panel functions for ‘levelplot’ can represent irregular (x,
    y) points with a color covariate.  ‘panel.levelplot.points’ simply
    draws color-coded points.  ‘panel.voronoi’ uses the ‘deldir’
    package to calculate the spatial extension of a set of points in 2
    dimensions.  This is known variously as a Voronoi mosaic, a
    Dirichlet tesselation, or Thiessen polygons.


PPS
есть еще вот это
CODE

fields::Krig

Fits a surface to irregularly spaced data. The Kriging model assumes that the unknown function is
a realization of a Gaussian random spatial processes. The assumed model is additive Y = P(x) +
Z(X) + e, where P is a low order polynomial and Z is a mean zero, Gaussian stochastic process with
a covariance that is unknown up to a scale constant. The main advantages of this function are the
flexibility in specifying the covariance as an R language function and also the supporting functions
plot, predict, predict.se, surface for subsequent analysis. Krig also supports a correlation model
where the mean and marginal variances are supplied.




Darid 13.11.2010 12:44
Спасибо большое.
Попробую.



Игорь Черниенко 01.12.2010 04:52
Я для подобных целей gstat использую. Там, кстати, рабочий пример есть как раз про почвы.





Powered by Invision Power Board (http://www.invisionboard.com)
© Invision Power Services (http://www.invisionpower.com)