Страницы

Поиск по вопросам

Показаны сообщения с ярлыком r. Показать все сообщения
Показаны сообщения с ярлыком r. Показать все сообщения

суббота, 11 апреля 2020 г.

Извлечение сплайнов из модели класса GAM (`mgcv::gam`)

#r

                    
Прошу подсказать по следующему вопросу, который является продолжением данного вопроса.
Строим аддитивную модель, как показано на примере из справки ?predict.gam:

 library(mgcv)
 n <- 200
 sig <- 2
 dat <- gamSim(1,n=n,scale=sig)

 b <- gam(y ~ s(x0) + s(I(x1^2)) + s(x2) + offset(x3), data = dat)

 newd <- data.frame(x0=(0:30)/30, x1=(0:30)/30, x2=(0:30)/30, x3=(0:30)/30)

 Xp <- predict(b, newd, type="lpmatrix")

 ##################################################################
 ## The following shows how to use use an "lpmatrix" as a lookup 
 ## table for approximate prediction. The idea is to create 
 ## approximate prediction matrix rows by appropriate linear 
 ## interpolation of an existing prediction matrix. The additivity 
 ## of a GAM makes this possible. 
 ## There is no reason to ever do this in R, but the following 
 ## code provides a useful template for predicting from a fitted 
 ## gam *outside* R: all that is needed is the coefficient vector 
 ## and the prediction matrix. Use larger `Xp'/ smaller `dx' and/or 
 ## higher order interpolation for higher accuracy.  
 ###################################################################

 xn <- c(.341,.122,.476,.981) ## want prediction at these values
 x0 <- 1         ## intercept column
 dx <- 1/30      ## covariate spacing in `newd'
 for (j in 0:2) { ## loop through smooth terms
   cols <- 1+j*9 +1:9      ## relevant cols of Xp
   i <- floor(xn[j+1]*30)  ## find relevant rows of Xp
   w1 <- (xn[j+1]-i*dx)/dx ## interpolation weights
   ## find approx. predict matrix row portion, by interpolation
   x0 <- c(x0,Xp[i+2,cols]*w1 + Xp[i+1,cols]*(1-w1))
 }
 dim(x0)<-c(1,28) 
 fv <- x0%*%coef(b) + xn[4];fv    ## evaluate and add offset
 se <- sqrt(x0%*%b$Vp%*%t(x0));se ## get standard error
 ## compare to normal prediction
 predict(b,newdata=data.frame(x0=xn[1],x1=xn[2],
         x2=xn[3],x3=xn[4]),se=TRUE)


Возможно ли каким-то образом извлечь из модели в явном виде саму формулу сплайна
s(x), которая бы соответствовала форме записи:

yi = β0 + β1 b1 (xi) + β2 b2 (xi) + · · · + βK+3 bK+3(xi) 

(James G. et al. - An Introduction to Statistical Learning with Applications in R)

Заметил, что путем умножения коэффициентов coef(mod) на predict(mod, type="lpmatrix")
(матрица модели) можно получить предсказанные значения, возвращаемые функцией predict(mod,
type="response"). Собственно, проблема в том, что я не могу поставить эти коэффициенты
в соответствие степеням переменной x и (x-knot), как это предполагается при записи
сплайна в виде комбинации базисных функций. 
Приведенный пример предсказания значений на основе интерполяции значений из "lpmatrix"
служит, согласно справке, для использования модели вне среды R. Означает ли это, что
данная реализация GAM не предполагает получения записи модели в явном виде? Есть ли
отличия в этом плане в функции gam() из пакета gam, написанного Хасти и Тибширани -
создателями методам обобщенных аддитивных моделей?
Спасибо.
    


Ответы

Ответ 1



В пакете rms, который является приложением к известной книге Ф. Харрелла, есть функция Function(), которая выдает уравнение поданной на нее модели. К сожалению, объекты класса gam эта функция не принимает. Но в состав rms входят другие функции, которые позволяют подгонять сплайн-модели - возможно, они подойдут и для Ваших целей. Примеры можно посмотреть здесь и здесь.

воскресенье, 15 марта 2020 г.

Подшить к 1му массиву непропущенные значения со 2го в R

#r


Основной массив

df1 <- data.frame(id = c(1,2,3,4,5), dig = c(2,3,NA,5,NA), let = c("a",NA,"c","g",NA))

  id dig  let
1  1   2    a
2  2   3 
3  3  NA    c
4  4   5    g
5  5  NA 


Массив с новыми значениями

df2 <- data.frame(id = c(2,3,5), dig = c(NA,100,200), let = c("letter1",NA,"letter2"))

  id dig     let
1  2  NA letter1
2  3 100    
3  5 200 letter2


Нужно по id подшить непустые значения из df2. То есть, результат должен выглядеть так:

  id dig     let
1  1   2       a
2  2   3 letter1
3  3 100       c
4  4   5       g
5  5 200 letter2

    


Ответы

Ответ 1



repl <- which(is.na(df1[1:3, ]), arr.ind=T) df1[repl] <- df2[repl] Главное - размерности таблиц (исходной и со значениями на замену) должны совпадать, тут я вручную укоротил df1. UPD Более универсальный вариант от автора вопроса: repl <- which(is.na(df1[df1$id %in% df2$id, ]), arr.ind=T) df1[df1$id %in% df2$id, ][repl] <- df2[repl]

четверг, 13 февраля 2020 г.

Множественный классификатор на несбалансированных данных

#r #статистика #машинное_обучение


Нужно построить множественный классификатор (5 классов) на сильно несбалансированной
выборке.

> table(d$class)
    0   0.3   0.5   0.7     1 
12385   736   733    25  1869 


Если просто запустить RandomForest, то ничего получается.
Значит, надо как-то её балансировать. А вот этого я делать-то и не умею.
Все найденные пакеты по балансировке предполагают лишь бинарную классификацию.
Смотрел пакеты:


unbalanced
ROSE


Читал эту статью.

Думал, может в самом RandomForest есть возможность задать cost или сэмплинг — тоже
не нашлось. Как можно решить эту проблему?
    


Ответы

Ответ 1



В пакете "caret" существуют две функции: upSample() и downSample(), которые решают эту проблему. Балансировать классы необходимо независимо от применяемых методов классификации. Обновление Если отвечать широко,то нужно указать, что балансировка классов один из многих важных этапов, которые нужно выполнить прежде чем начать обучать модель. Просто перечислю: выбор и оценка входных переменных, разделение на тренировочную и тестовую выборку(желательно стратифицированную), балансировка классов(только тренировочного набора), препроцессинг (нормализация,стандартизация и т.л.), перемешивание тренировочного набора и др. От качества проведения этих работ на 80% зависит качество получаемого результата моделирования. Если отвечать так широко потребуется статья хорошего объема. Необходимость балансировки классов подтверждается многочисленными экспериментами (не только моими) с многими моделями. Просто сравните результаты классификации одинаковых наборов с и без балансировки Вы убедитесь в этом сами..

Ответ 2



Вот такая функция решает эту проблему my.strata <- function(v) { tmp <- as.vector(table(v)); num_clases <- length(tmp); min_size <- tmp[order(tmp,decreasing=FALSE)[1]]; rep(min_size,num_clases); } randomForest(.... sampsize=my.strata( ) ... )

среда, 12 февраля 2020 г.

синхронизация трех датафреймов по времени

#r


Есть три датафрейма немного разной длинны потому что наблюдения велись начиная с
разного времени,

как их можно синхронизировать по времени чтоб оставить только те наблюдения которые
есть во всех трех фреймах и выкинуть те которые попадаются только в отдельных фреймах 

вот сами дата фреймы

> head(sec1)
        date  time   open   high    low  close vol
1 2016.09.06 08:45 3081.5 3082.5 3080.5 3080.5   6
2 2016.09.06 08:50 3081.5 3081.5 3079.5 3080.5   6
3 2016.09.06 08:55 3081.5 3082.5 3081.5 3082.5  19
4 2016.09.06 09:00 3083.5 3083.5 3081.5 3082.5  19
5 2016.09.06 09:05 3083.5 3085.5 3082.5 3085.5   8
6 2016.09.06 09:10 3086.5 3086.5 3084.5 3086.5  15
> head(sec2)
        date  time  open  high   low close vol
1 2016.09.13 13:00 95.34 95.40 95.33 95.39  36
2 2016.09.13 13:05 95.40 95.43 95.39 95.41  40
3 2016.09.13 13:10 95.42 95.44 95.40 95.42  37
4 2016.09.13 13:15 95.41 95.42 95.39 95.39  25
5 2016.09.13 13:20 95.40 95.41 95.38 95.38  21
6 2016.09.13 13:25 95.39 95.42 95.38 95.42  32
> head(sec3)
        date  time    open    high     low   close vol
1 2016.09.14 18:10 1.12433 1.12456 1.12431 1.12450 137
2 2016.09.14 18:15 1.12444 1.12459 1.12424 1.12455 139
3 2016.09.14 18:20 1.12454 1.12477 1.12446 1.12469 148
4 2016.09.14 18:25 1.12468 1.12474 1.12442 1.12453 120
5 2016.09.14 18:30 1.12452 1.12483 1.12442 1.12482 156
6 2016.09.14 18:35 1.12481 1.12499 1.12472 1.12474 126


Те на выходе должно получиться три датафрейма одинаковой длинны (nrow) и все строчки
датафреймов должны иметь одинаковую дату и время
    


Ответы

Ответ 1



Если я правильно понял задачу, то нужно определить пересекающиеся интервалы дат и времени и отфильтровать наблюдения, попадающие в эти интервалы. Омечу, что приведённые в качестве примера данные не пересекаются по датам. Определим границы для дат: min_date <- list(df1, df2, df3) %>% sapply(. %>% .subset2("date") %>% as.Date(format = "%Y.%m.%d") %>% min()) %>% max() max_date <- list(df1, df2, df3) %>% sapply(. %>% .subset2("date") %>% as.Date(format = "%Y.%m.%d") %>% max()) %>% min() Теперь то же самое для времени: min_time <- list(df1, df2, df3) %>% sapply(. %>% .subset2("time") %>% as.POSIXct(format = "%H:%M") %>% min()) %>% max() max_time <- list(df1, df2, df3) %>% sapply(. %>% .subset2("time") %>% as.POSIXct(format = "%H:%M") %>% min()) %>% min() Теперь можно отфильтровать наблюдения: df1 <- df1 %>% mutate(date = as.Date(date, format = "%Y.%m.%d")) %>% filter(date >= min_date & date <= max_date) %>% mutate(time = as.POSIXct(time, format = "%H:%M")) %>% filter(time >= min_time & time <= max_time) Чтобы код работал, нужно загрузить пакет dplyr.

Ответ 2



Насколько я понимаю, задача сводится к тому, чтобы оставить в каждом датасете только те наблюдения, для которых есть наблюдения с аналогичными значениями date и time в двух других датасетах. Мне видится самым простым решением такое: слить вместе все три датасета сгруппировать наблюдения по переменным date и time посчитать количество наблюдений в группах оставить только те сочетания date и time, которые встречаются 3 раза отфильтровать исходные датасеты по датасету пересекающихся наблюдений Код (не проверял - лень генерировать исходные датасеты; напишите, если где-то что-то упустил, и не работает) library(tidyverse) df_cross <- bind_rows(df1, df2, df3) %>% group_by(date,time) %>% summarise(occurance = n()) %>% ungroup() %>% filter(occurance == 3) %>% select(-occurance) df1_refined <- left_join(df_cross, df1, by = c('date', 'time')) UPD все оказалось еще проще df_cross <- intersect(df1 %>% select(date,time), df2 %>% select(date,time), df3 %>% select(date,time)) df1_refined <- left_join(df_cross, df1, by = c('date', 'time'))

понедельник, 3 февраля 2020 г.

Редактирование дата-фрейма, содержащего NA

#циклы #r


Проблема заключается в следующем:

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

В случае с вектором данная проблема решается прогонкой цикла по вектору, однако для
редактирования дата-фрейма данный способ не подходит.
    


Ответы

Ответ 1



В пакете zoo есть функция na.locf которая делает именно это: > df <- data.frame(a=c(1,NA,2,NA,NA), b=c(1.3,NA,2.4,NA,1.1)) > df a b 1 1 1.3 2 NA NA 3 2 2.4 4 NA NA 5 NA 1.1 > na.locf(df) a b 1 1 1.3 2 1 1.3 3 2 2.4 4 2 2.4 5 2 1.1

пятница, 31 января 2020 г.

Как перенести таблицу в удобный для публикации формат?

#r #markdown


Что делать если нужно не просто получить результат, но и презентовать его?

Я знаю про rmarkdown, но как оформить именно таблицу?
Поделитесь, кто как с этим справляется.
    


Ответы

Ответ 1



Для таблицы самый очевидный вариант - kable(). По другим вариантам из ответа выше: xtable и pander имеют методы для печати различных объектов, например, возвращаемых функциями lm(), t.test и многими другими. Это очень удобно и полезно, и kable() тут не справится. LaTeX изучать не обязательно, вывод в html также работает. xtable с R Markdown используется следующим образом: print(xtable(fit1), type = "html") При этом нужно не забыть указать в чанке ```{r, results='asis'} pander еще проще: pander(fit1). Есть еще относительно экзотический способ, полезный при создании инфографики или при подготовке материалов для печати в солидном журнале: https://github.com/baptiste/gridextra/wiki/tableGrob Позволяет получать таблицы в виде красивых картинок, а также комбинировать их с графиками, например

Ответ 2



Для оформления таблиц есть много пакетов. Эти пакеты позволяют генерировать код в форматах LaTeX, pandoc, HTML. Более менее систематизированный обзор можно посмотреть на CRAN Task Views в разделе Reproducible Research. Ниже приведу краткий список пакетов, которые предоставляют функции для работы с таблицами. knitr (функция kable()) htmlTable xtable pander Я работал с каждым из них и могу подтвердить их работоспособность. Также могут быть полезны пакеты broom, xtable и pander, которые предоставляет функции для оформления типичных объектов R, получаемых в ходе статистического анализа (например, результаты регрессионного или дисперсионного анализа). Посмотреть поддерживаемые пакетом методы можно следующим образом (аналогично для других пакетов): library(pander) methods(pander) Пример вывода kable() в формате markdown: knitr::kable(head(mtcars), format = "markdown") #> | | mpg| cyl| disp| hp| drat| wt| qsec| vs| am| gear| carb| #> |:-----------------|----:|---:|----:|---:|----:|-----:|-----:|--:|--:|----:|----:| #> |Mazda RX4 | 21.0| 6| 160| 110| 3.90| 2.620| 16.46| 0| 1| 4| 4| #> |Mazda RX4 Wag | 21.0| 6| 160| 110| 3.90| 2.875| 17.02| 0| 1| 4| 4| #> |Datsun 710 | 22.8| 4| 108| 93| 3.85| 2.320| 18.61| 1| 1| 4| 1| #> |Hornet 4 Drive | 21.4| 6| 258| 110| 3.08| 3.215| 19.44| 1| 0| 3| 1| #> |Hornet Sportabout | 18.7| 8| 360| 175| 3.15| 3.440| 17.02| 0| 0| 3| 2| #> |Valiant | 18.1| 6| 225| 105| 2.76| 3.460| 20.22| 1| 0| 3| 1|

пятница, 24 января 2020 г.

Сравнение данных с вводом в R и Python

#python #r #конвертация


Доброго времени суток. Учу python ну и собственно хочу перевести код с R и очень
интересует следующий вопрос: Как удобнее всего было бы реализовать сравнение данных
из таблицы excel с критериями, которые вводятся в скрипте, как это сделано в скрипте R?

    a <- as.numeric(readline(prompt="Введите количество рабочих на вашем предприятии:
1, если меньше 8,5 тыс.чел, 2 в противном случае: "))
    b<-as.numeric(readline(prompt="Введите тип вашей индустрии: 2, если высокодоходная,
1 в противном случае: "))
    c<-as.numeric(readline(prompt="Введите производительность труда вашего предприятия
в тыс. руб./чел: "))
    d<-as.numeric(readline(prompt="Введите рентабельность компании в %: "))
    e<-as.numeric(readline(prompt="Введите темп роста компании в %: "))
    tryCatch(
      { localenv <- environment()
      asde<-work_file[as.numeric(work_file$WORKER)==a & as.numeric(work_file$OTRASL)==b
& (work_file$PROISVOD>c-500 & work_file$PROISVODd/100-0.2
& work_file$RENTABe/100-0.3 & work_file$TEMP


Ответы

Ответ 1



Начать можно с такого варианта: import pandas as pd url = r'd:/download/data.xlsx' # читаем Excel в Pandas DataFrame df = pd.read_excel(url) # это нужно будет переделать на ввод текста # или можно читать это из другого CSV/Excel файла a = 1 b = 1 c = 1000 d = 50 e = 80 # query ... qry = ''' WORKER == @a & \ OTRASL == @b & \ PROISVOD > @c-500 & PROISVOD < @c+500 & \ RENTAB > @d/100-0.2 & RENTAB < @d/100+0.2 & \ TEMP > @e/100-0.3 & TEMP < @e/100+0.3 \ ''' print(df.query(qry)) Output: КОМП НОМЕР REAL TEMP RENTAB PROISVOD OTRASL WORKER 25 Башнефть 18 36948.0 0.833 0.329 838.1 1 1 112 Норильск 6 134617.0 0.878 0.394 1402.1 1 1

Как парсить вакансии с API Зарплата.ру в R

#json #парсер #api #r


Есть API джоб-сайта "Зарплата.ру"

https://api.zp.ru/v1/


Как парсить вакансии в R с пакетом jsonlite? То есть, как именно нужно обратиться?
К примеру, через API HeadHunter это корректно делается так:

  string <-"https://api.hh.ru/vacancies?text=\"'machine+learning\"&page="
for (pageNum in 0:5){ # Всего страниц
  data <- fromJSON(paste0(string,  pageNum))
  vacanciesdf <- rbind(vacanciesdf, data.frame(
    data$items$area$name, # Город
    data$items$salary$currency, # Валюта
    data$items$salary$from, # Минимальная оплата
    data$items$employer$name, # Название компании
    data$items$name,#Название должности
    data$items$snippet$requirement)) # Требуемые навыки
  print(paste0("Upload pages:", pageNum + 1))
  Sys.sleep(3)
}


Как решить аналогичную задачу через API "Зарплата.ру" чтобы была возможность задать
ключевое слово и рассортировать данные по столбцам data.frame?
    


Ответы

Ответ 1



Пример кода для выполнения запроса с использованием пакета crul. # HHTP клиент cl <- crul::HttpClient$new(url = "https://api.zp.ru") # Запрос к API resp <- cl$get(path = "v1/vacancies", query = list(scope = "public", q = "machine+learning", limit = 100L)) # Парсинг ответа ans <- jsonlite::fromJSON(resp$parse(encoding = "UTF-8")) # Количество записей в результате выдачи cat(ans$metadata$resultset$count) #> 2 # Извлекаем необходимые поля res <- ans$vacancies data.frame( header = res$header, published_at = as.Date(res$publication$published_at), salary = res$salary, education = res$education$title, experience_length = res$experience_length$title, schedule = res$schedule$title, working_type = res$working_type$title, requirements = res$requirements, url = paste0("https://www.zp.ru", res$url), company = res$company$title, address = paste(res$address$city$title, res$address$street, res$address$building) ) #> header published_at salary education #> 1 Senior, Middle Data scientist 2017-12-08 договорная высшее #> 2 Junior Data scientist 2017-12-08 договорная высшее #> experience_length schedule working_type #> 1 3-5 лет гибкий график полная занятость #> 2 без опыта гибкий график полная занятость #> requirements #> 1 Высшее образование, стаж работы 3-5 лет, полная занятость #> 2 Высшее образование, без опыта, полная занятость #> url #> 1 https://www.zp.ru/vacancy/Senior_Middle_Data_scientist?id=139080429 #> 2 https://www.zp.ru/vacancy/Junior_Data_scientist?id=139080474 #> company address #> 1 СКБ Контур Екатеринбург Малопрудная 5 #> 2 СКБ Контур Екатеринбург Малопрудная 5 Для выгрузки всех результатов в случае если их больше 100, необходимо использовать параметр offset. Описание полей возвращаемого результата приведено в документации по API. При необходимости можно извлекать только нужные поля при помощи параметра fields. Например: query = list(scope = "public", q = "machine+learning", limit = 100L, fields = "header,company.title") Пример постраничной выгрузки: #' @title Функция для выгрузки вакансий с сайта zp.ru #' @param cl HTTP клиент. Создаётся при помощи `crul::HttpClient`. #' @param query Строка, содержащая запрос. #' @param limit Целое число от 1 до 100, определяющее количество результатов. #' @return data.frame с результатами запроса fetch_vacancies <- function(cl, query) { limit <- 100L q <- list( scope = "public", q = query, limit = limit ) fetch_data <- function(query) { # Запрос к API resp <- cl$get(path = "v1/vacancies", query = query) # Проверка статуса ответа resp$raise_for_status() # Парсинг ответа jsonlite::fromJSON(resp$parse(encoding = "UTF-8")) } extract_data <- function(data) { res <- data$vacancies data.frame( header = res$header, published_at = res$publication$published_at, salary = res$salary, education = res$education$title, experience_length = res$experience_length$title, schedule = res$schedule$title, working_type = res$working_type$title, requirements = res$requirements, url = paste0("https://www.zp.ru", res$url), company = res$company$title, address = paste(res$address$city$title, res$address$street, res$address$building) ) } ans <- fetch_data(q) res <- extract_data(ans) # Орабатываем случай, если результатов больше 100 if (ans$metadata$resultset$count > limit) { e <- new.env() # Доабвляем уже полученные данные e[["0"]] <- res offset <- 101L count <- ans$metadata$resultset$count # Повторяем запросы и парсинг со смещением в 100 while (offset < count) { q$offset <- offset res <- extract_data(fetch_data(q)) e[[as.character(offset)]] <- res # Предотвращаем спам запросов Sys.sleep(0.4) # Выводим сообщение cat("\rFetch page ", (offset - 1L) / 100L) offset <- offset + 100L } # Собираем все результаты res <- do.call(rbind, as.list(e)) } # Добавляем обработанный запрос в атрибуты рзультата attr(res, "query") <- ans$metadata$query$searched_q return(res) } # HHTP клиент cl <- crul::HttpClient$new(url = "https://api.zp.ru") query <- "Аналитик" res <- fetch_vacancies(cl, query)

вторник, 7 января 2020 г.

Сохранение пропорциональности графиков в R (radial.pie)

#r #graph


Не в первой сталкиваюсь с такой проблемой, но решения пока не нашел.

Я строю четыре графика отображающих процентное соотношение значений четырех переменных. 

library(plotrix)

x1 <- c(0.09,0.187,0.269,0.441,0.013)
x2 <- c(0.207,0.262,0.259,0.242,0.03)
x3 <- c(0.147,0.237,0.249,0.339,0.027)
x4 <- c(0.052,0.11,0.242,0.581,0.015)

df <- cbind(x1, x2, x3, x4)
df1 <- data.frame(df)
df1 <- sqrt(df1)    #Беру корень по значениям, чтобы сектора выглядели более наглядно.

draw_rp <- function(dataf){
  for (i in 1:length(dataf)){
    xnew1 <- c(0.1, 0.1, 0.1, 0.1, 0.1, dataf[[i]])
    colors <- c(rep('white', nrow(dataf)), rainbow(nrow(dataf)))
    radial.pie(xnew1, show.grid = F, clockwise = T, sector.colors = colors)
  }
}

par(mfrow=c(2, 2))
draw_rp(df1)


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

df1
#         x1        x2        x3        x4
#1 0.3000000 0.4549725 0.3834058 0.2280351
#2 0.4324350 0.5118594 0.4868265 0.3316625
#3 0.5186521 0.5089204 0.4989990 0.4919350
#4 0.6640783 0.4919350 0.5822371 0.7622336
#5 0.1140175 0.1732051 0.1643168 0.1224745




Один из выходов это экспортировать графики в картинки. Оставить как есть график,
у которого самое большое значение и пропорционально размерам его секторов уменьшать
все остальные графики. В этом случае 4й график (0.7622336) оставить как есть, а картинки
остальных подгонять под него.

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


Ответы

Ответ 1



Попробуй добавить в белый сектор max значение из таблицы. ( мне кажется получается то что ты хочешь) draw_rp <- function(dataf){ for (i in 1:length(dataf)){ xnew1 <- c(max(dataf), 0.1, 0.1, 0.1, 0.1, dataf[[i]]) colors <- c(rep('white', nrow(dataf)), rainbow(nrow(dataf))) radial.pie(xnew1, show.grid = F, clockwise = T, sector.colors = colors) } }

понедельник, 6 января 2020 г.

Как в объекте типа data.frame отобрать переменные только одного типа?

#r #dataframe


Например, имеется большой массив данных >100 переменных. А мне необходимо отобрать
 лишь количественные. Как можно сделать подобное?
    


Ответы

Ответ 1



Предварительно стоит изучить структуру данных с помощью функции str(). Выполнить поставленную Вами задачу можно следующим образом: sapply(DF, is) # классы столбцов DF[, sapply(DF, is.numeric)] # все столбцы класса numeric DF[, sapply(DF, is.factor)] # все столбцы класса factor DF[, sapply(DF, is.character)] # все столбцы класса character Можно также использовать комбинацию sapply + which: DF[, which(sapply(DF, is.numeric))]. В некоторых случаях данный вариант показывает более высокую производительность. Также стоит отметить, что для числовых переменных могут применяться разные классы: integer (целые числа), double (числа с плавающей запятой) или numeric (включает в себя два предыдущих типа). Это может пригодиться, если, например, необходимо отфильтровать только столбцы, содержащие целые числа.

Ответ 2



df_numeric <- df[ , sapply(df, is.numeric)] Создание нового дата фрейма df_numeric только с количественными данными из исходного df. is.numeric проверяет столбцы на предмет того, являются ли они количественными

пятница, 27 декабря 2019 г.

R округляет time до даты

#postgresql #r


При загрузке таблицы из PostgreSQL (RPostgreSQL) R округляет до даты.

Вместо 2015-01-28 03:04:01 CET имею 2015-01-28 CET, но Class "POSIXct" "POSIXt",
не только показывает, но и на самом деле 2015-01-28 00:00:00 CET.

Причем это только на локальном маке, R нa сервере получает ту же таблицу без проблем.

Скорее всего, какие-то настройки. Может кто помочь?

> command3 <- "SELECT requested_at FROM rides  WHERE  city_id != 1;"
> riders3 <- dbSendQuery(con, command3)
> riders_total <- fetch(riders3, n = -1)
> riders_total$requested_at[1]
[1] "2015-04-19 CEST"     #####    is "2015-04-19 03:04:31 CEST" !
> riders_total$requested_at[1] + 1
[1] "2015-04-19 00:00:01 CEST"

    


Ответы

Ответ 1



Надо попробовать проверить time zone при чтении. Например base::format(dbGetQuery(con, command3)$requested_at, format="%Z") Также нужно быть уверенным, что в PostgreSQL тип данных http://www.postgresql.org/docs/9.4/static/datatype-datetime.html#DATATYPE-TIMEZONES.

Ответ 2



Проверьте, для начала, одинаковые ли версии R и RSQLite установлены на сервере (под Linux'ом?) и на Mac'е. Лучше везде обновиться до последних версий.

Поиск неповторяющихся значений

#математика #r


Необходимо вывести индексы нулей, но чтобы в строке и в столбце было не более одного
значения. Эта картинка расскажет лучше:



То, что зеленое - нужно вывести индексы. 
То, что красное - не нужно.

Необходимый ответ: 

1-4; 2-1; 3-6; 4-2; 5-7; 6-3; 7-5


Получилось добиться лишь такого:which(tabl_w==0, arr.ind = T)

      row col
 [1,]   2   1
 [2,]   4   1
 [3,]   5   1
 [4,]   6   1
 [5,]   2   2
 [6,]   4   2
 [7,]   6   3
 [8,]   1   4
 [9,]   7   4
[10,]   7   5
[11,]   3   6
[12,]   1   7
[13,]   5   7


А как дальше выбрать индексы, я не знаю.

В итоге получилось сделать:

zero_index <- which(tabl_w==0, arr.ind = T)

choose_unique_zero <- function(x,i=0,j=0,index = data.frame(row=numeric(),col=numeric())){
  zero_temp <- x[(!x[,1] %in% i) & (!x[,2] %in% j),]
  if(length(zero_temp)>2){
    i <- c(i,as.vector(zero_temp[1,1]))
    j <- c(j,as.vector(zero_temp[1,2]))
    index <- rbind(index,setNames(as.list(zero_temp[1,]), names(index)))
    choose_unique_zero(zero_temp,i,j,index)}
  else
    rbind(index,zero_temp)
}

unique_index <- choose_unique_zero(zero_index)


Спасибо всем огромное кто отозвался на мой вопрос!
    


Ответы

Ответ 1



Раз уж вопрос изначально был с тэгом r, то привожу решение с использованием R. Логика решения аналогична коду на PHP. find_ind_r <- function(x, value = 0L) { # Опередляем индексы idx <- which(x == value, arr.ind = TRUE) # Выделяем результирующий объект res <- list( row = rep(NA_integer_, nrow(x)), col = rep(NA_integer_, ncol(x)) ) # Счётчик count <- 1L # Считаем не повторяющиейся индексы for (i in seq_len(nrow(idx))) { r <- idx[i, ] if (!(r[1] %in% res$row) && !(r[2] %in% res$col)) { res$row[count] <- r[1] res$col[count] <- r[2] count <- count + 1L } } # Убираем незаполненные элементы length(res$row) <- count - 1L length(res$col) <- count - 1L return(res) } Если предполагается работа с большими массивами и требуется высокая производительность, то лучше реализовать алгоритм на C++. Ниже приведено решение с использованием Rcpp. Файл test.cpp: #include using namespace Rcpp; // [[Rcpp::export]] List find_ind_cpp(const NumericMatrix & x, double value) { std::size_t ncols = x.ncol(), nrows = x.nrow(); typedef std::unordered_set ind_set; ind_set rows, cols; ind_set::iterator rows_end = rows.end(), cols_end = cols.end(); for (std::size_t i = 0; i < nrows; ++i) { for (std::size_t j = 0; j < ncols; ++j) { if (x[i + nrows * j] == value) { if (rows.find(i + 1) == rows_end && cols.find(j + 1) == cols_end) { rows.insert(i + 1); cols.insert(j + 1); } } } } List res = List::create(rows, cols); res.attr("names") = CharacterVector::create("rows", "cols"); return res; } Выполняется в R-сессии или скрипте: x <- c(4,6,3,0,1,5,0, 0,0,4,3,3,1,12, 8,1,3,3,3,0,2, 0,0,10,5,3,5,12, 0,2,9,8,9,13,0, 0,15,0,1,2,6,10, 2,4,6,0,0,12,1) m <- matrix(x, nrow = 7, ncol = 7, byrow = TRUE) Rcpp::sourceCpp("/tmp/test.cpp") res <- find_ind_r(m, 0) paste(res$row, res$col, sep = "-", collapse = "; ") #> [1] "2-1; 4-2; 6-3; 1-4; 7-5; 3-6; 5-7" Небольшой бенчмарк: set.seed(42) create_mat <- function(n) { matrix(sample(0:15, size = n * n, replace = TRUE), nrow = n, ncol = n) } res <- bench::press( n = c(10, 100, 1000), { mat <- create_mat(n) value <- 0 bench:::mark( min_iterations = 100, find_ind_cpp(mat, value), find_ind_r(mat, value), check = FALSE ) } ) plot(res)

Ответ 2



Вроде это решает ваш вопрос: $arr_1) { foreach ($arr_1 as $key_2 => $val) { $v_tmp++; if ($val == 0 && !in_array($v_tmp, $v) && !in_array($h_tmp, $h)) { $h[] = $h_tmp; $v[] = $v_tmp; $result[] = array('key1' => $key_1, 'key2' => $key_2); } } $v_tmp = 0; $h_tmp++; } var_dump($result); ?> Однако начальные индексы тут начинаются с 0. Вы писали правильные ответы 1-4; 2-1, но на деле будет 0-3; 1-0. Если вам нужно именно 1-4, то просто сделайте инкремент тут: $result[] = array('key1' => $key_1+1, 'key2' => $key_2+1); Наш результат: Array ( [0] => Array ( [key1] => 0 [key2] => 3 ) [1] => Array ( [key1] => 1 [key2] => 0 ) [2] => Array ( [key1] => 2 [key2] => 5 ) [3] => Array ( [key1] => 3 [key2] => 1 ) [4] => Array ( [key1] => 4 [key2] => 6 ) [5] => Array ( [key1] => 5 [key2] => 2 ) [6] => Array ( [key1] => 6 [key2] => 4 ) )

Ответ 3



Задача неясно поставлена. Куда-то пропал критерий оптимальности, если он вообще был. Если цель состоит только в том, чтобы получить не более одного значения в каждой строке и в каждом столбце, то задача решается тривиально: просто хватай нули наугад и следи за тем, чтобы соблюдалось требование "не более одного". Все. Более того, можно вообще взять только один, первый попавшийся ноль и на этом остановиться. Условию задачи это удовлетворяет. По-видимому, требуется найти наибольший набор нулей, удовлетворяющих данному требованию. В идеале - в количестве, равном размеру квадратной матрицы (если она всегда квадратна). Именно это вы, возможно, и хотели показать своим "то, что зеленое". В таком виде это классическая задача на нахождение Максимального Паросочетания в двудольном графе. Для каждой строки матрицы заводится вершина в одной доле графа. Для каждого столбца матрицы - вершина в другой доле графа. Вершины соединяются ребром, если на пересечении соответствующих строки и столбца стоит ноль. Для решения этой задачи в двудольном графе существуют эффективные алгоритмы (Хопкрофт-Карп, например).

среда, 25 декабря 2019 г.

Язык программирования R: Есть ли годная литература на русском? [дубликат]

#книги #r


        
             
                
                    
                        
                            This question already has an answer here:
                            
                        
                    
                
                        
                            Книги и учебные ресурсы по языку R
                                
                                    (1 ответ)
                                
                        
                                Закрыт 3 года назад.
            
                    
Где-то год назад узнал от человека, работающего с геномом, про существование такого
крайне удобного для математических и статистических вычислений языка, как R. Возникла
идея попробовать применить его для сугубо практических целей, например -- в разработке
игр, так как он не только годится для расчётов, но и для удобного вывода обработанной
информации.
Однако язык создавался учёными и для учёных, что несколько отразилось на его синтаксисе
и принципах работы. Короче говоря, попытка решить с его помощью первую же задачу, отличающуюся
по сложности от "Hello, World!", закончилась откровенной неспособностью вкурить в документацию.
Если на Хешкоде есть люди, знакомые с этим языком, прошу у вас совета: Есть ли русскоязычные
книги и мануалы по R (Викиучебник я уже смотрел, но там только самые азы, и то отрывочно),
а также что стоит почитать, чтобы без проблем понимать, что означают мудрёные аргументы
тамошних функций (Текущая проблема, например, возникла с генератором случайных чисел,
в заданном диапазоне)?
Зимой будут лекции по R на Курсэре, но чтобы хоть что-то из них вынести, тем более
-- на английском, в любом случае надо понимать хотя бы терминологию.    


Ответы

Ответ 1



Я как-то не особо искал изданий на русском языке по R. Иногда находил небольшие pdf с описанием методов + наглядное применение. В основном справка в самой программе или в интернете дают исчерпывающий ответ на все вопросы. К сожалению, только на английском. Однако, есть место где есть кое-какое описание на русском. Довольно много методов и примеров разобрано. Буду рад, если это Вам поможет.

Создание функции ggplot

#r #ggplot2


Допустим у меня есть данные

data=data.frame(s=c(10,13,17,8),
                pr=c("a","b","a","b"),
                m=c(rep(as.Date('01.01.2015','%d.%m.%Y'),2), rep(as.Date('01.02.2015','%d.%m.%Y'),2)),
                pr2=c("c","d","d","c"))


И я пытаюсь создать функцию которая рисует ggplot в зависимости от col1 - колонка
по которой нужно делать fill.

plot_function=function(col1){...}


Функция должна возвращать ggplot

ggplot(data = data, aes(x = as.factor(m), y = s/2,fill=col1 ,ymax = max(s/2)*1.1)) + 
  geom_bar(position = "dodge", stat="identity") + 
  geom_text(aes(y=s/4,label=paste(round(s/2,3),"%")),position = position_dodge(.9)) + 
  scale_x_discrete(labels = function(x) format(as.Date(x), "%m/%y")) + 
  xlab("m")


где fill переменная внутри функции.

Пробовал сделать нечто подобное используя aes_string() вместо aes() но не знаю как
совместить операции над переменными например s/2 и aes_string().

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

в dplyr я делал вот так

group=function(data,...){

    dat1=group_by_(data,.dots = lazyeval::lazy_dots(...))
    return(dat1)
  }


Заранее спасибо за помощь.
    


Ответы

Ответ 1



На основании комментария @ArtemKlevtsov и свои догадок на этот счет в итоге сделал вот так gr_plot=function(data_12,nm){ i=which(colnames(data_12)==nm) data_12$var=data_12$s/2 data_12$m=as.factor(data_12$m) j=which(colnames(data_12)=="m") k=which(colnames(data_12)=="var") return( ggplot(data = data_12, aes_string(x = names(data_12)[j], y = names(data_12)[k],fill=names(data_12)[i]))+ geom_bar(position = "dodge",stat="identity")+ geom_text(aes(y=s/4,label=paste(round(s/2,3),"%")),position = position_dodge(.9)) + scale_x_discrete(labels = function(x) format(as.Date(x), "%m/%y")) + xlab("m") ) } соответственно достаточно просто нарисовать в различных разрезах gr_plot(data,"pr") gr_plot(data,"pr2") Если у кого то есть идеи как сделать нечто подобное через lazyeval было бы интересно увидеть.

пятница, 20 декабря 2019 г.

регистр слов в R

#r


Есть вектор с словами

 # текст
text <- c("R is a very essential tool for data analysis. While it 
          is regarded as domain specific, it is a very complete programming 
          language. Almost certainly, many people who would benefit from
          using R, do not use it")
# разбиваю текст на вектор сo словами с пом. пакета stringr
text <- unlist(    stringr::str_match_all(text , '\\w+\\b')   )



 text
 [1] "R"           "is"          "a"           "very"        "essential"   "tool"
       "for"        
 [8] "data"        "analysis"    "While"       "it"          "is"          "regarded"
   "as"         
[15] "domain"      "specific"    "it"          "is"          "a"           "very"
       "complete"   
[22] "programming" "language"    "Almost"      "certainly"   "many"        "people"
     "who"        
[29] "would"       "benefit"     "from"        "using"       "R"           "do" 
        "not"        
[36] "use"         "it"    


Я хочу найти в нем слово "using"

text[text=="using"]
[1] "using"


все нормально, все находит
но если немножко изменить регистр 

text[text=="Using"]
character(0)


то слово найти уже не получиться

Вопрос как сделать поиск слов не чувствительным к регистру?
    


Ответы

Ответ 1



Можно использовать функцию grep grep("Using",text,ignore.case=TRUE,value=TRUE) ignore.case=TRUE- игнор регистра value=TRUE - Возврат значения из вектора, а не позиции найденного слова UPD grep будет искать вхождения. Поэтому поиск "it" вернет не совсем верный результат: grep("it",text,ignore.case=TRUE,value=T,useBytes = T) [1] "it" "it" "benefit" "it" в этом случае функция regexpr отработает лучше: match<-regexpr("IT$",text,ignore.case=TRUE) text[match==1] [1] "it" "it" "it"

Ответ 2



Достаточно привести слова в векторе к одному регистру: text <- tolower(text)

понедельник, 16 декабря 2019 г.

Поддержка юникода в Shiny

#кодировка #unicode #кириллица #r


В процессе создания приложения при помощи пакета Shiny столкнулся с неполной поддержкой
юникода (если это можно так назвать) для кириллических символов.
Русский текст, содержащий букву "я", вызывает ошибку

unexpected INCOMPLETE-STRING


Если текст не содержит этой буковки, то все ок.  

Файл сохраняю в юникоде, как описано в Unicode characters in Shiny apps.  

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

Пример кода:

# ui.R

shinyUI(fluidPage(
    titlePanel("Название приложения"),

    sidebarLayout(
        sidebarPanel( "sidebar panel"),
        mainPanel("main panel")
    )
))

# server.R

shinyServer(function(input, output) {
})


sessionInfo():

R version 3.2.0 (2015-04-16)
Platform: i386-w64-mingw32/i386 (32-bit)
Running under: Windows 7 (build 7601) Service Pack 1

locale:
[1] LC_COLLATE=Russian_Russia.1251  LC_CTYPE=Russian_Russia.1251   
[3] LC_MONETARY=Russian_Russia.1251 LC_NUMERIC=C                   
[5] LC_TIME=Russian_Russia.1251    

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] shiny_0.11.1 zoo_1.7-12  

loaded via a namespace (and not attached):
 [1] R6_2.0.1        htmltools_0.2.6 tools_3.2.0     Rcpp_0.11.6     RJSONIO_1.3-0  
 [6] grid_3.2.0      digest_0.6.8    xtable_1.7-4    httpuv_1.3.2    mime_0.3       
[11] lattice_0.20-31

    


Ответы

Ответ 1



Почитал тут: http://anton-pribora.ru/articles/php/locales/ Немного поэкспериментировал. У меня сработало следующее: один раз в текущей сессии выполнил установку хитрой локали: Sys.setlocale("LC_ALL","Russian_Russia.20866") Все, можно запускать Shiny-приложение с буковками "я" R version 3.2.1 (2015-06-18) Platform: i386-w64-mingw32/i386 (32-bit) Running under: Windows 7 (build 7601) Service Pack 1 locale: [1] LC_COLLATE=Russian_Russia.20866 LC_CTYPE=Russian_Russia.20866 [3] LC_MONETARY=Russian_Russia.20866 LC_NUMERIC=C [5] LC_TIME=Russian_Russia.20866

среда, 27 ноября 2019 г.

Книги и учебные ресурсы по языку R


Рекомендуемая литература, курсы и документация по языку R.


Не создавайте новых ответов — редактируйте общий ответ.
Не размещайте ссылки на нелегальный контент, вроде торрент-трекеров.
Старайтесь сохранять разделение по категориям.





  Данный перечень входит в поддерживаемый сообществом Сборник учебных ресурсов по программированию.

    


Ответы

Ответ 1



Литература Русский А.Б. Шипунов, Е.М. Балдин, П.А. Волкова и др.: Наглядная статистика. Используем R! (ISBN: 978-5-97060-094-8) Бесплатная электронная версия (Книга передана в общественное достояние) Исходные тексты LaTeX Роберт И. Кабаков: R в действии. Анализ и визуализация данных на языке R (перевод с английского, ISBN: 978-1-93518-239-9, 978-5-94074-912-7, 978-5-97060-077-1) Дуглас Люк: Анализ сетей (графов) в среде R. Руководство пользователя. Цветное издание. (перевод с английского, ISBN: 978-5-97060-428-1 ) Джеймс Г., Уиттон Д., Хасти Т., Тибширани Р.: Введение в статистическое обучени с примерами на языке R. Цветное издание. (перевод с английского, ISBN: 978-5-97060-293-5 ) Мастицкий С.Э., Шитиков В.К.: Статистический анализ и визуализация данных с помощью R (ISBN: 978-5-97060-301-7) Бесплатная электронная версия (CC BY-NC-SA 4.0) Репозиторий автора на GitHub, с примерами к книге Блог автора Шитиков В.К., Мастицкий С.Э. Классификация, регрессия и другие алгоритмы Data Mining с использованием R Бесплатная электронная версия (CC BY-NC-SA 4.0) Репозиторий автора на GitHub, с примерами к книге English: Список англоязычной литературы на официальном сайте R-project Сетевые ресурсы Документация и онлайн литература Русский: Викиучебник: Язык программирования R Статьи по R на Хабре А. Б. Шипунов, Е. М. Балдин: Анализ данных с R (сборник статей, изданных в журнале Linux Format, электронная версия с сайта автора) Девид Мертц, Бред Хантинг: Статистическое программирование на R (электронный перевод, части 1, 2, 3) Перевод документации по пакетам dplyr, tidyr, data.table (руководство по data.table на bookdown.org) English: CRAN Task View: обзоры пакетов R по определеным направлениям от лучших разработчиков. Официальный мануал от разработчиков Шпаргалки (примеры кода) от RStudio Cookbook - примеры решения распространенных задач Quick-R (Блог) R-Bloggers: сайт, куда стекаются посты из сотен блогов, общее - использование R. The R Inferno: описание того, как не стоит делать в R и как стоит Bookdown: сайт, который хостит книги, написанные с помощью пакета bookdown. Разумеется, первые книги - сплошные сокровища по освоению и использованию R. Awesome R Advanced R by Hadley Wickham R Data Science Tutorials R for Data Science (Garrett Grolemund, Hadley Wickham) Онлайн курсы Русский: Stepic: Анализ данных в R Stepic: Основы программирования на R Репозиторий автора на GitHub, с кодом и данными Stepic: Анализ данных в R. Часть 2 English: Coursera: R Programming (входит в специализацию Data Science) Data Camp: Introduction to R Udacity: Data Analysis with R. Visually Analyze and Summarize Data Sets Code School: Try R Интерактивное обучение swirl: R пакет для интерактивного обучения (см. русскоязычный материал на хабре)

понедельник, 8 июля 2019 г.

Кластерный анализ методом k-means

Здравствуйте. Имеется две случайная величины. Задача попробовать класстеризовать данные с помощью метода k-means. Разбивал на три кластера данные. Результаты меня удивили.

Почему часть данных, принадлежащая 3 кластеру (синее точки) окружена точками из 2 кластера (зеленные точки)? На картинке это можно увидеть в левом нижнем углу если по-хорошему присмотреться.


Ответ

Это нормальная ситуация. То есть не совсем нормальная, но данный метод кластеризации иногда приводит к таким странным результатам. Дело в том, что расстояния считаются от центроидов, а не от соседних точек. Подробнее об этом и других нюансах можно почитать по ссылке http://dungba.org/the-strange-effect-of-k-means/

воскресенье, 30 июня 2019 г.

Извлечение сплайнов из модели класса GAM (`mgcv::gam`)

Прошу подсказать по следующему вопросу, который является продолжением данного вопроса. Строим аддитивную модель, как показано на примере из справки ?predict.gam
library(mgcv) n <- 200 sig <- 2 dat <- gamSim(1,n=n,scale=sig)
b <- gam(y ~ s(x0) + s(I(x1^2)) + s(x2) + offset(x3), data = dat)
newd <- data.frame(x0=(0:30)/30, x1=(0:30)/30, x2=(0:30)/30, x3=(0:30)/30)
Xp <- predict(b, newd, type="lpmatrix")
################################################################## ## The following shows how to use use an "lpmatrix" as a lookup ## table for approximate prediction. The idea is to create ## approximate prediction matrix rows by appropriate linear ## interpolation of an existing prediction matrix. The additivity ## of a GAM makes this possible. ## There is no reason to ever do this in R, but the following ## code provides a useful template for predicting from a fitted ## gam *outside* R: all that is needed is the coefficient vector ## and the prediction matrix. Use larger `Xp'/ smaller `dx' and/or ## higher order interpolation for higher accuracy. ###################################################################
xn <- c(.341,.122,.476,.981) ## want prediction at these values x0 <- 1 ## intercept column dx <- 1/30 ## covariate spacing in `newd' for (j in 0:2) { ## loop through smooth terms cols <- 1+j*9 +1:9 ## relevant cols of Xp i <- floor(xn[j+1]*30) ## find relevant rows of Xp w1 <- (xn[j+1]-i*dx)/dx ## interpolation weights ## find approx. predict matrix row portion, by interpolation x0 <- c(x0,Xp[i+2,cols]*w1 + Xp[i+1,cols]*(1-w1)) } dim(x0)<-c(1,28) fv <- x0%*%coef(b) + xn[4];fv ## evaluate and add offset se <- sqrt(x0%*%b$Vp%*%t(x0));se ## get standard error ## compare to normal prediction predict(b,newdata=data.frame(x0=xn[1],x1=xn[2], x2=xn[3],x3=xn[4]),se=TRUE)
Возможно ли каким-то образом извлечь из модели в явном виде саму формулу сплайна s(x), которая бы соответствовала форме записи:
yi = β0 + β1 b1 (xi) + β2 b2 (xi) + · · · + βK+3 bK+3(xi)
(James G. et al. - An Introduction to Statistical Learning with Applications in R)
Заметил, что путем умножения коэффициентов coef(mod) на predict(mod, type="lpmatrix") (матрица модели) можно получить предсказанные значения, возвращаемые функцией predict(mod, type="response"). Собственно, проблема в том, что я не могу поставить эти коэффициенты в соответствие степеням переменной x и (x-knot), как это предполагается при записи сплайна в виде комбинации базисных функций. Приведенный пример предсказания значений на основе интерполяции значений из "lpmatrix" служит, согласно справке, для использования модели вне среды R. Означает ли это, что данная реализация GAM не предполагает получения записи модели в явном виде? Есть ли отличия в этом плане в функции gam() из пакета gam, написанного Хасти и Тибширани - создателями методам обобщенных аддитивных моделей? Спасибо.


Ответ

В пакете rms, который является приложением к известной книге Ф. Харрелла, есть функция Function(), которая выдает уравнение поданной на нее модели. К сожалению, объекты класса gam эта функция не принимает. Но в состав rms входят другие функции, которые позволяют подгонять сплайн-модели - возможно, они подойдут и для Ваших целей. Примеры можно посмотреть здесь и здесь

вторник, 4 июня 2019 г.

Подшить к 1му массиву непропущенные значения со 2го в R

Основной массив
df1 <- data.frame(id = c(1,2,3,4,5), dig = c(2,3,NA,5,NA), let = c("a",NA,"c","g",NA))
id dig let 1 1 2 a 2 2 3 3 3 NA c 4 4 5 g 5 5 NA
Массив с новыми значениями
df2 <- data.frame(id = c(2,3,5), dig = c(NA,100,200), let = c("letter1",NA,"letter2"))
id dig let 1 2 NA letter1 2 3 100 3 5 200 letter2
Нужно по id подшить непустые значения из df2. То есть, результат должен выглядеть так:
id dig let 1 1 2 a 2 2 3 letter1 3 3 100 c 4 4 5 g 5 5 200 letter2


Ответ

repl <- which(is.na(df1[1:3, ]), arr.ind=T) df1[repl] <- df2[repl]
Главное - размерности таблиц (исходной и со значениями на замену) должны совпадать, тут я вручную укоротил df1. UPD Более универсальный вариант от автора вопроса:
repl <- which(is.na(df1[df1$id %in% df2$id, ]), arr.ind=T) df1[df1$id %in% df2$id, ][repl] <- df2[repl]