 * |   |
|
Molbiol.ru | О проекте | Справочник | Методы | Растворы | Расчёты | Литература | Орг.вопросы Web | Фирмы | Coffee break | Картинки | Работы и услуги | Биржа труда |  Zbio-wikiNG SEQUENCING · ЖИЗНЬ РАСТЕНИЙ · БИОХИМИЯ · ГОРОДСКИЕ КОМАРЫ · А.А.ЛЮБИЩЕВ · ЗООМУЗЕЙ Темы за 24 часа [ Вход* | Регистрация* ] Форум: |
|
Правила
* R Help -- Давайте составим русский FAQ --
Кураторы темы:* plantago
Операции: Хочу стать куратором* · Подписаться на тему* · Отправить страницу по e-mail · Версия для печати*
Внешний вид:* Схема · [ Стандартный ] · +Перв.сообщ.
   |
 PS2004RПостоянный участник  |
 15.04.2020 10:57
(Алекс3212 @ 15.04.2020 00:52)  Да я и сам написал, а потом уже понял, что сам бы не понял, что нужно на месте другого человека....извиняюсь.исходный файл вот такой > d <- read.table("data.txt") > d V1 V2 V3 1 1111 2222 2222 2 1111 2222 2223 3 1111 2222 1111 4 1111 2222 1112 Те значения что нужны min и max содержатся в V3, проблема (для меня пока что), в том, что этот столбец удаляется (оно и ясно иначе не выйдет) для получения файла такого вида X1 X2 X3 X4 x 1 1111 2222 1 1 2 2 1111 2222 1 2 2 Если на реальном примере и используя код который в моем сообщении был то исходник такой V1 <int> V2 <int> V3 <int> delta <dbl> 170509 5337 1545 1 170509 5337 1546 1 170509 5337 1547 1 170509 5337 1548 1 170509 5337 1549 1 Отлично считается количество файлов, которые идут подряд, но V3 столбец исчезает, а именно там содержатся файлы, которые нужны теперь(( min и max значение. И получается фрейм такого вида X1 <fctr> X2 <fctr> X3 <fctr> x <int> 170509 5337 1 160 170509 5385 0 1 170509 5401 1 5 А нужно, чтобы был примерно следующий 170509 5337 1 160 1545 1705 170509 5385 0 1 5385 5385 170509 5401 1 5 5401 5406 т.е. последние 2 столбца это значения минимум и максимум того что суммируется...пока не придумал, как сделать то что нужно мне, но в голове представляется пока прям огромный путь с манипуляцией данными, т.к. придется это делать с использованием каких то пакетов и использовать какие то ухищрения, кажется мне, что с помощью base пакета можно сделать все в пару строчек, но я его плохо знаю, от чего и такой вопрос)) обещаю выложить код с садомазо ухищрениями как придумаю, но не уверен что кому то нужен будет он, повторюсь, мне кажется там будет пару десятков строк. Поэтому буду рад за любую инфу, пусть не готовый код а просто куда копать)) решение простое, нанимаете r-программиста на почасовую ставку, и интерактивно задаете ему вводные, оплачивая каждый этап исполнение по затратам времени программиста, до тех пор пока вас не удовлетворит очередной результат. |
|
plantago Постоянный участник |
15.04.2020 12:37
 но все же попробую  170509 5337 1 160 1545 1705 170509 5385 0 1 5385 5385 170509 5401 1 5 5401 5406 т.е. последние 2 столбца это значения минимум и максимум того что суммируется.. Шестой столбец -- это сумма четверного и пятого, правильно?
|
|
Алекс3212 |
17.04.2020 20:38
(plantago @ 15.04.2020 13:37) В общем-то согласен с PS2004R но все же попробую Шестой столбец -- это сумма четверного и пятого, правильно? Да все верно 4+5 =6)) Проблема вытащить 4 столбец у меня...пока что...если так будет справедливее, пусть решение выложит сюда тот кто сможет, после того как я выложу свой код с костылями..я в любом случае буду благодарен. Насчет программиста я согласен полностью, любой труд должен быть оплачен и я не прошу выкладывать готовый код, прошу просто подсказать где может быть решение, хотя б в какой стороне) Можно было бы и нанять и заплатить, но это не коммерческая задача в целом, и мне нравится думать своей головой, что касается r, хоть иногда я попадаю в тупик из за недостатка знаний в r, все таки я планирую использовать r как можно больше, т.к. это дает преимущества в выполнении разных задач)) |
|
plantago Постоянный участник |
17.04.2020 22:00
> Да все верно ... Тогда нужен только пятый, шестой Вы сами легко посчитаете. > ... мне нравится думать своей головой ... Тогда придумайте, пожалуйста, минимальный воспроизводимый пример, как требует stackoverflow: или Если Вы его придумаете, тогда Вам, скорее всего помогут. Если нет, то вряд ли. |
|
Алекс3212 |
17.04.2020 22:39
(plantago @ 17.04.2020 23:00) >> Шестой столбец -- это сумма четверного и пятого, правильно?> Да все верно ... Тогда нужен только пятый, шестой Вы сами легко посчитаете. > ... мне нравится думать своей головой ... Тогда придумайте, пожалуйста, минимальный воспроизводимый пример, как требует stackoverflow: или Если Вы его придумаете, тогда Вам, скорее всего помогут. Если нет, то вряд ли. Вот готовый код, но я не уверен, что все сделал верно...буду рад если укажите на ошибки или поможете сократить его..всем спасибо. d.y <- read.table("G://badchannel0416.asc") %>% # по ссылке текстовик mutate(prpk = (V1*10000+V2)*10000+V3) %>% arrange(prpk) %>% select(V1, V2, V3) delta.d.y <- c(diff(d.y[,3]), 0) s <- which(!(delta.d.y == 1)) delta.d.y[which(delta.d.y==1)+1] <- 1 delta.d.y[which(!(delta.d.y==1))] <- 0 sec <- rep(1:(length(s)), times=c(s[1],diff(s))) d.y <- data.frame(d.y, delta.d.y, sec) res.y <- rle(paste(d.y[,1], d.y[,2], d.y[,4], d.y[,5])) FR1 <- data.frame(do.call(rbind, strsplit(res.y$values," ")),x=res.y$lengths) FR1 <- FR1 %>% filter(X3 == 1) #------ dfg <- d.y %>% mutate(prpkff = (V2*10000+V3)) dfg$dif <- c(dfg$prpkff[1], diff(dfg$prpkff)) dfg01 <- dfg %>% mutate(dif2 = ifelse(dif != 1, 0, 1)) dfg02 <- dfg01 %>% mutate(row1 = 1:nrow(dfg01)) %>% mutate(minV3 =ifelse(dif2 == 0, V3, 0)) dfg_first <- d.y %>% mutate(row1 = 1:nrow(d.y)) dfg_full <- full_join(dfg_first, dfg02, by=c("row1"="row1")) %>% select(V1.x, V2.x, V3.x, delta.d.y.x, minV3) res1 <- rle(paste(dfg_full[,1], dfg_full[,2], dfg_full[,4], dfg_full[,5])) FFF1 <- data.frame(do.call(rbind, strsplit(res1$values," ")),x=res1$lengths) FFF1$X4 <- as.numeric(FFF1$X4) FFF1$x <- as.numeric(FFF1$x) FFF1 <- FFF1 %>% filter(X3 ==1) %>% filter(X4 != 0) FR1$min <- FFF1$X4 FR2 <- FR1 %>% mutate(max = (x + min)-1) %>% select(-X3, -X4) %>% filter(x >= 10) Файл не прикрепляется, загрузил в облако Сообщение было отредактировано Алекс3212 - 18.04.2020 01:41 |
|
plantago Постоянный участник |
18.04.2020 10:15
Error in `$<-.data.frame`(`*tmp*`, min, value = c(146, 1, 8, 1, 18, 1, : replacement has 7497 rows, data has 3765 так и должно быть?
|
|
Алекс3212 |
18.04.2020 17:46
(plantago @ 18.04.2020 11:15) > FR1$min <- FFF1$X4Error in `$<-.data.frame`(`*tmp*`, min, value = c(146, 1, 8, 1, 18, 1, : replacement has 7497 rows, data has 3765 так и должно быть? Нет, поменял немного на вот так d.y <- read.table("G:/badchannel0416.asc") %>% mutate(prpk = (V1*10000+V2)*10000+V3) %>% arrange(prpk) %>% select(V1, V2, V3) delta.d.y <- c(diff(d.y[,3]), 0) s <- which(!(delta.d.y == 1)) delta.d.y[which(delta.d.y==1)+1] <- 1 delta.d.y[which(!(delta.d.y==1))] <- 0 sec <- rep(1:(length(s)), times=c(s[1],diff(s))) d.y <- data.frame(d.y, delta.d.y, sec) res.y <- rle(paste(d.y[,1], d.y[,2], d.y[,4], d.y[,5])) FR1 <- data.frame(do.call(rbind, strsplit(res.y$values," ")),x=res.y$lengths) FR1 <- FR1 %>% filter(X3 == 1) #------ dfg <- d.y %>% mutate(prpkff = (V2*10000+V3)) dfg$dif <- c(dfg$prpkff[1], diff(dfg$prpkff)) dfg01 <- dfg %>% mutate(dif2 = ifelse(dif != 1, 0, 1)) dfg02 <- dfg01 %>% mutate(row1 = 1:nrow(dfg01)) %>% mutate(minV3 =ifelse(dif2 == 0, V3, 0)) dfg_first <- d.y %>% mutate(row1 = 1:nrow(d.y)) dfg_full <- full_join(dfg_first, dfg02, by=c("row1"="row1")) %>% select(V1.x, V2.x, V3.x, delta.d.y.x, minV3) res1 <- rle(paste(dfg_full[,1], dfg_full[,2], dfg_full[,4], dfg_full[,5])) FFF1 <- data.frame(do.call(rbind, strsplit(res1$values," ")),x=res1$lengths) FFF1$X4 <- as.character(FFF1$X4) FFF1$X4 <- as.numeric(FFF1$X4) FFF1$x <- as.numeric(FFF1$x) FFF1 <- FFF1 %>% filter(X3 ==1) %>% filter(X4 > 1) FR1 FFF1 FR1$min <- FFF1$X4 FR2 <- FR1 %>% mutate(max = (x + min)-1) %>% select(-X3, -X4) %>% filter(x >= 10) |
|
plantago Постоянный участник |
19.04.2020 08:56
library(dplyr) вначале. Теперь, пожалуйста, сформулируйте Ваш вопрос _в связи с этим примером_. Какую именно переменную надо сделать, где она должна появиться, откуда ее брать, и т.д.
|
|
Алекс3212 |
19.04.2020 19:34
(plantago @ 19.04.2020 09:56) Теперь работает, только Вы забыли указать library(dplyr) вначале. Теперь, пожалуйста, сформулируйте Ваш вопрос _в связи с этим примером_. Какую именно переменную надо сделать, где она должна появиться, откуда ее брать, и т.д. Да, библиотеку не указал... Это готовый код, который дает нужный мне результат, я обещал выложить, я выложил как только сделал...то, что он длинный и с костылями будет я предупреждал...вопросы только в части сокращения кода, если это возможно...какую переменную откуда и куда понятно мне кажется из кода. Вообще я надеюсь, что это недопонимание между нами)) ато как то странно, выложить тут готовый код, а потом писать что мне нужно помочь в его написании...он ведь готовый как бы, вроде работает...изначально да, был другой вопрос, но разобрался сам... Есть другие вопросы по другим задачам, но честно сказать я уже примерно знаю что услышу, поэтому буду сам пытаться, чтобы не нервировать никого)) |
|
plantago Постоянный участник |
20.04.2020 07:26
> Проблема вытащить 4 столбец у меня > разобрался сам... Так Вы уже решили эту проблему самостоятельно?
|
|
plantago Постоянный участник |
22.04.2020 21:10
Что касается улучшения скрипта, то здесь лучше применить принцип "работает -- не трогай". Но если есть какие-то _конкретные_ вопросы, задавайте.
|
|
-Эн- |
09.07.2020 03:41
|
|
plantago Постоянный участник |
09.07.2020 05:56
|
|
plantago Постоянный участник |
09.07.2020 11:59
|
|
ПолинаШ Участник |
09.07.2020 20:57
Хочется на картосхеме показать пространственное распределение показателя Z , построенное по выборке точек с нерегулярно заданными координатами. Для этого надо выполнить интерполяцию на равноточечный грид. Начала с неспециализированых методов loess и gam: eastseq <- seq(44.532, 55.648, by = .2) northseq <- seq(48.65536, 55.11243, by = .2) R.grid <- expand.grid(Долгота = eastseq, Широта = northseq) ### Сглаживаеие loess m <- loess(z ~ Долгота * Широта, span = 0.85,degree = 2, data = df) R.fit <- predict(m, R.grid) image(eastseq, northseq, R.fit, col=colorRampPalette(c("white", "gray10"))(100), xlab="Долгота",ylab="Широта") contour(eastseq, northseq, R.fit, nlevels=10 , col="green", add=T) ### Сглаживаеие gam require("mgcv") mod <- gam(z ~ te(Долгота , Широта), data = df) G.fit <- matrix(predict(mod, R.grid), ncol = length(northseq)) image(eastseq, northseq, G.fit, col=colorRampPalette(c("white", "gray10"))(100), xlab="Долгота",ylab="Широта") contour(eastseq, northseq, G.fit, nlevels=10 , col="green", add=T) Получила на редкость отвратительные картинки Вопрос легко решился использованием интерполяции из пакета gstat. Но осадочек остался: может я что не так делаю. Ведь в пакете vegan функция ordisurf как раз проводит сглаживание сплайнами, обращаясь к gam. Хочу посмотреть, как это делается. Но как посмотреть R-код этого метода: UseMethod("ordisurf") |
|
plantago Постоянный участник |
10.07.2020 07:24
Ну вот кто мешает приложить этот самый df, чтобы и другие могли посмотреть на картинки и понять, как Вам можно помочь? > Вопрос легко решился использованием интерполяции из пакета gstat. И этот код тоже было бы хорошо привести, ведь тогда будет лучше понятно, что именно Вам нужно. И другим читателям форума полезно. > Но как посмотреть R-код этого метода === > library(vegan) Loading required package: permute Loading required package: lattice This is vegan 2.5-6 > UseMethod("ordisurf") Error in UseMethod("ordisurf") : 'UseMethod' used in an inappropriate fashion > methods("ordisurf") [1] ordisurf.default* ordisurf.formula* see '?methods' for accessing help and source code > getAnywhere(ordisurf.default) A single object matching ‘ordisurf.default’ was found It was found in the following places registered S3 method for ordisurf from namespace vegan namespace:vegan with value function (x, y, choices = c(1, 2), knots = 10, family = "gaussian", ... === Вот так А можно скачать исходники отсюда Исходники Открыть их, найти файл "ordisurf.R" и посмотреть содержимое, на этот раз в авторском оформлении и, что важно, с авторскими комментариями в коде. |
|
ПолинаШ Участник |
10.07.2020 17:55
Еще чуток поразбираюсь и все выложу в конфу
|
|
ПолинаШ Участник |
07.08.2020 15:53
Вот как я получал прилагаемые рисунки: df1 <- read.table("PRODIAMESINAE.txt",sep ="\t") df1 <- df[df$z > 0,] eastseq <- seq(44.532, 55.648, by = .2) northseq <- seq(48.65536, 55.11243, by = .2) R.grid <- expand.grid(Долгота = eastseq, Широта = northseq) nrow(R.grid) 1848 ### Сглаживаеие четырьмя различными методами m <- loess(z ~ Долгота * Широта, span = 0.85,degree = 2, data = df) R.fit <- predict(m, R.grid) image(eastseq, northseq, R.fit, col=colorRampPalette(c("white", "gray10"))(100), xlab="Долгота",ylab="Широта") contour(eastseq, northseq, R.fit, nlevels=10 , col="green", add=T) points(df1$Долгота,df1$Широта, pch=20, col="yellow") require("mgcv") mod <- gam(z ~ te(Долгота , Широта), data = df) G.fit <- matrix(predict(mod, R.grid), ncol = length(northseq)) image(eastseq, northseq, G.fit, col=colorRampPalette(c("white", "gray10"))(100), xlab="Долгота",ylab="Широта") contour(eastseq, northseq, G.fit, nlevels=10 , col="green", add=T) points(df1$Долгота,df1$Широта, pch=20, col="yellow") library(akima) akima.spl <- interp(df$Долгота,df$Широта,df$z,eastseq,northseq, linear=FALSE) image(akima.spl, col=colorRampPalette(c("white", "gray10"))(100), xlab="Долгота",ylab="Широта") contour(akima.spl, nlevels=10 , col="green", add=T) points(df1$Долгота,df1$Широта, pch=20, col="yellow") library(sp) library(gstat) df_xy <- df colnames(df_xy) <- c("x","y","z") coordinates(df_xy) = ~x + y xy.grid <- expand.grid(x = eastseq, y = northseq) coordinates(xy.grid) <- ~x + y gridded(xy.grid) <- TRUE idw <- idw(formula = z ~ 1, locations = df_xy, newdata = xy.grid) idw.output = as.data.frame(idw) names(idw.output)[1:3] <- c("Долгота", "Широта", "z") library(ggplot2) # start needed libraries library(metR) My_map_cont + ggplot() + geom_tile(data = idw.output, aes(x = Долгота, y = Широта, fill=z), alpha = "0.75") + scale_fill_gradient(low = "white", high = "green") + geom_point(data=df1, aes(x=Долгота, y=Широта), shape=21, colour="red") + geom_contour(data = idw.output,aes(x = Долгота, y = Широта, z=z), bins = 7) + theme_bw() Вроде бы мне понятно, что в данном случае использование геостатистики естественно. Но есть много задач, не относящихся к геоинформатике. Как тут быть? Картинки:  loes1.jpg — (61.46к) Файл/ы:
|
PRODIAMESINAE.txt|
plantago Постоянный участник |
08.08.2020 00:32
Но в целом я так и не понял, в чем проблема. Поправьте, пожалуйста, скрипт, чтобы он работал и объясните, что именно нужно. |
|
ПолинаШ Участник |
10.08.2020 08:51
Просто уберите единичку в первой строке (извините, недоглядела из-за разных экспериментов, а таблицу грузила из другого места): df <- read.table("PRODIAMESINAE.txt",sep ="\t") и закомментируйте My_map_cont + (это подгрузка карты подложки, здесь не нужной). И все заработает без каких-либо дополнительных усилий  Но, еще раз простите за невнимательность). А проблема в том, что сглаживание LOES дает совершенно абсурдный результат, сглаживание GAM - чуть получше, но тоже получается несъедобная поверхность, Аkima - в целом верно, но эстетически скучно, да и пропорции исказились... Сообщение было отредактировано ПолинаШ - 10.08.2020 09:01 |
|
PS2004R Постоянный участник |
10.08.2020 10:16
(ПолинаШ @ 10.08.2020 08:51) А проблема в том, что сглаживание LOES дает совершенно абсурдный результат, сглаживание GAM - чуть получше, но тоже получается несъедобная поверхность, Аkima - в целом верно, но эстетически скучно, да и пропорции исказились... Это называется "кригинг". Он проводиться опираясь на тот или иной вариант вариограммы. Если вариограмма построена чисто формально, то результат бессмысленный. Начинайте с вариограммы, доказывайте (да хоть бутстрепом) что она хоть что то реально отражает. И потом кригингом "сглаживайте" (интерполируйте-экстраполируйте). Сообщение было отредактировано PS2004R - 11.08.2020 09:14 |
|
ПолинаШ Участник |
11.08.2020 18:16
Но есть много задач, не относящихся к геоинформатике. Например, построить поверхность зависимости концентрации озона от температуры и скорости ветра. В частности, в функции ordisurf() из Veganа используют сглаживание сплайнами через gam() и никаких вариограмм не строят. В принципе, нельзя отбрасывать и возможность трехмерного сглаживания LOES, так хорошо зарекомендовавшее себя в двумерном случае. Но у меня ничего не получилось и я просто спросила: это что, принципиальная "ущербность" сглаживающих функций в трехмерном случае, или я что-то не так делала. Но, как говорят в нашей богадельне, на нет и судна нет. 
|
|
PS2004R Постоянный участник |
11.08.2020 18:27
(ПолинаШ @ 11.08.2020 18:16) Как я писала выше, мне понятно, что в данном случае использование геостатистики естественно и не спрашивала о различных вариантах кригинга - там есть своя теория и ее надо использовать. Но есть много задач, не относящихся к геоинформатике. Например, построить поверхность зависимости концентрации озона от температуры и скорости ветра. В частности, в функции ordisurf() из Veganа используют сглаживание сплайнами через gam() и никаких вариограмм не строят. В принципе, нельзя отбрасывать и возможность трехмерного сглаживания LOES, так хорошо зарекомендовавшее себя в двумерном случае. Но у меня ничего не получилось и я просто спросила: это что, принципиальная "ущербность" сглаживающих функций в трехмерном случае, или я что-то не так делала. Но, как говорят в нашей богадельне, на нет и судна нет. Я ответил на ваш вопрос ранее, другого ответа для вас у меня увы нет. |
|
plantago Постоянный участник |
12.08.2020 17:04
Вот исправленный вариант на случай, если кто-то хочет выполнить его без больших проблем. (Только мне еще пришлось поставить 11 пакетов). CODE df <- read.table("PRODIAMESINAE.txt",sep ="\t", fileEncoding="windows-1251") ## corrected !! names(df) <- c("Dolgota", "SHirota", "z") ## added !! df1 <- df[df$z > 0,] eastseq <- seq(44.532, 55.648, by = .2) northseq <- seq(48.65536, 55.11243, by = .2) R.grid <- expand.grid(Dolgota = eastseq, SHirota = northseq) nrow(R.grid) 1848 ### Sglazhivaeie chetihrjmya razlichnihmi metodami m <- loess(z ~ Dolgota * SHirota, span = 0.85,degree = 2, data = df) R.fit <- predict(m, R.grid) image(eastseq, northseq, R.fit, col=colorRampPalette(c("white", "gray10"))(100), xlab="Dolgota",ylab="SHirota") contour(eastseq, northseq, R.fit, nlevels=10 , col="green", add=T) points(df1$Dolgota,df1$SHirota, pch=20, col="yellow") require("mgcv") mod <- gam(z ~ te(Dolgota , SHirota), data = df) G.fit <- matrix(predict(mod, R.grid), ncol = length(northseq)) image(eastseq, northseq, G.fit, col=colorRampPalette(c("white", "gray10"))(100), xlab="Dolgota",ylab="SHirota") contour(eastseq, northseq, G.fit, nlevels=10 , col="green", add=T) points(df1$Dolgota,df1$SHirota, pch=20, col="yellow") library(akima) akima.spl <- interp(df$Dolgota,df$SHirota,df$z,eastseq,northseq, linear=FALSE) image(akima.spl, col=colorRampPalette(c("white", "gray10"))(100), xlab="Dolgota",ylab="SHirota") contour(akima.spl, nlevels=10 , col="green", add=T) points(df1$Dolgota,df1$SHirota, pch=20, col="yellow") library(sp) library(gstat) df_xy <- df colnames(df_xy) <- c("x","y","z") coordinates(df_xy) = ~x + y xy.grid <- expand.grid(x = eastseq, y = northseq) coordinates(xy.grid) <- ~x + y gridded(xy.grid) <- TRUE idw <- idw(formula = z ~ 1, locations = df_xy, newdata = xy.grid) idw.output = as.data.frame(idw) names(idw.output)[1:3] <- c("Dolgota", "SHirota", "z") library(ggplot2) # start needed libraries library(metR) #My_map_cont + ## corrected !! ggplot() + geom_tile(data = idw.output, aes(x = Dolgota, y = SHirota, fill=z), alpha = 0.75) + ## corrected !! scale_fill_gradient(low = "white", high = "green") + geom_point(data=df1, aes(x=Dolgota, y=SHirota), shape=21, colour="red") + geom_contour(data = idw.output,aes(x = Dolgota, y = SHirota, z=z), bins = 7) + theme_bw() (данные "PRODIAMESINAE.txt" нужно скачать выше по треду).
|
|
plantago Постоянный участник |
12.08.2020 17:08
Да, и данные никак не объяснены. Что это за данные? Для чего проводится анализ? И что такое "z"? |
|
Den-N Постоянный участник |
13.08.2020 20:36
(ПолинаШ @ 11.08.2020 21:16) ...Но есть много задач, не относящихся к геоинформатике. Например, построить поверхность зависимости концентрации озона от температуры и скорости ветра. В частности, в функции ordisurf() из Veganа используют сглаживание сплайнами через gam() и никаких вариограмм не строят. В принципе, нельзя отбрасывать и возможность трехмерного сглаживания LOES, так хорошо зарекомендовавшее себя в двумерном случае. Но у меня ничего не получилось и я просто спросила: это что, принципиальная "ущербность" сглаживающих функций в трехмерном случае, или я что-то не так делала. Но, как говорят в нашей богадельне, на нет и судна нет. Нужные вам негеостатистические интерполирующие функции есть и их много. Действительно, геостатистика - специфическое направление и рассматривается некоторыми авторами сродни искусству. Например, какую функцию выбрать для построения вариаграммы, учитывать ли и как анизотропию, какой метод кригинга использовать? - многовато субъективизма. И раньше, и сейчас более жёсткие и простые интерполяторы поверхностей широко распространены и даже встраиваются в некоторые ГИС. Помимо метода обратных расстояний (inverse distance weighting interpolation, IDW) это: триангуляция (triangulation), многоуровневое сглаживание сплайнами (multilevel B-spline approximation, MBA), радиальная интерполяция (radial basis function interpolation - для округлых объектов), метод ближайшего соседа (nearest neighbour interpolation), может ещё что-то более экзотичное. Все эти функции есть в R в качестве самостоятельных пакетов и/или входят в крупные. Также читал про симуляционные процедуры, многократно интерполирующие с использованием части данных и затем усредняющие результат для каждой ячейки грида - в R это можно реализовать для любого интерполятора и получить не только предметную карту по средним, но и карту неопределённости по дисперсиям. Из подобных методов за простоту и результативность мне нравятся MBA и наименьшего искривления (minimum curvature, не нашёл в R). Но в MBA, как и в любых сплайнах есть проблема излишней волнистости, поэтому я обрадовался увидев в этой ветке так любимые мной в 2D сплайны Акимы. В 3D я пользовался ими только раз в технической задаче (функциональная зависимость практически без шума) и тогда они показали себя лучше LOESS. Но конкретно для вашего набора данных они всё-таки ушли в минуса подобно всяким сплайнам (на юге центральной части). Думаю, можно искусственно приравнять отрицательные значения в наборе нулю, но это уже не вполне comme il faut. У меня получилось подобно вашей картинке, только я шаг поменьше сделал, вышло более гладко (по умолчанию akima делит область по горизонтали и вертикали на 40 частей), а также использовал другую цветовую палитру. #c русскими осями не заморачивался (просто x, y, z) т.к. обычно довожу и объединяю рисунки из R в векторном редакторе TpX. library(akima) ak1 <- with(df,interp(x,y,z, linear=FALSE, nx=500, ny=500)) library(RColorBrewer) with(ak1,image(x,y,z, col=rev(brewer.pal(n=10, name="Spectral")))) points(df$x,df$y, pch=20, col="black") C GAM в 3D экспериментировать пока некогда, но в 2D я остался ей очень доволен... Сообщение было отредактировано Den-N - 13.08.2020 21:30 Картинки:  akima.png — (16.73к)
|
|
PS2004R Постоянный участник |
13.08.2020 23:41
Да и вариограмма вполне себе "перспективно" выглядит (в смысле "видали и похужее" ) ). Картинки:  _________________2020_08_13_23_39_48.png — (21.45к)  _________________2020_08_13_23_32_06.png — (30.87к) |
|
ПолинаШ Участник |
15.08.2020 11:09
Но в моем скрипте формальных ошибок нет (кроме злополучной df1), у меня на Win7 все работает. Если у Вас проблемы с кириллицей - я тут не причем. PS2004R - Вы прекрасно отвечали, но, к сожалению, совсем не на мой вопрос. Den-N - большущее спасибо. Все прекрасно объяснили и вернули обсуждение в конструктивное русло. Спасибо и за наводку на MBA. Прекрасный пакет и думаю, что проблему излишней волнистости можно как-то решить подбором параметров. Вот что у меня получилось и меня это удовлетворило: df <- read.table("PRODIAMESINAE.txt",sep ="\t") df1 <- df[df$z > 0,] library(MBA) mba.int <- mba.surf(df, 300, 300, extend=TRUE)$xyz.est ##Image plot image(mba.int, xaxs="r", yaxs="r") contour(mba.int, nlevels=10 , col="green", add=T) points(df1$Долгота,df1$Широта, pch=20, col="blue") А Акима что-то смещает результат относительно системы координат Сообщение было отредактировано ПолинаШ - 15.08.2020 11:26 Картинки:  MBA.png — (40.41к)
|
|
ИНО Постоянный участник Донецк |
15.08.2020 11:33
Как пример, решается задача распределения численности личинок комара-звонца по некой территории. Это стохастическая система, потому строить поверхность, точно проходящую через все точки выборки неверно. Тут регрессия нужна.
|
|
PS2004R Постоянный участник |
15.08.2020 12:03
(ПолинаШ @ 15.08.2020 11:09) plantago PS2004R - Вы прекрасно отвечали, но, к сожалению, совсем не на мой вопрос. Мало ли что вам там показалось. Лучше займитесь оценкой модели выборочной вариограммы своих данных, без этого все эти "сглаживания" ничтожны с точки зрения статистики. PS Меня вообще всё это "ненавязчивое натягивание совы на глобус" жутко раздражает, донатягивались уже до края по моему. |
|
PS2004R Постоянный участник |
15.08.2020 12:32
Берем произвольный метод сглаживания и перебираем "параметры гладкости" в нем до момента когда функция "качества подгонки" достигает оптимума. Функция "качество подгонки" это почти обычная кроссвалидация, только с учетом пространственной природы данных. Исключаемый в кроссвалидации участок в этом случае должен быть "связанный регион". Вот посмотреть на пространство параметров подгонки заполенное качеством подгонки и найти область где есть оптимум качества. Но зачем весь этот "закат солнца вручную" делать, если есть _ясный_ и _понятный_ способ с вариограммой в которой все понятно и который по сути делает все тоже самое? PS если же хочется некого surface reconstruction то берите и делайте. Это не статистическая процедура получения гладкой поверхности. PPS сайт лежит, поэтому сразу отсюда Сообщение было отредактировано PS2004R - 15.08.2020 12:43 |
|
ПолинаШ Участник |
16.08.2020 15:01
Я же задала очень локальный и простой вопрос: как правильно построить интерполирующую 3-мерную поверхность с использованием функций loess() и gam(). Исключительно в методических целях. Например, loess() дает глубокую впадину посередине участка, тогда как функции idw(), akima(), mba.surf() да и ваш кригинг (не знаю, какую функцию Вы использовали) дали вполне нормальную идентичную картину с характерным трехгорбием. Вероятно, стоило бы поиграться со span и degree, но предикция с использованием параметров по умолчанию так шокировала, что захотелось спросить... И прошу прощения за беспокойство
|
|
PS2004R Постоянный участник |
16.08.2020 20:14
(ПолинаШ @ 16.08.2020 15:01) Однако саму оптимальность (в том числе, и кросс-проверкой) можно оценить только косвенным путем. Поскольку самих критериев оптимальности можно придумать много. 1. Тут прямой критерий (прямее не бывает), "насколько совпадает прогноз всей модели с фактическими данными исключенной из обучения области экстраполяции-интерполяции". Лично я затрудняюсь придумать более общий критерий качества проверяющий модель по выборочным данным (ну если принять, что кроме выборки ничего больше нет). Если что то есть, то с удовольствием выслушаю. 2. Простая интерполяция в ваших данных невозможна как решение, наиболее близко к высказанному желанию получить интерполирующую поверхность в 3д пространстве это упомянутый surface reconstruction который принципиально допускает ошибку измерения в данных (и даже грубые выбросы). И имеет вмонтированный в модель приор о существовании некой "поверхности". |
|
Guest IP-штамп: fretc2UtzYVf6 гость |
16.08.2020 23:48
(ПолинаШ @ 16.08.2020 15:01) ...Я же задала очень локальный и простой вопрос: как правильно построить интерполирующую 3-мерную поверхность с использованием функций loess() и gam(). Исключительно в методических целях. Попробовал gam, тоже в растерянности: очень много настроек. Вы использовали взаимодействие координат типа te, а там есть ещё три, причём каждый со своими настройками "по умолчанию", которые можно менять. С наскока не понял, какое решение всё-таки ближе к искомому "простому" интерполятору  mod2d <- gam(z~s(x, y), data = df, method = "REML") vis.gam(mod2d, view = c("x", "y"), n.grid=300, plot.type = "contour", color="topo", nCol=50, too.far = 0.25) points(df$x,df$y, pch=20, col="black")
|
|
Den-N Постоянный участник |
16.08.2020 23:50
Картинки:  gam4.png — (187.5к)
|
|
plantago Постоянный участник |
17.08.2020 08:05
(Почти) Оккамовское решение: === library(hexbin) pp <- read.table("PRODIAMESINAE.txt", h=TRUE) zz <- ifelse(pp$z == 0, 1, round(pp$z)) plot(hexbin(pp[rep(1:nrow(pp), times=zz), -3])) ===  (можно обойтись без round(), но rep() использует as.integer(), а не round(), и вообще это ключевое место лучше описывать явно)
|
|
ПолинаШ Участник |
17.08.2020 16:12
PS2004R - иногда важна не только точное и устойчивое воспроизведение моделью выборочных данных (оцениваемое кросс-проверкой), но и отсутствие аномалий поверхности (в первую очередь, неприятные горбы или впадины по краям). Den-N Ага!, а еще можно и gam(z ~ te(x) + te(y), ... и сплайнов штук 6 типов, да и разные типы взять для x и y, да и кнотом поварьировать. Комбинаторики намного хватит. plantago - Все правильно. Все зависит от задачи. "Предсказание", "Объяснение", "Глобальный тренд", "Локальный тренд", "Изоляция расстоянием" и проч. Здесь я к проблеме отнеслась еще проще. Аппроксимируем плоскостью данные z относительно географических координат xy: z = -45.19 - 0.0495x + 0.9437y Эмпирическое значение тангенса угла наклона плоскости k = 0,9450 Строим распределение угла наклона при справедливости нулевой гипотезы. Для этого формируем 1000 рандомизированных выборок, в каждой из которых значения z случайно перемешаны относительно x и y. Среднее имитируемое 0.268, стандартное отклонение 0.171 и только 3 значения из 1000 оказались больше, чем эмпирическое значение. Таким образом, нулевая гипотеза об отсутствии пространственного тренда отклонена с достигнутым уровнем значимости 0.003 Далее показаем на графике "направление главного удара", т.е. проводим стрелку в направлении, совпадающем с проекцией нормали к полученной плоскости на плоскость XOY. Делаем вывод, что животных больше в северных лесах. НаВашем графике это отчетливо видно. И извините за некоторые недоразумения. Сообщение было отредактировано ПолинаШ - 17.08.2020 16:15
|
|
ПолинаШ Участник |
06.01.2021 12:32
Нужна функция, возвращающая содержимое таблицы по имени. get_df <- function(table_name) get(table_name) library(vegan) сhoices = c("BCI", "dune", "mite", "pyrifos", "sipoo","varespec") data(mite) df <- get_df(сhoices[3]) str(df) 'data.frame': 70 obs. of 35 variables: $ Brachy : int 17 2 4 23 5 19 17 5 3 22 ... $ PHTH : int 5 7 3 7 8 7 3 4 3 4 ... $ HPAV : int 5 16 1 10 13 5 8 8 2 5 ... get_df(сhoices[2]) Ошибка в get(сhoices[2]) :объект 'dune' не найден for (obj in ls()) { if(class(get(obj)) == "data.frame") print(obj) } [1] "df" [1] "mite" Как открыть объект данных, не используя предварительно data() |
|
ПолинаШ Участник |
06.01.2021 15:19
Кнопочная благодарность не срабатывает |
|
Guest IP-штамп: frt6Z9EOlgm.A гость |
06.01.2021 17:44
(ПолинаШ @ 06.01.2021 12:32) Подскажите в довольно простом вопросе.Нужна функция, возвращающая содержимое таблицы по имени. get_df <- function(table_name) get(table_name) library(vegan) сhoices = c("BCI", "dune", "mite", "pyrifos", "sipoo","varespec") data(mite) df <- get_df(сhoices[3]) str(df) 'data.frame': 70 obs. of 35 variables: $ Brachy : int 17 2 4 23 5 19 17 5 3 22 ... $ PHTH : int 5 7 3 7 8 7 3 4 3 4 ... $ HPAV : int 5 16 1 10 13 5 8 8 2 5 ... get_df(сhoices[2]) Ошибка в get(сhoices[2]) :объект 'dune' не найден for (obj in ls()) { if(class(get(obj)) == "data.frame") print(obj) } [1] "df" [1] "mite" Как открыть объект данных, не используя предварительно data() Use ‘getAnywhere’ for searching for an object anywhere, including in other namespaces, and ‘getFromNamespace’ to find an object in a specific namespace. |
|
plantago Постоянный участник |
06.01.2021 18:41
Я сначала написал ответ, в котором создавалось новое окружение (environment), но потом понял, что это не решает проблемы и стер его. Похоже, что надо изменить код таким образом: === get_df <- function(table_name) get(table_name) library(vegan) choices <- c("BCI", "dune", "mite", "pyrifos", "sipoo", "varespec") data(list=choices) get_df(choices[3]) === Без data(), кажется, не обойтись, но зато все эти таблицы можно загрузить скопом. Кстати говоря, в оригинальном коде неприятная проблема: в слове "choices" первая буква кириллическая. |
|
Guest IP-штамп: frt6Z9EOlgm.A гость |
06.01.2021 19:28
(plantago @ 06.01.2021 18:41) Похоже, что get* функции не будут работать. Наверное, потому, что механизм загрузки датасетов очень специфический и к тому же в разных пакетах разный. Например, "vegan" не использует lazy-loading, поэтому после его загрузки 'str(dune)' или 'getAnywhere(dune)' работать не будет, в то время как в других пакетах такое работает.Я сначала написал ответ, в котором создавалось новое окружение (environment), но потом понял, что это не решает проблемы и стер его. Похоже, что надо изменить код таким образом: === get_df <- function(table_name) get(table_name) library(vegan) choices <- c("BCI", "dune", "mite", "pyrifos", "sipoo", "varespec") data(list=choices) get_df(choices[3]) === Без data(), кажется, не обойтись, но зато все эти таблицы можно загрузить скопом. Кстати говоря, в оригинальном коде неприятная проблема: в слове "choices" первая буква кириллическая. Доктор, меня все игнорируют (С) (но в хорошей компании, вместе со всем "?" скопом )))CODE latform: x86_64-pc-linux-gnu (64-bit) R -- это свободное ПО, и оно поставляется безо всяких гарантий. Вы вольны распространять его при соблюдении некоторых условий. Введите 'license()' для получения более подробной информации. R -- это проект, в котором сотрудничает множество разработчиков. Введите 'contributors()' для получения дополнительной информации и 'citation()' для ознакомления с правилами упоминания R и его пакетов в публикациях. Введите 'demo()' для запуска демонстрационных программ, 'help()' -- для получения справки, 'help.start()' -- для доступа к справке через браузер. Введите 'q()', чтобы выйти из R. [Загружено ранее сохраненное рабочее пространство] > get_df <- function(table_name) getAnywhere(table_name) > df <- get_df(сhoices[3]) > str(df) List of 5 $ name : chr "mite" $ objs :List of 1 ..$ .GlobalEnv:'data.frame': 70 obs. of 35 variables: .. ..$ Brachy : int [1:70] 17 2 4 23 5 19 17 5 3 22 ... .. ..$ PHTH : int [1:70] 5 7 3 7 8 7 3 4 3 4 ... .. ..$ HPAV : int [1:70] 5 16 1 10 13 5 8 8 2 5 ... .. ..$ RARD : int [1:70] 3 0 1 2 9 9 2 2 2 3 ... .. ..$ SSTR : int [1:70] 2 6 2 2 0 3 3 1 1 0 ... .. ..$ Protopl : int [1:70] 1 0 0 0 13 2 0 2 1 0 ... .. ..$ MEGR : int [1:70] 4 4 3 4 0 3 3 3 12 0 ... .. ..$ MPRO : int [1:70] 2 2 0 0 0 0 0 0 0 0 ... .. ..$ TVIE : int [1:70] 2 0 0 1 0 0 0 0 0 0 ... .. ..$ HMIN : int [1:70] 1 0 0 2 3 20 19 1 0 11 ... .. ..$ HMIN2 : int [1:70] 4 1 6 10 14 16 3 4 3 4 ... .. ..$ NPRA : int [1:70] 1 3 3 0 3 2 0 10 9 5 ... .. ..$ TVEL : int [1:70] 17 21 20 18 32 13 22 12 20 16 ... .. ..$ ONOV : int [1:70] 4 27 17 47 43 38 27 25 8 33 ... .. ..$ SUCT : int [1:70] 9 12 10 17 27 39 37 26 19 29 ... .. ..$ LCIL : int [1:70] 50 138 89 108 5 3 0 0 30 0 ... .. ..$ Oribatl1: int [1:70] 3 6 3 10 1 5 4 6 7 8 ... .. ..$ Ceratoz1: int [1:70] 1 0 0 1 0 0 0 0 0 2 ... .. ..$ PWIL : int [1:70] 1 1 2 0 5 1 0 3 0 1 ... .. ..$ Galumna1: int [1:70] 8 3 1 1 2 1 1 1 2 5 ... .. ..$ Stgncrs2: int [1:70] 0 9 8 2 1 8 3 3 1 2 ... .. ..$ HRUF : int [1:70] 0 1 0 1 0 0 0 2 0 0 ... .. ..$ Trhypch1: int [1:70] 0 1 3 2 1 4 0 0 0 0 ... .. ..$ PPEL : int [1:70] 0 1 0 1 0 0 0 0 0 0 ... .. ..$ NCOR : int [1:70] 0 2 2 3 0 1 0 0 0 0 ... .. ..$ SLAT : int [1:70] 0 2 0 2 0 0 0 0 0 3 ... .. ..$ FSET : int [1:70] 0 2 8 12 12 10 9 6 0 10 ... .. ..$ Lepidzts: int [1:70] 0 1 0 0 2 0 0 1 0 0 ... .. ..$ Eupelops: int [1:70] 0 0 0 0 0 0 1 0 0 1 ... .. ..$ Miniglmn: int [1:70] 0 0 0 0 0 0 0 1 0 2 ... .. ..$ LRUG : int [1:70] 0 0 0 0 0 0 0 0 0 0 ... .. ..$ PLAG2 : int [1:70] 0 0 0 0 0 0 0 0 0 0 ... .. ..$ Ceratoz3: int [1:70] 0 0 0 0 0 0 0 0 0 0 ... .. ..$ Oppiminu: int [1:70] 0 0 0 0 0 0 0 0 0 0 ... .. ..$ Trimalc2: int [1:70] 0 0 0 0 0 0 0 0 0 0 ... $ where : chr ".GlobalEnv" $ visible: logi TRUE $ dups : logi FALSE - attr(*, "class")= chr "getAnywhere" > |
|
plantago Постоянный участник |
 06.01.2021 21:18
> library(vegan) Loading required package: permute Loading required package: lattice This is vegan 2.5-7 > > choices <- c("BCI", "dune", "mite", "pyrifos", "sipoo", "varespec") > getAnywhere("dune") no object named ‘dune’ was found > getAnywhere(choices[2]) no object named ‘dune’ was found > get_df(choices[2]) no object named ‘dune’ was found > sessionInfo() R version 4.0.3 (2020-10-10) Platform: x86_64-pc-linux-gnu (64-bit) Running under: Ubuntu 16.04.7 LTS Matrix products: default BLAS: /usr/lib/libblas/libblas.so.3.6.0 LAPACK: /usr/lib/lapack/liblapack.so.3.6.0 locale: [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8 [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8 [7] LC_PAPER=en_US.UTF-8 LC_NAME=C [9] LC_ADDRESS=C LC_TELEPHONE=C [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C attached base packages: [1] stats graphics grDevices utils datasets methods base other attached packages: [1] vegan_2.5-7 lattice_0.20-41 permute_0.9-5 colorout_1.2-2 loaded via a namespace (and not attached): [1] MASS_7.3-53 compiler_4.0.3 Matrix_1.3-0 parallel_4.0.3 tools_4.0.3 [6] mgcv_1.8-33 splines_4.0.3 nlme_3.1-151 grid_4.0.3 cluster_2.1.0 |
|
Guest IP-штамп: fr7rsGcRlbyvg гость |
07.01.2021 02:59
CODE > data(package="vegan")$result[,3] [1] "BCI" "BCI.env" "dune" "dune.env" [5] "dune.phylodis" "dune.taxon" "mite" "mite.env" [9] "mite.pcnm" "mite.xy" "pyrifos" "sipoo" [13] "sipoo.map" "varechem" "varespec" dune это имя объекта файловой системы (там в result остальные элементы пути лежат к нему), и data() все это разбирает |
|
ИНО Постоянный участник Донецк |
08.01.2021 22:16
), поэтому прошу помощи. Имеется датафрэйм подобного вида: CODE dataframe1<-data.frame(Pred1=c(1,2,3), Pred2=c(30, 20, 10), RespA=c(0, 1, 2), RespB=c(2, 1, 0)) dataframe1 Pred1 Pred2 RespA RespB 1 1 30 0 2 2 2 20 1 1 3 3 10 2 0 Где последние два столбца - таблица сопряженности. Требуется ее развернуть, чтоб получилось что-то типа этого: 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"))) dataframe2 Pred1 Pred2 Resp 1 1 30 B 2 1 30 B 3 2 20 A 4 2 20 B 5 3 10 A 6 3 10 A Сообщение было отредактировано ИНО - 08.01.2021 22:17 |
|
Guest IP-штамп: fr7rsGcRlbyvg гость |
09.01.2021 01:59
(ИНО @ 08.01.2021 22:16) Простая вроде бы задачка, но с R я далеко не "на ты" (раз в годик встречаемся обычно ), поэтому прошу помощи. Имеется датафрэйм подобного вида:CODE dataframe1<-data.frame(Pred1=c(1,2,3), Pred2=c(30, 20, 10), RespA=c(0, 1, 2), RespB=c(2, 1, 0)) dataframe1 Pred1 Pred2 RespA RespB 1 1 30 0 2 2 2 20 1 1 3 3 10 2 0 Где последние два столбца - таблица сопряженности. Требуется ее развернуть, чтоб получилось что-то типа этого: 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"))) dataframe2 Pred1 Pred2 Resp 1 1 30 B 2 1 30 B 3 2 20 A 4 2 20 B 5 3 10 A 6 3 10 A Непонятный вопрос, но вот так например: CODE > reshape2::melt(dataframe1, id.var=c("Pred1","Pred2")) Pred1 Pred2 variable value 1 1 30 RespA 0 2 2 20 RespA 1 3 3 10 RespA 2 4 1 30 RespB 2 5 2 20 RespB 1 6 3 10 RespB 0 |
|
ИНО Постоянный участник Донецк |
09.01.2021 03:09
Попробую объяснить подробнее: есть данные, собранные в форме первой таблицы (а ля dataframe2). Pred1 и Pred2 - некие независимые переменные A и B - варианты бинарного отклика. Ситуация вполне жизненная: измеряем значения Pred1 и Pred2, подсчитываем количество объектов исследования, находящихся в состоянии A, заносим полученное число в столбец RespA, подсчитываем количество объектов, находящихся в состоянии B - заносим число в столбец RespB. Далее идем в другое место, снова измеряем Pred1 и Pred2 (получаем уже другие значения) и снова подсчитаем объекты находящиеся в каждом их двух альтернативных состояний. И т. д. Пока храним, все хорошо и компактно, но вот вздумалось применить например, логистическую регрессию, чтобы определить влияние Pred1 и Pred2 на вероятность нахождения объекта в состоянии через A или Б, а glm() требует данные в "длинном формате" (а ля dataframe2) , по одной строчке на каждый объект (или же можно и таблицу сопряженности ей как-то скормить?). Естественно, приведенный пример вымышленный, в реальном и предикатов, и наблюдений много больше. Руками переписывать доооолго. Думал, ситуация достаточно обычная
|
|
Guest IP-штамп: fr7rsGcRlbyvg гость |
09.01.2021 15:00
(ИНО @ 09.01.2021 03:09) Увы, не то. Надо, чтобы сток с сочетанием Pred1=1, Pred2=30 и Resp=B должно стало две штуки, а Pred1=1, Pred2=30 и Resp=A - ни одной и т. д. Попробую объяснить подробнее: есть данные, собранные в форме первой таблицы (а ля dataframe2). Pred1 и Pred2 - некие независимые переменные A и B - варианты бинарного отклика. Ситуация вполне жизненная: измеряем значения Pred1 и Pred2, подсчитываем количество объектов исследования, находящихся в состоянии A, заносим полученное число в столбец RespA, подсчитываем количество объектов, находящихся в состоянии B - заносим число в столбец RespB. Далее идем в другое место, снова измеряем Pred1 и Pred2 (получаем уже другие значения) и снова подсчитаем объекты находящиеся в каждом их двух альтернативных состояний. И т. д. Пока храним, все хорошо и компактно, но вот вздумалось применить например, логистическую регрессию, чтобы определить влияние Pred1 и Pred2 на вероятность нахождения объекта в состоянии через A или Б, а glm() требует данные в "длинном формате" (а ля dataframe2) , по одной строчке на каждый объект (или же можно и таблицу сопряженности ей как-то скормить?). Естественно, приведенный пример вымышленный, в реальном и предикатов, и наблюдений много больше. Руками переписывать доооолго. Думал, ситуация достаточно обычная Ну так в последнем столбце число таких повторов и записано. Вот так например CODE > d <- reshape2::melt(dataframe1, id.var=c("Pred1","Pred2")) > do.call(rbind, lapply(apply(d[d[,4]!=0,], 1, function(x) replicate(x[4],x[1:3])),t) ) Pred1 Pred2 variable [1,] "2" "20" "RespA" [2,] "3" "10" "RespA" [3,] "3" "10" "RespA" [4,] "1" "30" "RespB" [5,] "1" "30" "RespB" [6,] "2" "20" "RespB"
|
|
ИНО Постоянный участник Донецк |
09.01.2021 18:11
|
|
plantago Постоянный участник |
09.01.2021 20:05
=== > dataframe1 <- data.frame(Pred1=c(1,2,3), Pred2=c(30, 20, 10), RespA=c(0, 1, 2), RespB=c(2, 1, 0)) > df2 <- rbind(as.matrix(dataframe1[, 1:3]), as.matrix(dataframe1[, c(1:2, 4)])) > df2 <- as.data.frame(df2) > df2$Resp <- rep(c("A", "B"), each=nrow(dataframe1)) > df3 <- df2[rep(1:nrow(df2), times=df2[, 3]), -3] > df3 Pred1 Pred2 Resp 2 2 20 A 3 3 10 A 3.1 3 10 A 4 1 30 B 4.1 1 30 B 5 2 20 B ===
|
| « Предыдущая тема · Биофизика и матметоды в биологии · Следующая тема » |
 |
· Викимарт - все интернет-магазины в одном месте · Доска объявлений Board.com.ua ·
--- сервер арендован в компании Hetzner Online, Германия ---
--- администрирование сервера: Intervipnet ---
Хеликон · Диаэм · ИнтерЛабСервис · Beckman Coulter · SkyGen · ОПТЭК · BIOCAD · Евроген · Синтол · БиоЛайн · Sartorius · Химэксперт · СибЭнзим · Tecan · Даниес · НПП "ТРИС" · Биалекса · ФизЛабПрибор · Genotek · АТГ Сервис Ген · Биоген-Аналитика
Powered by Invision Power Board(Trial) v2.0.0 © 2024 IPS, Inc.