Molbiol.ru | О проекте | Справочник | Методы | Растворы | Расчёты | Литература | Орг.вопросы Web | Фирмы | Coffee break | Картинки | Работы и услуги | Биржа труда | Zbio-wiki NG SEQUENCING · ЖИЗНЬ РАСТЕНИЙ · БИОХИМИЯ · ГОРОДСКИЕ КОМАРЫ · А.А.ЛЮБИЩЕВ · ЗООМУЗЕЙ Темы за 24 часа [ Вход* | Регистрация* ] Форум: | |
ИНО Постоянный участник Донецк |
|
ИНО Постоянный участник Донецк |
CODE df3$Resp<-as.factor(df3$Resp) #Способ рекомендованный Щитиковым и Мастицким: resp_<-cbind(dataframe1$RespA, dataframe1$RespB) model1<-glm(resp_~Pred2, family = binomial("logit"), data = dataframe1) #Способ, который я встречал в любых других публикациях и всегда применял ранее: model2<-glm(Resp~Pred2, family = binomial("logit"), data = df3) Понятно, что на этом высосанном из пальца наборе данных обе модели одинаково бессмысленны, но проблема в том, что почему-то они разные! Какая из них задана верно? Сообщение было отредактировано ИНО - 10.01.2021 05:28 |
plantago Постоянный участник |
?glm: "... A typical predictor has the form ‘response ~ terms’ ... For ‘binomial’ and ‘quasibinomial’ families the response can also be specified as ... a two-column matrix with the columns giving the numbers of successes and failures. > А не могли бы Вы завернуть этот код в одну функцию Не думаю, что это разумно, потому что данные могут быть очень разными. Проще понять способ, записать себе его с комментариями и применять, модифицируя.
|
ИНО Постоянный участник Донецк |
Свои таблицы я уже по разворачивал с использованием кода, предложенного гостем (вашего ответа еще не было тогда), но теперь уже стал сомневаться: а то ли я по ним считаю, потому как по свернутым и поданным процитированным Вами способом результат другой CODE dataframe1<-data.frame(Pred1=c(1,2,3), Pred2=c(30, 20, 10), RespA=c(0, 1, 2), RespB=c(2, 1, 0) ) resp_<-cbind(dataframe1$RespB, dataframe1$RespA) model1<-glm(resp_~Pred2, family = binomial("logit"), data = dataframe1) summary(model1) Call: glm(formula = resp_ ~ Pred2, family = binomial("logit"), data = dataframe1) Deviance Residuals: 1 2 3 1.29e-05 0.00e+00 -1.29e-05 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -47.80 94013.91 -0.001 1 Pred2 2.39 4700.70 0.001 1 (Dispersion parameter for binomial family taken to be 1) Null deviance: 5.5452e+00 on 2 degrees of freedom Residual deviance: 3.3297e-10 on 1 degrees of freedom AIC: 5.3863 Number of Fisher Scoring iterations: 22 CODE dataframe2<-data.frame(Pred1=c(1,1,2,2,3,3), Pred2=c(30,30,20,20,10,10), Resp=factor(c("B","B","A","B","A","A"))) model2<-glm(Resp~Pred2, family = binomial("logit"), data = dataframe2) summary(model2) Call: glm(formula = Resp ~ Pred2, family = binomial("logit"), data = dataframe2) Deviance Residuals: 1 2 3 4 5 6 0.00005 0.00005 -1.17741 1.17741 -0.00005 -0.00005 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -41.132 17730.370 -0.002 0.998 Pred2 2.057 886.518 0.002 0.998 (Dispersion parameter for binomial family taken to be 1) Null deviance: 8.3178 on 5 degrees of freedom Residual deviance: 2.7726 on 4 degrees of freedom AIC: 6.7726 Number of Fisher Scoring iterations: 19 Сообщение было отредактировано ИНО - 10.01.2021 16:59 |
ПолинаШ Участник |
АВТОР <- c("Аветисян Самвел", "Самуил Маршак") search <- "Сам" tibble(АВТОР) %>% filter(str_detect(АВТОР, pattern = paste("^",search,sep = ""))) # A tibble: 1 x 1 АВТОР <chr> 1 Самуил Маршак tibble(АВТОР) %>% filter(str_detect(АВТОР, pattern = paste("\<", search,sep = ""))) Ошибка: '\<' -- нераспознаваемый префикс в строке в строке начинающейся с ""\<" tibble(АВТОР) %>% filter(str_detect(АВТОР, pattern = paste("\\<", search,sep = ""))) # A tibble: 0 x 1 # ... with 1 variable: АВТОР <chr> И простите за ламерский вопрос: у меня перестали работать все иконки панели редактора сообщений кроме кнопок "Пред.просмотр" и "Ответить" |
Guest IP-штамп: frD0EXpX.eS.Q гость |
(ПолинаШ @ 13.01.2021 12:27) Как сделать выборку строк из таблицы с полем char, содержащим набор слов. Поиск должен идти с начала каждого слова. В примере ниже должны выбраться обе строки: АВТОР <- c("Аветисян Самвел", "Самуил Маршак") search <- "Сам" tibble(АВТОР) %>% filter(str_detect(АВТОР, pattern = paste("^",search,sep = ""))) # A tibble: 1 x 1 АВТОР <chr> 1 Самуил Маршак tibble(АВТОР) %>% filter(str_detect(АВТОР, pattern = paste("\<", search,sep = ""))) Ошибка: '\<' -- нераспознаваемый префикс в строке в строке начинающейся с ""\<" tibble(АВТОР) %>% filter(str_detect(АВТОР, pattern = paste("\\<", search,sep = ""))) # A tibble: 0 x 1 # ... with 1 variable: АВТОР <chr> И простите за ламерский вопрос: у меня перестали работать все иконки панели редактора сообщений кроме кнопок "Пред.просмотр" и "Ответить" Можно вот так: CODE > lapply(lapply(strsplit(АВТОР, " "), grepl, pattern="^Сам"), any) [[1]] [1] TRUE [[2]] [1] TRUE
|
ИНО Постоянный участник Донецк |
Я же развернул dreissena_infection.txt в длинный вид и задал модель по нему. Коэффициенты в точности сошлись с посчитанными по таблицы частот, уровни значимости предикторов тоже сошлись, уровень значимости свободного члена сошелся не совсем точно, степени свободы, девианс и AIC совсем не сошлись. CODE #Модель, построенная по таблице частот: summary(M1) Call: glm(formula = inf.tbl ~ Day + Depth, family = binomial(link = "logit"), data = infection) Deviance Residuals: Min 1Q Median 3Q Max -2.1129 -0.9595 -0.1563 0.7182 2.0214 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -3.836627 0.651911 -5.885 3.98e-09 *** Day 0.011039 0.004623 2.388 0.01695 * Depth4m 1.543597 0.537679 2.871 0.00409 ** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for binomial family taken to be 1) Null deviance: 34.394 on 11 degrees of freedom Residual deviance: 18.146 on 9 degrees of freedom AIC: 45.338 Number of Fisher Scoring iterations: 5 #Модель, построенная по таблице, развернутой в длинный вид: summary(M3) Call: glm(formula = variable ~ Day + Depth, family = binomial(link = "logit"), data = infection3) Deviance Residuals: Min 1Q Median 3Q Max -0.9297 -0.5024 -0.4034 -0.2380 2.7778 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -3.836627 0.651858 -5.886 3.96e-09 *** Day 0.011039 0.004623 2.388 0.01695 * Depth4m 1.543597 0.537638 2.871 0.00409 ** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for binomial family taken to be 1) Null deviance: 138.38 on 182 degrees of freedom Residual deviance: 122.14 on 180 degrees of freedom AIC: 128.14 Number of Fisher Scoring iterations: 5 Также совершенно иные остатки, их количество равно не количеству наблюдений, а количеству строк. Остаток для каждой строки близок (но не равен) среднему арифметическому из остатков для содержащихся в этой строке наблюдений, посчитанным по длинному формату. Как использовать такие остатки в диагностике модели совсем непонятно, ведь в случаях если в одной строке встречались оба варианта отклика, получается некая неинформативная "средняя температура по больнице". predict() c опцией response=TRUE возвращает абсолютно тождественные результаты для обоих моделей. Вывод: задание таблицы частот и таблицы наблюдений в длинном виде - отнюдь не эквивалентные способы подачи данных на glm(), обрабатывает их по-разному. Я не шибко смыслю в GLM, но помню, что при подгонке распределений случайной величины по группированным данным иногда применяют поправки Шепарда и при построении моделей имеют лишь по одной точки для каждого класс-интервала. Возможно, и здесь используется похожий принцип? Но терзают меня сомнения в оправданности такого подхода, ведь в данном случае это естественная группировка, а не искусственная. В общем, прошу просветить меня по этому вопросу. Сообщение было отредактировано ИНО - 14.01.2021 16:43 |
plantago Постоянный участник |
|
ИНО Постоянный участник Донецк |
|
guest: taras IP-штамп: frBr/JzPmtwyY гость |
помогите разобраться с вопросом. Имеются многолетние данные по продуктивности вида в разных условиях местообитания. Я сделал базовую статистики, сравнение по годам и анализ временных изменений в каждом из них по отдельности. Однако не пойму как сделать сравнение между несколькими временными рядами. Т.е. сказать, что в таких-то условиях виду хорошо на протяжении длительного времени несмотря на изменение условий, а в других ему "хуже" |
ИНО Постоянный участник Донецк |
|
guest: taras IP-штамп: frBr/JzPmtwyY гость |
Графики 1- изменение продуктивности от глубины воды 2- сводный ts.plot |
guest: taras IP-штамп: frBr/JzPmtwyY гость |
Прямые ссылки |
ИНО Постоянный участник Донецк |
(guest: taras @ 18.01.2021 00:57) Потому что вы их неправильно вставляли. Надо так: |
ИНО Постоянный участник Донецк |
Насчет графика - видна ярко выраженная сезонность, выраженного многолетнего трэнда не видно, но с этим, наверное Вы уже разобрались при анализе каждого ряда отдельно. А что именно хотите от сравнения трех? Видно что динамика в целом похожая. Видно, что "ph" было значительно лучше других в 17 г., а С.el аномально поплохело между второй и третьей точками 19 г., но очень сомневаюсь, что к этому можно прикрутить статистическую значимость. Вот если б Вы в каждый момент показатель продуктивности оценивали не точечно, а интервально, то можно было бы посмотреть на перекрытие ДИ и не только... Как измеряли продуктивность эту? Если на графике - все имеющиеся данные, то лично я бы ограничился словесным описанием. Можно конечно спрогнозировать, что и дальше в теплое время года показатель будет расти, а в холодное падать, но не уверен что в этом есть какая-либо научная новизна. А многолетнюю динамику тут особо не спргнозируешь. И главное: на графике метки времени совсем не те, которые Вы описали, а январь, май и по одной между ними (всего по 4 за год). Как так? |
Guest IP-штамп: frBr/JzPmtwyY гость |
первый какое имеет отношение к описанной вами проблеме? Пытался отобразить изменение продуктивности в зависимости от уровня воду за период наблюдения. И главное: на графике метки времени совсем не те, которые Вы описали, а январь, май и по одной между ними (всего по 4 за год). Как так? Это как-то странно подписи к графику построились. Данные первого квартала -0. А что именно хотите от сравнения трех? Мы в целом знаем, что наблюдается снижение параметров продуктивности вида в зависимости от снижения уровня воды. Но хотелось найти интервалы для различных местообитаний. Во вторых посмотреть взаимосвязь между продуктивностью постоянно встречающихся групп видов. Например тростник, осоки, болотное разнотравье. Тренды для них выявлены, но как их сравнивать? Вот если б Вы в каждый момент показатель продуктивности оценивали не точечно, а интервально, то можно было бы посмотреть на перекрытие ДИ и не только... Как измеряли продуктивность эту? Брался укос растений с 0,25м2 с 2-3 кратной повторностью. Вот если б Вы в каждый момент показатель продуктивности оценивали не точечно, а интервально, то можно было бы посмотреть на перекрытие ДИ и не только... Можно по подробнее? В целом хочется научится вытаскивать как можно больше данных за командировки (их мало и как правило за свой счет). Возможно надо менять схему сбора данных? |
ИНО Постоянный участник Донецк |
Что-то никаких трендов на ваших графиках я не увидел. Или есть еще данные? Про интервальную оценку. Если у Вас для каждой точки времени есть выборка значений продуктивности, то можно построить ДИ для метематического ожидания этого параметра в генеральной совокупности. Увы, двух-трехкратная повторность тут вряд ли сильно поможет, интервалы выйдут очень широкими. Но все равно держите в уме, что для каждого момента у Вас на самом деле есть не одна, а 2-3 точки, авось пригодятся. Но таки схему поменять не помешало бы: делать для каждого вида больше покосов в каждый момент времени (лучше одинаковое число). И не редуцировать этот набор данных до одного среднего арифиметического. Сколько именно укосов, зависит от величины различий в этом показателе между видами, биотопами или что у вас там. Если она очень велика, то может даже хватит и имеющихся двух-трех. Попробуйте построить ДИ для среднего по имеющимся наблюдениям (тем самым 2-3 укосам) в одной временной точке для разных видов (биотопов или что у Вас там), если они не перекроются, то Вам крупно повезло. Но при таком примитивном подходе, по каждой точке отдельно, естественно, сильно страдает мощность. Более мощной альтернативной может быть построение нелинейной (можно непараметрической, ядерной, например) кривой регрессии от времени. Встречал я такой подход к коротким временным рядам. Более того, если нет явных погодичных изменений (а я их не вижу, кроме, может быть, "рн" в 2017), то можно сделать ход конем и замкнуть кривые в подобие окружностей, соответствующих одному году. Тогда на каждый квартал придется у же не по 2-3 укоса, а по 6-9, что уже кое-что. Для такого финта используются методы статистики угловых величин (англ. "circular statistics"). Но, наверное, подобное можно и через традиционные методы анализа временных рядов решить, я с ними очень плохо знаком. То немногое, из арсенала этих методов, с чем мне доводилось сталкиваться, заточено сугубо на прогноз и не представляет интереса в данном случае. Знаю, более опытные пользователи данного форума разбираются во временных рядах намного лучше, но они что-то не торопятся отвечать... И все еж настоятельно рекомендую разобраться с кодировкой данных, в норме значение временной переменной не должно браться программой из ниоткуда по собственной воле. Если забьете на эту маленькую техническую неувязку, то где гарантия, что все остальное она вам посчитает верно? |
guest: taras IP-штамп: frBr/JzPmtwyY гость |
|
ИНО Постоянный участник Донецк |
Во анализе же временных рядов обычно после выделения сезонного компонента, в причины его обуславливающие глубоко не копают, потому как для прогноза это не требуется, а цель там практически всегда в прогнозе. Не думаю, что это то, что нужно Вам. Еще один важный экологический момент - значение 0 в приросте биомассы имеет глубокий смысл. Во-первых на модель полагаются дополнительные ограничение, чтоб регрессионная кривая (или прямая) не ушла даже капельку ниже него, потому как это будет с точки зрения экологии нонсенс. Во-вторых, этот ноль на самом деле не точечное значение, он держится на протяжении довольно значительно периода времени. То есть, если Вы будете строить регрессию от времени, то правильная с экологической точки зрения кривая должна быть похожей на синусоиду с горбами, подрезанными снизу (в этих местах она будет превращаться в отрезки прямой со значением 0). Если же будете строить линейную регрессию от температуры воздуха то лучше взять не простую, а эффективную температуру и исключить свободный член (прямая должна проходить через начало координат). Надеюсь, теперь я действительно дал достаточно пищи для размышлений. |
Guest IP-штамп: frBr/JzPmtwyY гость |
Пообщался с народом и они то же предложили использовать линейные регрессии. CODE # p/p_Car_acut — проективное покрытие осоки # p/p_Ph - проективное покрытие тростника # p/p_other_sp - проективное покрытие других встречающихся видов # kol-vo_Ph — количество стеблей тростника на 1 м2 # mass_Ph — вес стеблей тростника на 1 м2 # mass_other_sp — вес остальной растительности # total_mass — общий вес # wat_depth — глубина воды (см). Является главным фактором существования рассматриваемого участка. p<-read.table("plot_1.dat",header=TRUE,dec=",") # year Q p.p_Car_acut p.p_Ph p.p_other_sp kol.vo_Ph mass_Ph mass_other_sp total_mass wat_depth p$fQ=factor(p$Q,levels=c("I","II","III","IV")) p$fYear=factor(p$year) i.sel=3 if(i.sel == 1) M1=lm(p.p_Ph~fQ/wat_depth,data=p) if(i.sel == 2) M1=lm(kol.vo_Ph~fQ/wat_depth,data=p) if(i.sel == 3) M1=lm(mass_Ph~fQ/wat_depth,data=p) summary(M1) old.par=par(mfrow=c(2,2)) plot(M1,which=1) plot(M1,which=2) plot(M1,which=3) plot(M1,which=4) par(old.par) i.sel=1 if(i.sel == 1) M2=lm(p.p_Ph~fQ*wat_depth,data=p) if(i.sel == 2) M2=lm(kol.vo_Ph~fQ*wat_depth,data=p) if(i.sel == 3) M2=lm(mass_Ph~fQ*wat_depth,data=p) summary(M2) Теперь разбираюсь. В частности не понятна команда i.sel= и почему их используется два типа? И почему вывод объединен? |
ИНО Постоянный участник Донецк |
Похоже на кусок кода из недр некой пользовательской функции, самостоятельно смысла не имеет. А кварталы куда делись? Или это уже другая задача не связанная временными рядами на графике №2. Вообще столько разных показаателей, характерезующих одно и то же являение (продуктивность), при наличие всего лишь одного предиктора не есть хорошо. Есть опасность обнаружить статистически значимое влияния этого вашего уровня воды на северное сияние Выберете лучше что-то одно, что легче будет интерпретировать в контексте вашего исследования. Подозреваю, что это должна быть биомасса. И да, настаиваю, что температура важнее уровня воды, потому что при 0 градусов, какой бы там уровень ни был, а расти нифига не будет Так что учитывать придется или ее непосредственно, или через время года. Архивы метеостанций есть в открытом доступе, проблем со сбором данных не будет, если только участки эти расположены не совсем в заднице мира. То же самое, наверное, с фотопериодом, там с получением данных все просто даже для самой задницы мира Кроме того, если тут участвуют данные с рис. 1, то lm() с одним предиктором им не поможет. Думаю, лучше обойтись без этих коллег, если только они не возьмут на себя анализ и интерпретацию результатов в полном объеме. |
Guest IP-штамп: frBr/JzPmtwyY гость |
Кварталы врде стали фактором CODE p$fQ=factor(p$Q,levels=c("I","II","III","IV"))[quote].
Наверно последую Вашему совету и сделаю регрессию по отдельности. [quote]И да, настаиваю, что температура важнее уровня воды, потому что при 0 градусов, какой бы там уровень ни был, а расти нифига не будет umnik.gif Так что учитывать придется или ее непосредственно, или через время года. Архивы метеостанций есть в открытом доступе, проблем со сбором данных не будет, если только участки эти расположены не совсем в заднице мира.[/quote] С этим очень сложно. Суть в том, что температура на метеостанции значительно отличается от температур в различных ценозах. влажность воздуха так же очень сильно отличается (например болота и степь). Я пробую сейчас два варианта. Первый - измерение локальных микроклиматических показателей с помощью микрометеостанции. Даже патент получил :shuffle: Выявлено, что вертикальный градиент температуры, влажности в различных ценозах стат. значимо отличается и связан с его состоянием. Второй - Использование космоснимков и различных индексов. [quote]Думаю, лучше обойтись без этих коллег, если только они не возьмут на себя анализ и интерпретацию результатов в полном объеме.[/quote] Сам буду делать. Попробую так по каждому параметру [code]p<-read.table("plot_1.dat",header=TRUE,dec=",") p$fQ=factor(p$Q,levels=c("I","II","III","IV")) p$fYear=factor(p$year) M1=lm(p.p_Ph~fQ/wat_depth,data=p) summary(M1) old.par=par(mfrow=c(2,2)) plot(M1,which=1) plot(M1,which=2) plot(M1,which=3) plot(M1,which=4) par(old.par) |
ИНО Постоянный участник Донецк |
Данные метеостанции, конечно отличаются, от имеющихся в биотопе, но весьма тесно кореллируют с ними. Этого хватит. Микрометеостанция - это интересно, расскажите подробнее. Мне такой штуки очень не хватает. У нас мало того что с 14 года метеостанции нету вообще, так еще и частота обновления для мох задач нужна гораздо более высокая, чем обычно используют метеорологи. Но то на будущее, а сейчас Вам же ж надо как-то обработать уже собранные данные, не выкидывать же их. Тут вполне и данные городской метеостанции сойдут. Без данных о температуре (хотя бы приблизительных) линейную модель там просто не построить. |
ИНО Постоянный участник Донецк |
|
Guest IP-штамп: frBr/JzPmtwyY гость |
Микрометеостанция - это интересно, расскажите подробнее патент UA130334U. Но этот девайс уронили с квадрокоптера. Сейчас пробуем создать беспроводную версию с передачей данный по Wi-Fi. Есть еще 2 статьи, где использовал данные, полученные с помощью прибора. Первая THE TECHNOLOGY NEW INSTRUMENTAL ESTIMATION METHOD OF THE VEGETATION COVER MICROCLIMATIC CHARACTERISTICS (есть в открытом доступе), вторая ОЦІНКА ЗМІН МІКРОКЛІМАТИЧНИХ УМОВ В УГРУПОВАННЯХ PHRAGMITETUM AUSTRALIS В ДЕЛЬТІ ДНІСТРА ПІД ДІЄЮ ІНВАЗІЇ ZIZANIA LATIFOLIA (пока закрыта) |
ИНО Постоянный участник Донецк |
|
Guest IP-штамп: frBr/JzPmtwyY гость |
Особенностью является одновременный опрос датчиков. Сейчас делаем клиент-серверную структуру. клиентами выступают модуль управления датчиками (записывает и сохраняет данные) и по команде сервера передает накопленные данные. Сервер позволяет опрашивать и задавать алгоритм работы клиентов. Скорость ветра очень важный параметр, но мне только в прыжке его замерить Высота тростника 3-4 метра. Хотя если можно будет ставить датчики стационарно (с надеждой, что их не сукраинят), то можно и анемометр поставить. Все это делаем за свой счет и работа идет туго.
|
ПолинаШ Участник |
|
plantago Постоянный участник |
|
ПолинаШ Участник |
Но, видимо, придется доваивать дендрограммы в фотошопе... |
plantago Постоянный участник |
Есть еще ggcorrplot, можно посмотреть, что он умеет. |
d-taras |
Считаю индекс Жаккара для ряда описаний используя следующие команды CODE Jaccard = function (x, y) { M.11 = sum(x == 1 & y == 1) M.10 = sum(x == 1 & y == 0) M.01 = sum(x == 0 & y == 1) return (M.11 / (M.11 + M.10 + M.01)) } input.variables = data.frame(tab) m = matrix(data = NA, nrow = length(input.variables), ncol = length(input.variables)) for (r in 1:length(input.variables)) { for (c in 1:length(input.variables)) { if (c == r) { m[r,c] = 1 } else if (c > r) { m[r,c] = Jaccard(input.variables[,r], input.variables[,c]) } } } variable.names = sapply(input.variables, attr, "label") colnames(m) = variable.names rownames(m) = variable.names jaccards = m data.frame получаю командой CODE tab<-as.data.frame(mytable) Результаты правильные, но почему-то названия теряются CODE NULL NULL NULL NULL NULL 1 0.7032967 0.2407407 0.4715447 NULL NA 1.0000000 0.2551020 0.5740741 NULL NA NA 1.0000000 0.1825397 NULL NA NA NA 1.0000000 Хотя названия столбцов присутствует и должно быть Carpinus, Oak, e.t.c. Как прямо задать название столбцов и ограничить вывод 2 знаками после запятой? |
plantago Постоянный участник |
> mynames <- c("Carpinus", "Oak", "Orpinus", "Boak") > set.seed(0) > mytable <- data.frame(sample(0:1, 10, replace=TRUE), + sample(0:1, 10, replace=TRUE), + sample(0:1, 10, replace=TRUE), + sample(0:1, 10, replace=TRUE)) > names(mytable) <- mynames > str(mytable) 'data.frame': 10 obs. of 4 variables: $ Carpinus: int 0 1 0 0 1 1 0 1 1 1 $ Oak : int 1 1 1 1 0 1 1 0 0 0 $ Orpinus : int 0 0 1 0 0 0 1 1 1 0 $ Boak : int 1 0 0 1 1 1 0 1 1 1 > tab <- as.data.frame(mytable) > Jaccard <- function (x, y) { + M.11=sum(x == 1 & y == 1) + M.10=sum(x == 1 & y == 0) + M.01=sum(x == 0 & y == 1) + return (M.11 / (M.11 + M.10 + M.01)) + } > input.variables <- data.frame(tab) > m <- matrix(data=NA, nrow=length(input.variables), ncol=length(input.variables)) > for (r in 1:length(input.variables)) { + for (c in 1:length(input.variables)) { + if (c == r) { + m[r,c] = 1 + } else if (c > r) { + m[r,c] = Jaccard(input.variables[,r], input.variables[,c]) + } + } + } > variable.names <- names(input.variables) > colnames(m) <- variable.names > rownames(m) <- variable.names > jaccards <- m > jaccards Carpinus Oak Orpinus Boak Carpinus 1 0.2857143 0.0 0.4285714 Oak NA 1.0000000 0.2 0.3750000 Orpinus NA NA 1.0 0.1666667 Boak NA NA NA 1.0000000 Но только какой-то странный у Вас Jaccard (я не проверял, и не оптимизировал ничего). Потому что: > library(vegan) Loading required package: permute Loading required package: lattice This is vegan 2.5-7 > vegdist(t(mytable), method="jaccard") Carpinus Oak Orpinus Oak 0.7142857 Orpinus 1.0000000 0.8000000 Boak 0.5714286 0.6250000 0.8333333 |
d-taras |
(plantago @ 24.02.2021 15:35) А где же данные? Пришлось мне их за Вас синтезировать. Но только какой-то странный у Вас Jaccard (я не проверял, и не оптимизировал ничего). Потому что: > library(vegan) > vegdist(t(mytable), method="jaccard") Carpinus Oak Orpinus Oak 0.8000000 Orpinus 0.7500000 0.7500000 Boak 0.3750000 0.7000000 0.7777778 Спасибо за ответ! Данные, если ещё интересно, приложил. Я думал, что что-то в тексте не так, что название не выводятся. Что же касается индекса Жаккара, имеется множество вариаций и пользуются тем к чему привыкли. Разобрался.... команда получается такая CODE 1 - vegdist(t(mydata), method = "jaccard") Почему-то файл не цепляется... Сообщение было отредактировано d-taras - 24.02.2021 18:59
|
ИНО Постоянный участник Донецк |
|
plantago Постоянный участник |
|
Алекс3212 |
A B C 10 3 13 в R C <- ((A*100)+B) в excel C= A&B |
Алекс3212 |
paste () , где в кавычках перечисляете столбцы которые надо присоединить |
plantago Постоянный участник |
|
ИНО Постоянный участник Донецк |
|
ИНО Постоянный участник Донецк |
CODE library(spatstat) load("сн_2017Pgppp.Rda") load("сн_2017Pg_lL.Rda") plot(сн_2017Pgppp) Задача - залить их цветом в зависимости от значений сн_2017Pg_lL, и приделать сбоку цветовую шкалу в виде б. м. непрерывного градиента и подписями значений соответствующего показателя через равные интервалы. Файл/ы:
|
plantago Постоянный участник |
(ИНО @ 06.09.2022 11:37) Прицепом еще один вопрос. Есть у меня точки на карте (файлы прилагаются): CODE library(spatstat) load("сн_2017Pgppp.Rda") load("сн_2017Pg_lL.Rda") plot(сн_2017Pgppp) Задача - залить их цветом в зависимости от значений сн_2017Pg_lL, и приделать сбоку цветовую шкалу в виде б. м. непрерывного градиента и подписями значений соответствующего показателя через равные интервалы. CODE library(spatstat) load("сн_2017Pgppp.Rda") load("сн_2017Pg_lL.Rda") # plot(сн_2017Pgppp) # error! aa <- (get(ls()[1])) # does not always work; if not, use rm(list=ls()) _before_ load() bb <- (get(ls()[2])) # and please do not name R objects with anything except [A-z0-9_.] rr <- rank(aa, ties.method="last") cc <- terrain.colors(max(rr)) plot(bb, cols=0) points(bb$x, bb$y, pch=1, col=cc[rr]) # next step: deal with legend Про шкалу можно, наверное, посмотреть сюда |
plantago Постоянный участник |
(ИНО @ 06.09.2022 10:49) Как заставить ctree() из {party} работать не с положением (средними), а с масштабом (дисперсией) зависимой переменной? Насколько понял из руководства, для этой цели можно использовать трансформацию при помощи аргумента ytrafo, но как именно им пользоваться так и не уразумел. Что только не перепробовал и с готовыми ytrafo из {coin}, на которые дается ссылка в том же руководстве и с самопиской функцией, считающей абсолютное отклонение от среднего, - либо матерится разными матюками, либо молча считает различия в средних, точно так же как с настройками по умолчанию, будто бы никаких манипуляций с ytrafo и не было. К сожалению, без конкретных примеров и ссылок не смог ничего понять. |
ИНО Постоянный участник Донецк |
CODE > plot(bb, cols=0) Warning messages: 1: In plot.window(...) : "cols" is not a graphical parameter 2: In plot.xy(xy, type, ...) : "cols" is not a graphical parameter 3: In axis(side = side, at = at, labels = labels, ...) : "cols" is not a graphical parameter 4: In axis(side = side, at = at, labels = labels, ...) : "cols" is not a graphical parameter 5: In box(...) : "cols" is not a graphical parameter 6: In title(...) : "cols" is not a graphical parameter > points(bb$x, bb$y, pch=1, col=cc[rr]) # next step: deal with legend Error in bb$x : $ operator is invalid for atomic vectors При том выдается не карта, а обычный скаттерплот с незакрашенными точками. Возможно играют роль версии ПО, в моем случае это R 3.3.1 и spatstat 1.46-1. В то время как plot(сн_2017Pgppp) работает вполне корректно у меня, никаких error. И проблемы с кириллицей видел только в заголовках столбцов таблиц. Наверное, дело в операционке. Ну, уберите из названия объекта эти две незнакомые буржуям буквы, если мешают, это ж просто пример, мне б принцип понять. С деревом можно никаких файлов не грузить, а обойтись простым синтетическим примером: CODE y<-c(rnorm(10, 0, 1), rnorm(10, 0, 100)) x<-factor(c(rep(1,10), rep(2,10))) xy<-cbind(x, y) ctree(y~x, data=xy)#ветвление дерева нет, так как проверяется разница в средних, а надо чтобы проверялась разница в дисперсиях (или любом другом параметре масштаба). xy$z<-c(abs(y[1:10]-mean(y[1:10])), abs(y[1:20]-mean(y[1:20])))#модуль отклонений от среднего в каждой группе ctree(z~x)#теперь различия уловлены правильно Надо чтоб преобразование y в z проходило автоматически в зависимости от формулы, переданной ctree(). По идее для этого есть в ctree() аргумент ytrafo, но как его правильно задавать, я так и не понял. |
23432432432432dd IP-штамп: fryfA2WSUbN5E гость |
|
« Предыдущая тема · Биофизика и матметоды в биологии · Следующая тема » |