Rambler's Top100
Лёгкая версия форума* Виртуальная клавиатура  English  
Molbiol.ru | О проекте | Справочник | Методы | Растворы | Расчёты | Литература | Орг.вопросы
Web | Фирмы | Coffee break | Картинки | Работы и услуги | Биржа труда | Междисциплинарный биологический онлайн-журналZbio-wiki

NG SEQUENCING · ЖИЗНЬ РАСТЕНИЙ · БИОХИМИЯ · ГОРОДСКИЕ КОМАРЫ · А.А.ЛЮБИЩЕВ · ЗООМУЗЕЙ


Темы за 24 часа  [ Вход* | Регистрация* ]  
   



Форум: 
 

Щёлкните, чтобы внести в Избранные Темы* R Help -- Давайте составим русский FAQ --
Кураторы темы:* plantago
Операции: Хочу стать куратором* · Подписаться на тему* · Отправить страницу по e-mail · Версия для печати*
Внешний вид:* Схема · [ Стандартный ] · +Перв.сообщ.


Добавить сообщение в темуСоздать новую темуСоздать голосование
Участник оффлайн! ИНО
Постоянный участник
Донецк



 прочитанное сообщение 09.01.2021 21:23     Сообщение для модератора  Сообщение для куратора темы       Личное письмо

Plantago, спасибо, очень хитро! А не могли бы Вы завернуть этот код в одну функцию, которую можно было бы переменить для таблицы с любым числом столбцов, передав ей имена двух столбцов, содержащих частоты? А то можно легко ошибиться, меняя номера столбцов вручную в нескольких местах каждый раз при переходе на таблицу с другим числом столбцов.
Участник оффлайн! ИНО
Постоянный участник
Донецк



 прочитанное сообщение 10.01.2021 05:27     Сообщение для модератора  Сообщение для куратора темы       Личное письмо

Прочитал я на с.123 книги Мастицкий С.Э., Шитиков В.К."Статистический анализ и визуализация данных с помощью R" пример скармливания glm() для построения модели логистической регрессии непосредственно таблицы частот в виде отдельной матрицы (что сделало бы ненужной предыдущую возню с реструктуризацией данных). Опробовал - считает что-то и не ругается (хотя при этом получается, что левая часть формулы - один объект объекта, а правая - часть другого объекта, что для меня несколько странно). Однако же считает нечто совсем не то, что в случае ее применения к развернутой в длинный вид таблице. Для примера возьмем созданные в предыдущих постах dataframe1 и df3. В последней только надо тип данных поменять сначала:
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
Постоянный участник



 прочитанное сообщение 10.01.2021 14:45     Сообщение для модератора  Сообщение для куратора темы       Личное письмо  Отправить e-mail  Web-адрес

> ... таблицы частот в виде отдельной матрицы ...

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

> А не могли бы Вы завернуть этот код в одну функцию

Не думаю, что это разумно, потому что данные могут быть очень разными. Проще понять способ, записать себе его с комментариями и применять, модифицируя.

Всего благодарностей: 1Поблагодарили (1): ИНО
Участник оффлайн! ИНО
Постоянный участник
Донецк



 прочитанное сообщение 10.01.2021 15:47     Сообщение для модератора  Сообщение для куратора темы       Личное письмо

Ну так вроде я именно такую матрицу и подал, в левом столбце "numbers of successes", в правом "numbers of failures". Почему же при двух, по идее, эквивалентных способах подачи данных модели получились разные? В чем ошибка?

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

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


confused.gif

Сообщение было отредактировано ИНО - 10.01.2021 16:59
Участник оффлайн! ПолинаШ
Участник



 прочитанное сообщение 13.01.2021 12:27     Сообщение для модератора  Сообщение для куратора темы       Личное письмо  Отправить e-mail

Как сделать выборку строк из таблицы с полем 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>

И простите за ламерский вопрос: у меня перестали работать все иконки панели редактора сообщений кроме кнопок "Пред.просмотр" и "Ответить"
Guest
IP-штамп: frD0EXpX.eS.Q
гость



 прочитанное сообщение 13.01.2021 13:13     Сообщение для модератора  Сообщение для куратора темы     

(ПолинаШ @ 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


Всего благодарностей: 1Поблагодарили (1): ПолинаШ
Участник оффлайн! ИНО
Постоянный участник
Донецк



 прочитанное сообщение 14.01.2021 16:39     Сообщение для модератора  Сообщение для куратора темы       Личное письмо

Никто мне не ответил, так что я исследовал по мере своих весьма скромных возможностей владения R описанную проблему самостоятельно. Повторил пример из Щитикова и Мастицкого - загрузил файл dreissena_infection.txt, там данные в виде таблицы частот, как в моей dataframe1. Задал модель логистической регрессии функцией glm() - результаты в точности сошлись с книгой. Там же есть файл с данными в длинном виде, как в моей dataframe1, dreissena_infection_raw_data.txt, но, как оказалось, в нем одного наблюдения не хватает! Хотя в книге постулируется, будто это те же данные и по ним построена вторая модель, мало похожая на первую, что авторов почему-то совсем не смутило eek.gif

Я же развернул 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
Постоянный участник



 прочитанное сообщение 14.01.2021 18:56     Сообщение для модератора  Сообщение для куратора темы       Личное письмо  Отправить e-mail  Web-адрес

Правильнее всего, мне кажется, связаться по этому поводу с самими Шитиковым и/или Мастицким. По-моему, они вполне доступны для вопросов, например, через их блоги.
https://r-analytics.blogspot.com/p/blog-page.html
https://stok1946.blogspot.com/
Участник оффлайн! ИНО
Постоянный участник
Донецк



 прочитанное сообщение 14.01.2021 19:32     Сообщение для модератора  Сообщение для куратора темы       Личное письмо

Так у меня вопрос вовсе не о странностях в их книгах (там их на самом деле сильно больше одной), а о нюансах задании аргументов функции glm(). Неужто никто на форуме не пользуется логистической или пробит регрессией? Это ж, вроде бы, один из самых распространенных видов регрессионного анализа.
guest: taras
IP-штамп: frBr/JzPmtwyY
гость



 прочитанное сообщение 17.01.2021 20:15     Сообщение для модератора  Сообщение для куратора темы     

Приветствую всех!
помогите разобраться с вопросом.
Имеются многолетние данные по продуктивности вида в разных условиях местообитания. Я сделал базовую статистики, сравнение по годам и анализ временных изменений в каждом из них по отдельности.
Однако не пойму как сделать сравнение между несколькими временными рядами. Т.е. сказать, что в таких-то условиях виду хорошо на протяжении длительного времени несмотря на изменение условий, а в других ему "хуже"
Участник оффлайн! ИНО
Постоянный участник
Донецк



 прочитанное сообщение 17.01.2021 21:47     Сообщение для модератора  Сообщение для куратора темы       Личное письмо

Для начала посмотреть глазами на графики, далее - в зависимости от увиденного. Или Вам надо просто какое-нибудь p для галочки?
guest: taras
IP-штамп: frBr/JzPmtwyY
гость



 прочитанное сообщение 17.01.2021 23:54     Сообщение для модератора  Сообщение для куратора темы     

Данные собирались в течение 3-х лет 3 раза за вегетацию (2-4 квартал т.е. Июнь, Сентябрь, Ноябрь).
Графики
1- изменение продуктивности от глубины воды
2- сводный ts.plot
user posted image
user posted image
guest: taras
IP-штамп: frBr/JzPmtwyY
гость



 прочитанное сообщение 17.01.2021 23:57     Сообщение для модератора  Сообщение для куратора темы     

Почему-то рисунки не вставились
Прямые ссылки
https://ibb.co/zHGVhcj][img]https://i.ibb.co/vvZ4x5g/depth-massa.gif
https://ibb.co/k0sdYX1][img]https://i.ibb.co/pZBmCdy/massa-ts.gif
Участник оффлайн! ИНО
Постоянный участник
Донецк



 прочитанное сообщение 18.01.2021 00:05     Сообщение для модератора  Сообщение для куратора темы       Личное письмо

(guest: taras @ 18.01.2021 00:57)
Ссылка на исходное сообщение  Почему-то рисунки не вставились


Потому что вы их неправильно вставляли. Надо так:

user posted image
user posted image
Участник оффлайн! ИНО
Постоянный участник
Донецк



 прочитанное сообщение 18.01.2021 00:47     Сообщение для модератора  Сообщение для куратора темы       Личное письмо

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

Насчет графика - видна ярко выраженная сезонность, выраженного многолетнего трэнда не видно, но с этим, наверное Вы уже разобрались при анализе каждого ряда отдельно. А что именно хотите от сравнения трех? Видно что динамика в целом похожая. Видно, что "ph" было значительно лучше других в 17 г., а С.el аномально поплохело между второй и третьей точками 19 г., но очень сомневаюсь, что к этому можно прикрутить статистическую значимость. Вот если б Вы в каждый момент показатель продуктивности оценивали не точечно, а интервально, то можно было бы посмотреть на перекрытие ДИ и не только... Как измеряли продуктивность эту?

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

И главное: на графике метки времени совсем не те, которые Вы описали, а январь, май и по одной между ними (всего по 4 за год). Как так?
Guest
IP-штамп: frBr/JzPmtwyY
гость



 прочитанное сообщение 18.01.2021 01:41     Сообщение для модератора  Сообщение для куратора темы     

Спасибо за ответ.
первый какое имеет отношение к описанной вами проблеме?

Пытался отобразить изменение продуктивности в зависимости от уровня воду за период наблюдения.
И главное: на графике метки времени совсем не те, которые Вы описали, а январь, май и по одной между ними (всего по 4 за год). Как так?

Это как-то странно подписи к графику построились. Данные первого квартала -0.
А что именно хотите от сравнения трех?

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

Брался укос растений с 0,25м2 с 2-3 кратной повторностью.
Вот если б Вы в каждый момент показатель продуктивности оценивали не точечно, а интервально, то можно было бы посмотреть на перекрытие ДИ и не только...

Можно по подробнее? В целом хочется научится вытаскивать как можно больше данных за командировки (их мало и как правило за свой счет). Возможно надо менять схему сбора данных?
Участник оффлайн! ИНО
Постоянный участник
Донецк



 прочитанное сообщение 18.01.2021 21:23     Сообщение для модератора  Сообщение для куратора темы       Личное письмо

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

Что-то никаких трендов на ваших графиках я не увидел. Или есть еще данные?

Про интервальную оценку. Если у Вас для каждой точки времени есть выборка значений продуктивности, то можно построить ДИ для метематического ожидания этого параметра в генеральной совокупности. Увы, двух-трехкратная повторность тут вряд ли сильно поможет, интервалы выйдут очень широкими. Но все равно держите в уме, что для каждого момента у Вас на самом деле есть не одна, а 2-3 точки, авось пригодятся. Но таки схему поменять не помешало бы: делать для каждого вида больше покосов в каждый момент времени (лучше одинаковое число). И не редуцировать этот набор данных до одного среднего арифиметического. Сколько именно укосов, зависит от величины различий в этом показателе между видами, биотопами или что у вас там. Если она очень велика, то может даже хватит и имеющихся двух-трех. Попробуйте построить ДИ для среднего по имеющимся наблюдениям (тем самым 2-3 укосам) в одной временной точке для разных видов (биотопов или что у Вас там), если они не перекроются, то Вам крупно повезло. Но при таком примитивном подходе, по каждой точке отдельно, естественно, сильно страдает мощность. Более мощной альтернативной может быть построение нелинейной (можно непараметрической, ядерной, например) кривой регрессии от времени. Встречал я такой подход к коротким временным рядам. Более того, если нет явных погодичных изменений (а я их не вижу, кроме, может быть, "рн" в 2017), то можно сделать ход конем и замкнуть кривые в подобие окружностей, соответствующих одному году. Тогда на каждый квартал придется у же не по 2-3 укоса, а по 6-9, что уже кое-что. Для такого финта используются методы статистики угловых величин (англ. "circular statistics"). Но, наверное, подобное можно и через традиционные методы анализа временных рядов решить, я с ними очень плохо знаком. То немногое, из арсенала этих методов, с чем мне доводилось сталкиваться, заточено сугубо на прогноз и не представляет интереса в данном случае. Знаю, более опытные пользователи данного форума разбираются во временных рядах намного лучше, но они что-то не торопятся отвечать...

И все еж настоятельно рекомендую разобраться с кодировкой данных, в норме значение временной переменной не должно браться программой из ниоткуда по собственной воле. Если забьете на эту маленькую техническую неувязку, то где гарантия, что все остальное она вам посчитает верно?
guest: taras
IP-штамп: frBr/JzPmtwyY
гость



 прочитанное сообщение 18.01.2021 22:13     Сообщение для модератора  Сообщение для куратора темы     

Спасибо за советы, буду пробовать. С кол-вом укосов согласен, что мало, но стараюсь делать максимум за отведенное время... Как у Владимира Семеновича пелось "приехал на две недели, работы на восемь лет..."
Участник оффлайн! ИНО
Постоянный участник
Донецк



 прочитанное сообщение 19.01.2021 17:21     Сообщение для модератора  Сообщение для куратора темы       Личное письмо

Обмозговал я вашу проблему на досуге с экологической точки зрения и пришел к эта цикличноть с нулевым ростом зимой - она ж не спроста. Почему бы не включить в модель в качестве предикторов температуру воздуха и длину светового дня? Ну и уровень воды до кучи. А время убрать. Тогда авось все удастся свести к многофакторной линейной регрессии. А это дает большие плюшки в плане интерпритации результатов. Единственное, видится мне из графика, скорость роста в первой половине года и во второй неодинакова, причем неодинковость эта у разных кривых разная. Связано ли это с аналогичной асимметрией изменения уровня воды? Если нет, то, возможно с фазами вегетации растений? Например, до цветения, или до созревания плодов продуктивность плавно нарастает, а после - резко падает. Или наоборот. Понять в чем дело модно будет по регрессионным остаткам. Если там обнаружатся намекающие на такую ситуацию аномалии, то следует ввести еще один предиктор - фаза жизненного цикла. Кстати, возможно они сменяются асинхронно у разных видов, отсюда и разный характер асимметрии для них. В любом случае предикторов для линейной модели будет много, а наблюдений - не очень, посему придется либо отбирать лучшее подмноджество (что позволит одновременно сделать вывод о их экологической важности), либо редуцировать размерность, например при помощи PCA. Последнее хоть и сложнее, но тоже вполне интерпретируемо: например зависимость от первой главной компоненты, имеющую максимальные нагрузки температуры воздуха и фотопериода (почти уверен, что это будет именно так), и второй компоненты, которую будет нагружать, допустим, уровень воды, находит вполне логичное экологическое объяснение объяснение.

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

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

Надеюсь, теперь я действительно дал достаточно пищи для размышлений.
Guest
IP-штамп: frBr/JzPmtwyY
гость



 прочитанное сообщение 19.01.2021 18:15     Сообщение для модератора  Сообщение для куратора темы     

Оу, что за пару дней массив данных не соберешь это трудно объяснить начальству, когда командировку выбиваешь.... А то у них на море, гад, хочешь....
Пообщался с народом и они то же предложили использовать линейные регрессии.
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= и почему их используется два типа? И почему вывод объединен?
Участник оффлайн! ИНО
Постоянный участник
Донецк



 прочитанное сообщение 19.01.2021 20:04     Сообщение для модератора  Сообщение для куратора темы       Личное письмо

А это не команда, это числовой вектор единичной длины, то бишь скаляр smile.gif Вы ж его там сами назначили: сначала 3, потом 1, зачем х. з.

Похоже на кусок кода из недр некой пользовательской функции, самостоятельно смысла не имеет.

А кварталы куда делись? Или это уже другая задача не связанная временными рядами на графике №2.

Вообще столько разных показаателей, характерезующих одно и то же являение (продуктивность), при наличие всего лишь одного предиктора не есть хорошо. Есть опасность обнаружить статистически значимое влияния этого вашего уровня воды на северное сияние smile.gif Выберете лучше что-то одно, что легче будет интерпретировать в контексте вашего исследования. Подозреваю, что это должна быть биомасса. И да, настаиваю, что температура важнее уровня воды, потому что при 0 градусов, какой бы там уровень ни был, а расти нифига не будет umnik.gif Так что учитывать придется или ее непосредственно, или через время года. Архивы метеостанций есть в открытом доступе, проблем со сбором данных не будет, если только участки эти расположены не совсем в заднице мира. То же самое, наверное, с фотопериодом, там с получением данных все просто даже для самой задницы мира smile.gif

Кроме того, если тут участвуют данные с рис. 1, то lm() с одним предиктором им не поможет. Думаю, лучше обойтись без этих коллег, если только они не возьмут на себя анализ и интерпретацию результатов в полном объеме.
Guest
IP-штамп: frBr/JzPmtwyY
гость



 прочитанное сообщение 19.01.2021 20:27     Сообщение для модератора  Сообщение для куратора темы     

Спасибо, буду смотреть.
Кварталы врде стали фактором
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)
Участник оффлайн! ИНО
Постоянный участник
Донецк



 прочитанное сообщение 19.01.2021 20:49     Сообщение для модератора  Сообщение для куратора темы       Личное письмо

Что именно "по отдельности"? Если отдельно для каждого предиктора строить регрессии то из этого в вашем случае ничего хорошего не выйдет, лучше сначала построить со всеми возможными, благо немного их, а потом исключать поочередно. А если под отдельностями Вы имеете в виду разные показатели биомассы, то огорчу. Надо строить для одного, выбранного априори из теоретических соображений. Иначе придется корректировать ошибку первого рода, а с вашим-то количеством наблюдений это заведомо проигрышный путь. Вариант же "построить пр регрессии для каждого от каждого, вдруг что найдется, про то и напишем, а об остальном забудем" рассматривать не будем. Так "влияние на северное сияние" обычно и получают. Стройте модели с умом.

Данные метеостанции, конечно отличаются, от имеющихся в биотопе, но весьма тесно кореллируют с ними. Этого хватит.

Микрометеостанция - это интересно, расскажите подробнее. Мне такой штуки очень не хватает. У нас мало того что с 14 года метеостанции нету вообще, так еще и частота обновления для мох задач нужна гораздо более высокая, чем обычно используют метеорологи. Но то на будущее, а сейчас Вам же ж надо как-то обработать уже собранные данные, не выкидывать же их. Тут вполне и данные городской метеостанции сойдут. Без данных о температуре (хотя бы приблизительных) линейную модель там просто не построить.
Участник оффлайн! ИНО
Постоянный участник
Донецк



 прочитанное сообщение 19.01.2021 20:57     Сообщение для модератора  Сообщение для куратора темы       Личное письмо

существуют еще методы регрессии с многомерным откликом, где можно было бы корректно учесть все показатели продуктивности в одной модели, но там во-первых весьма сложная интерпретация результатов, а во-вторых нужно много больше наблюдений, чем в случае одномерного отклика. Мой совет как эколога: выбирайте показатели, основанные на биомассе, про проективное покрытие и количество стеблей забудьте. А категориальная переменная "осока или тростник" пойдет в качестве предиктора.
Guest
IP-штамп: frBr/JzPmtwyY
гость



 прочитанное сообщение 19.01.2021 22:09     Сообщение для модератора  Сообщение для куратора темы     

Микрометеостанция - это интересно, расскажите подробнее

патент UA130334U. Но этот девайс уронили с квадрокоптера. Сейчас пробуем создать беспроводную версию с передачей данный по Wi-Fi.
Есть еще 2 статьи, где использовал данные, полученные с помощью прибора. Первая THE TECHNOLOGY NEW INSTRUMENTAL ESTIMATION METHOD OF THE VEGETATION COVER MICROCLIMATIC CHARACTERISTICS (есть в открытом доступе), вторая ОЦІНКА ЗМІН МІКРОКЛІМАТИЧНИХ УМОВ В УГРУПОВАННЯХ PHRAGMITETUM AUSTRALIS В ДЕЛЬТІ ДНІСТРА ПІД ДІЄЮ ІНВАЗІЇ ZIZANIA LATIFOLIA (пока закрыта)
Участник оффлайн! ИНО
Постоянный участник
Донецк



 прочитанное сообщение 19.01.2021 22:48     Сообщение для модератора  Сообщение для куратора темы       Личное письмо

Ну, хотя бы датчики перечислите, которые туда входят. В патенте только фотка "стандартного датчика". Что за датчик такой, какие показатели измеряет? Вообще конкретно с температурой есть решение очень простое и компактное - термохрон. Бывает еще и с гигрометром - гигрохрон. А вот с анемометром не бывает frown.gif
Guest
IP-штамп: frBr/JzPmtwyY
гость



 прочитанное сообщение 19.01.2021 23:25     Сообщение для модератора  Сообщение для куратора темы     

К прибору подключается 3 одинаковых модуля содержащие датчики температуры и влажности (SHT31-D IIC), прямой и обратной освещенности (GY-49 MAX44009 I2C) 2 шт на верхней и нижней стороне платы. В самом приборе стоит датчик давления (BMP180). Так же есть дополнительные датчики СО2 и пыли (можно и другие). И подключаем если надо в городе проводить измерения.
Особенностью является одновременный опрос датчиков.
Сейчас делаем клиент-серверную структуру. клиентами выступают модуль управления датчиками (записывает и сохраняет данные) и по команде сервера передает накопленные данные. Сервер позволяет опрашивать и задавать алгоритм работы клиентов.
Скорость ветра очень важный параметр, но мне только в прыжке его замерить smile.gif Высота тростника 3-4 метра. Хотя если можно будет ставить датчики стационарно (с надеждой, что их не сукраинят), то можно и анемометр поставить.
Все это делаем за свой счет и работа идет туго.

Всего благодарностей: 1Поблагодарили (1): ИНО
Участник оффлайн! ПолинаШ
Участник



 прочитанное сообщение 17.02.2021 14:56     Сообщение для модератора  Сообщение для куратора темы       Личное письмо  Отправить e-mail

Можно ли как-то объединить функционал heatmap и corrplot в том смысле, что вместо заливки ячеек разным цветом, как это в тепловой карте, были бы кружочки разного размера (method="circle" в corrplot)
Участник оффлайн! plantago
Постоянный участник



 прочитанное сообщение 18.02.2021 05:43     Сообщение для модератора  Сообщение для куратора темы       Личное письмо  Отправить e-mail  Web-адрес

Не понимаю. Что нужно здесь https://cran.r-project.org/web/packages/cor...plot-intro.html и как именно поменять?
Участник оффлайн! ПолинаШ
Участник



 прочитанное сообщение 22.02.2021 12:34     Сообщение для модератора  Сообщение для куратора темы       Личное письмо  Отправить e-mail

Есть исходная неквадратная матрица T[n, m]. Функция heatmap заливает каждую ячейку цветом, пропорциональным значению числа в этой ячейке. И строит вверху и слева две дендрограммы. При публикации в журналах черно-серо-белая палитра сливается и выглядит неряшливо. Другое дело - серые яица разного размера как в corrplot(). Сами яица нарисовать - не проблема - это обычный scatterplot(). Но хочется добавить дендрограммы, а не боксплоты.
Но, видимо, придется доваивать дендрограммы в фотошопе...
Участник оффлайн! plantago
Постоянный участник



 прочитанное сообщение 22.02.2021 14:53     Сообщение для модератора  Сообщение для куратора темы       Личное письмо  Отправить e-mail  Web-адрес

Если нужно просто добавить дендрограммы, используйте layout(). Кстати, corrplot() должен уметь произвольные матрицы с "is.corr=FALSE, order='original'"
Есть еще ggcorrplot, можно посмотреть, что он умеет.
Участник оффлайн! d-taras




 прочитанное сообщение 24.02.2021 10:04     Сообщение для модератора  Сообщение для куратора темы       Личное письмо  Отправить e-mail

Приветствую всех!
Считаю индекс Жаккара для ряда описаний используя следующие команды
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
Постоянный участник



 прочитанное сообщение 24.02.2021 14:35     Сообщение для модератора  Сообщение для куратора темы       Личное письмо  Отправить e-mail  Web-адрес

А где же данные? Пришлось мне их за Вас синтезировать.

> 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




 прочитанное сообщение 24.02.2021 18:56     Сообщение для модератора  Сообщение для куратора темы       Личное письмо  Отправить e-mail

(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

Всего благодарностей: 1Поблагодарили (1): plantago
Участник оффлайн! ИНО
Постоянный участник
Донецк



 прочитанное сообщение 07.03.2021 21:58     Сообщение для модератора  Сообщение для куратора темы       Личное письмо

Есть ли пакеты с функцией определения продолжительности светового дня по дате и географическим координатам?
Участник оффлайн! plantago
Постоянный участник



 прочитанное сообщение Сообщение на английском  08.03.2021 02:54     Сообщение для модератора  Сообщение для куратора темы       Личное письмо  Отправить e-mail  Web-адрес

https://www.google.com/search?q=R+package+daytime+length

Всего благодарностей: 1Поблагодарили (1): ИНО

*




Кнопка "Транслит" перекодирует
текст из транслита в кирилицу.
Правила перекодировки здесь;
текст в квадратных скобках'[]'
не преобразуется.
Имя:

 преобразовывать смайлики · показать смайлики
Назначение кнопок:

   Поблагодарить автора сообщения — поблагодарить автора
   Удалить сообщение — удалить
   Редактировать сообщение — редактировать
   Поместить сообщение в колонку новостей — поместить в колонку новостей
   Цитировать — цитировать сообщение
   не входит в цитирование/входит в цитирование — цитировать несколько
   Отметить СПАМ-сообщение — обозначить спам
   Сообщение для модератора — связь с модератором
   Участник онлайн!/Участник оффлайн! — автор онлайн/оффлайн
   Фотография — фотография автора

   - остальные обозначения -
 
   *
« Предыдущая тема · Биофизика и матметоды в биологии · Следующая тема »
Быстрый ответДобавить сообщение в темуСоздать новую тему

Rambler   molbiol.ru - методы, информация и программы для молекулярных биологов              

 ·  Викимарт - все интернет-магазины в одном месте  ·  Доска объявлений Board.com.ua  · 
--- сервер арендован в компании Hetzner Online, Германия ---
--- администрирование сервера: Intervipnet ---

Хеликон · Диаэм · ИнтерЛабСервис · Beckman Coulter · SkyGen · ОПТЭК · BIOCAD · Евроген · Синтол · БиоЛайн · Sartorius · Химэксперт · СибЭнзим · Tecan · Даниес · НПП "ТРИС" · Биалекса · ФизЛабПрибор · Genotek · АТГ Сервис Ген · Биоген-Аналитика
Ваш форум  ·  redactor@molbiol.ru  ·  реклама  ·  Дата и время: 21.10.21 12:48
Bridged By IpbWiki: Integration Of Invision Power Board and MediaWiki © GlobalSoft