GNU R: Вопросы и ответы

[Hogfather]

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

> head(MyD)
  A B   C     Y
1 1 1 0.1   6.6
2 2 1 0.2   9.2
3 3 2 0.3  32.3
4 4 2 0.4  40.4
5 5 3 0.5  88.0
6 6 3 0.6 105.6

Тогда линейна модель МНК Y от всех без исключения параметров выглядит так:

> summary(mdl<-lm(Y~.,data=MyD))

Call:
lm(formula = Y ~ ., data = MyD)

Residuals:
    Min      1Q  Median      3Q     Max 
-56.588 -26.540  -0.223  20.860  79.223 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept)   -21.69      52.98  -0.409    0.691
A             -76.39      44.15  -1.730    0.114
B              64.29      77.05   0.834    0.424
C             548.60     761.95   0.720    0.488

Residual standard error: 46.28 on 10 degrees of freedom
Multiple R-squared: 0.987,	Adjusted R-squared: 0.9831 
F-statistic: 252.6 on 3 and 10 DF,  p-value: 1.01e-09

Если нам нужны еще произведения параметров, то делаем вот так:

> summary(mdl<-lm(Y~(.)^2,data=MyD))

Call:
lm(formula = Y ~ (.)^2, data = MyD)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.68607 -0.04158 -0.00208  0.11227  0.68607 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)    8.6861     1.4453   6.010 0.000537 ***
A              6.2334     0.6792   9.177 3.76e-05 ***
B             -4.1133     2.1816  -1.885 0.101355    
C           -121.0582     9.2645 -13.067 3.58e-06 ***
A:B            0.4813     0.5852   0.822 0.437927    
A:C           -2.2817     2.3619  -0.966 0.366206    
B:C           79.8960     0.2492 320.614 7.58e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 0.4427 on 7 degrees of freedom
Multiple R-squared:     1,	Adjusted R-squared:     1 
F-statistic: 1.398e+06 on 6 and 7 DF,  p-value: < 2.2e-16

Знаки "*" и "^" имеют в описании модели сакральный смысл.
Результат, аналогичный вышеописанному можно также получить

> summary(mdl<-lm(Y~(.)*(.),data=MyD))

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

> MyD$A2=MyD$A^2
> MyD$AB=MyD$A*MyD$B
> MyD$AC=MyD$A*MyD$C
> MyD$BC=MyD$B*MyD$C
> MyD$B2=MyD$B^2
> MyD$C2=MyD$C^2
> MyD$ABC=MyD$A*MyD$B*MyD$C
> head(MyD)
  A B   C     Y A2 AB  AC  BC B2   C2  ABC
1 1 1 0.1   6.6  1  1 0.1 0.1  1 0.01  0.1
2 2 1 0.2   9.2  4  2 0.4 0.2  1 0.04  0.4
3 3 2 0.3  32.3  9  6 0.9 0.6  4 0.09  1.8
4 4 2 0.4  40.4 16  8 1.6 0.8  4 0.16  3.2
5 5 3 0.5  88.0 25 15 2.5 1.5  9 0.25  7.5
6 6 3 0.6 105.6 36 18 3.6 1.8  9 0.36 10.8

Теперь строим модель, способом показанным выше

> summary(mdl<-lm(Y~.,data=MyD))

Call:
lm(formula = Y ~ ., data = MyD)

Residuals:
       Min         1Q     Median         3Q        Max 
-7.240e-14 -1.505e-14  0.000e+00  2.387e-14  6.347e-14 

Coefficients: (3 not defined because of singularities)
              Estimate Std. Error    t value Pr(>|t|)    
(Intercept) -6.533e-13  5.138e-13 -1.271e+00    0.251    
A            1.725e+01  5.911e-13  2.918e+13   <2e-16 ***
B            2.500e+00  4.752e-13  5.261e+12   <2e-16 ***
C           -1.740e+02  3.122e-12 -5.573e+13   <2e-16 ***
A2          -1.250e+00  6.484e-14 -1.928e+13   <2e-16 ***
AB          -7.663e-14  8.777e-14 -8.730e-01    0.416    
AC                  NA         NA         NA       NA    
BC           5.000e+01  1.619e-12  3.089e+13   <2e-16 ***
B2                  NA         NA         NA       NA    
C2                  NA         NA         NA       NA    
ABC          5.000e+00  2.707e-13  1.847e+13   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 6.341e-14 on 6 degrees of freedom
Multiple R-squared:     1,	Adjusted R-squared:     1 
F-statistic: 5.843e+31 on 7 and 6 DF,  p-value: < 2.2e-16

Чисто чтоб поржать, сделаем автоматический подбор наилучшей модели.

> (newmod<-step(mdl))
Start:  AIC=-846.76
Y ~ A + B + C + A2 + AB + AC + BC + B2 + C2 + ABC


Step:  AIC=-846.76
Y ~ A + B + C + A2 + AB + AC + BC + B2 + ABC


Step:  AIC=-846.76
Y ~ A + B + C + A2 + AB + AC + BC + ABC


Step:  AIC=-846.76
Y ~ A + B + C + A2 + AB + BC + ABC

       Df Sum of Sq     RSS     AIC
- AB    1    0.0000  0.0000 -852.25
<none>               0.0000 -846.76
- B     1    0.1113  0.1113  -53.68
- ABC   1    1.3721  1.3721  -18.52
- A2    1    1.4944  1.4944  -17.32
- A     1    3.4249  3.4249   -5.71
- BC    1    3.8362  3.8362   -4.12
- C     1   12.4892 12.4892   12.40

Step:  AIC=-852.25
Y ~ A + B + C + A2 + BC + ABC

       Df Sum of Sq     RSS     AIC
<none>               0.0000 -852.25
- B     1    0.3263  0.3263  -40.63
- ABC   1    1.5047  1.5047  -19.23
- A2    1    1.5742  1.5742  -18.59
- A     1    3.5163  3.5163   -7.34
- BC    1    4.1983  4.1983   -4.86
- C     1   15.0305 15.0305   12.99

Call:
lm(formula = Y ~ A + B + C + A2 + BC + ABC, data = MyD)

Coefficients:
(Intercept)            A            B            C           A2           BC  
 -3.950e-13    1.725e+01    2.500e+00   -1.740e+02   -1.250e+00    5.000e+01  
        ABC  
  5.000e+00  

Предупреждения
1: попытка выбора модели для в целом хорошей подгонки -- глупость 
2: попытка выбора модели для в целом хорошей подгонки -- глупость 
3: попытка выбора модели для в целом хорошей подгонки -- глупость 
4: попытка выбора модели для в целом хорошей подгонки -- глупость 
5: попытка выбора модели для в целом хорошей подгонки -- глупость 

А теперь "на бис", то же самое, но без создания фиктивных переменных. Поскольку символы "^" и "*" имеют особый смысл, посмотрите как это рекомендуется делать (можно было начать с "." через + добавить недостающие параметры, но я сделал в явном виде)

> (summary(mdl<-lm(Y~A+B+C+A:B+B:C+A:C+A:B:C+I(A^2)+I(B^2)+I(C^2),data=MyD)))

Call:
lm(formula = Y ~ A + B + C + A:B + B:C + A:C + A:B:C + I(A^2) + 
    I(B^2) + I(C^2), data = MyD)

Residuals:
       Min         1Q     Median         3Q        Max 
-7.500e-14 -1.737e-15  0.000e+00  2.102e-15  6.022e-14 

Coefficients: (3 not defined because of singularities)
              Estimate Std. Error    t value Pr(>|t|)    
(Intercept)  5.194e-13  3.875e-13  1.340e+00    0.229    
A            1.000e+00  1.183e-12  8.454e+11   <2e-16 ***
B           -6.000e-13  4.321e-13 -1.389e+00    0.214    
C            1.000e+00  1.052e-11  9.503e+10   <2e-16 ***
I(A^2)       8.061e-14  1.113e-13  7.240e-01    0.496    
I(B^2)       5.000e+00  1.551e-13  3.225e+13   <2e-16 ***
I(C^2)      -6.361e-12  6.620e-12 -9.610e-01    0.374    
A:B                 NA         NA         NA       NA    
B:C                 NA         NA         NA       NA    
A:C                 NA         NA         NA       NA    
A:B:C        5.000e+00  2.042e-13  2.449e+13   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 4.783e-14 on 6 degrees of freedom
Multiple R-squared:     1,	Adjusted R-squared:     1 
F-statistic: 1.027e+32 on 7 and 6 DF,  p-value: < 2.2e-16

Обратите внимание, получилось несколько иное решение, также с Эр Квадрат равное 1. А что Вы хотели, кто обещал, что будет легко?

Теперь о тесте Рамсея RESET.
Естественно, что он есть в GNU R и находится в пакете lmtest.

> library(lmtest)
> x <- c(1:30)
> y1 <- 1 + x + x^2 + rnorm(30)
> y2 <- 1 + x + rnorm(30)
> resettest(y1 ~ x , power=2, type="regressor")

	RESET test

data:  y1 ~ x 
RESET = 108553.5, df1 = 1, df2 = 27, p-value < 2.2e-16

> resettest(y2 ~ x , power=2, type="regressor")

	RESET test

data:  y2 ~ x 
RESET = 2.976, df1 = 1, df2 = 27, p-value = 0.09594

Ну, а в нашем случае, это выглядит так

> mdl<-lm(Y~A+B+C,data=MyD)
> resettest(mdl)

	RESET test

data:  mdl 
RESET = 204.7904, df1 = 2, df2 = 8, p-value = 1.347e-07

> mdl<-lm(Y~A+B+C+A:B+B:C+A:C+A:B:C+I(A^2)+I(B^2)+I(C^2),data=MyD)
> resettest(mdl)

	RESET test

data:  mdl 
RESET = 0.925, df1 = 2, df2 = 1, p-value = 0.5923

[Hogfather]

Допустим, стоит задача обрабатывать большие массивы данных. Работа с электронными таблицами, конечно, хороша, но иногда размер становится такой, что фильтрация и обработка пересечений и объединений таблиц становится нереальной.

Вот тут есть статья
Перенос табличных данных в SQLite с дальнейшим использованием в GNU R

Всем хороша, но содержит некоторые неточности.

Во-первых, пакет называется "RSQLite" (в принципе, там об этом сказано, но в примере кода стоит "require(SQLite)", что не работает)

Пример работающего кода

> # Подключаем библиотеку
> require(RSQLite) 
Загрузка требуемого пакета: RSQLite
Загрузка требуемого пакета: DBI
> # Подключаемся к БД
> drv <- dbDriver("SQLite")
> setwd("C:/Users/Hogfather/R project")
> fname<-"test.sqlite"
> con <- dbConnect(drv, dbname = fname)
> # Выполняем запрос. Все данные оказались в переменной my.data
> my.data<-dbGetQuery(con,"select * from MyTable where ID=666")
> # Отключаемся от БД
> dbDisconnect(con)
[1] TRUE
> dbUnloadDriver(drv)
[1] TRUE

Во вторых, на мой взгляд, гораздо удобнее импортировать в SQLite csv файлы с помощью дополнения SQLite Manager для Mozilla FireFox. Получается быстрее и сразу виден результат.

[will]

как сделать так, чтобы команда install.packages() загружала файлы не в директорию по умолчанию (на диске С), а в другую, на диске D? Сама R установлена на диске D
setwd() не меняет ситуации

[Hogfather]

will, спасибо за интересный вопрос. Ответ к нему легко обраружить в Help

?install.packages

Имеем

Install Packages from Repositories or Local Files

Description

Download and install packages from CRAN-like repositories or from local files.

Usage

install.packages(pkgs, lib, repos = getOption("repos"),
contriburl = contrib.url(repos, type),
method, available = NULL, destdir = NULL,
dependencies = NA, type = getOption("pkgType"),
configure.args = getOption("configure.args"),
configure.vars = getOption("configure.vars"),
clean = FALSE, Ncpus = getOption("Ncpus", 1L),
libs_only = FALSE, INSTALL_opts, ...)
Arguments

pkgs
character vector of the names of packages whose current versions should be downloaded from the repositories.

If repos = NULL, a character vector of file paths of ‘.tar.gz’ files. These can be source archives or binary package archive files (as created by R CMD build --binary). On a CRAN build of R for Mac OS X these can be ‘.tgz’ files containing binary package archives. Tilde-expansion will be done on the file paths.

If this is missing or a zero-length character vector, a listbox of available packages is presented where possible in an interactive R session.

lib
character vector giving the library directories where to install the packages. Recycled as needed. If missing, defaults to the first element of .libPaths().

Попробуйте

install.packages(pkgs, lib="d:\\Разная фигня\\R")

или исправьте .libpath

> .libPaths()
[1] "/Library/Frameworks/R.framework/Versions/2.15/Resources/library"
[2] "/Applications/RStudio.app/Contents/Resources/R/library"

The library search path is initialized at startup from the environment variable R_LIBS (which should be a colon-separated list of directories at which R library trees are rooted) followed by those in environment variable R_LIBS_USER. Only directories which exist at the time will be included.

By default R_LIBS is unset, and R_LIBS_USER is set to directory ‘R/R.version$platform-library/x.y’ of the home directory (or ‘Library/R/x.y/library’ for Mac OS X AQUA builds), for R x.y.z.

[will]

а деинсталлировать ранее установленные пакеты (дополнительные, не основные) средствами R можно?
И про обновление версий R. Только полная деинсталляция, и установка по-новой? Или как-то это можно упростить?

install.packages(pkgs, lib="d:\\Разная фигня\\R")

не срабатывает, устанавливается на С
хотя .libPaths() дает правильный путь

[Hogfather]

Сообщение от will

а деинсталлировать ранее установленные пакеты (дополнительные, не основные) средствами R можно?

> remove.packages(c("qcc","MASS"))

Сообщение от will

И про обновление версий R. Только полная деинсталляция, и установка по-новой? Или как-то это можно упростить?

Ну, смотря какая версия: мажорная или минорная. Для минорных можно просто обновлять, для мажорных -- с танцами и бубном.

P.S. На большинство вопросов есть ответ вот тут http://cran.rstudio.com/bin/windows/...o-upgrade_003f
Достаточно быстро проходит операция и все пакеты переносятся. Проверял.

[Hogfather]

R Cookbook (O'Reilly Cookbooks)


Книга из серии Must Have. Постоянно пользуюсь в качестве настольного справочника
Оценка: 5 из 5

Data Analysis and Graphics Using R: An Example-Based Approach (Cambridge Series in Statistical and Probabilistic Mathematics)



Весьма толковая книга. Если я когда-нибудь надумаю писать учебник по статистике, то я использую именно такой способ изложения. Рассматриваются основные подходы к статистическому анализу, много примеров кода, но книжка более серьёзная чем предыдущая, Пользуюсь достаточно часто,

Оценка: 5 из 5

Data Mining with R: Learning with Case Studies (Chapman & Hall/CRC Data Mining and Knowledge Discovery Series)



Книга разочаровала. Много воды и ответ на вопросы, ради которой её покупал не нашёл. Сама идея книги весьма неплоха (взять некий набор данных и показать на них Data Mining), но реализация немного подкачала.

Оценка: 3 из 5

[Hogfather]

R Graphics Cookbook [Paperback]



Толковая книжка. Открыл для себя много нового.

Оценка: 5 из 5


Doing Bayesian Data Analysis: A Tutorial with R and BUGS [Hardcover]



Зело интересно написанная книжка по байесовской статистике. С примерами, моделями, объясняется зачем это нужно и что с этим делать.

Оценка: 5 из 5

[Hogfather]

Задал мне тут вопрос Дмитрий В., говорит, Хогфазер, а есть ли в R лица Чернова, зело они мне на конференции понравились, хочу научному руководителю рожу гнусную состроить, но на научной основе?

Полез посмотреть, есть.

Без комментария, код по трем словарям

> library(aplpack)
> faces(rbind(c(2.95,3.4,1.25,9.1),c(2.64,2.34,1.3,8.5),c(2.39,1.51,1.20,7.10)),
+       labels=c("Словарь 1","Словарь 2","Словарь 3"),face.type=1)
effect of variables:
 modified item       Var   
 "height of face   " "Var1"
 "width of face    " "Var2"
 "structure of face" "Var3"
 "height of mouth  " "Var4"
 "width of mouth   " "Var1"
 "smiling          " "Var2"
 "height of eyes   " "Var3"
 "width of eyes    " "Var4"
 "height of hair   " "Var1"
 "width of hair   "  "Var2"
 "style of hair   "  "Var3"
 "height of nose  "  "Var4"
 "width of nose   "  "Var1"
 "width of ear    "  "Var2"
 "height of ear   "  "Var3"

Под спойлером -- более наглядный пример с моими данными.

Более наглядный пример

Добавлено через 13 часов 52 минуты
Как выяснилось, морды можно рисовать также с помощью библиотеки TeachingDemos. Что характерно, там две функции: face и face2. Первая практически совпадает с вышеописанной, но более убого по настройкам. А образец вывода второй функции прилагаю.

> library(TeachingDemos)
> faces2(rbind(c(2.95,3.4,1.25,9.1),c(2.64,2.34,1.3,8.5),c(2.39,1.51,1.20,7.10)),labels=c("Словарь 1","Словарь 2","Словарь 3"))

[Hogfather]

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

Итак, что нужно для счастья.
1. RStudio -- крайне желательна.
2. Установленный пакет knitr в R

В меню RStudio вызываем File->New->R Markdown

В принципе, образец виден сразу. Но мы сделаем свой пример:

Пример (см. MarkdownDemo.Rmd)

Обзор возможностей Markdown в RStudio
========================================================

Введение
--------

Одной из интересных воможностей исследования с использованием R является создание самодокументирующихся скриптов с использованием пакета **Knitr** и языка разметки [markdown](http://daringfireball.net/projects/markdown/syntax).

Рассмотрим простую задачу построения регрессионной модели для данных из пакета  **MASS** по весам кошачьих сердец

```{r}
library(MASS)
head(cats)
```

Как мы видим, имеются три колонки, с числом строк `r nrow(cats)` и числом столбцов - `r ncol(cats)` (обратите внимание, что цифры взяты прямо из данных). 

Разведочный анализ данных
-------------------------

Проведем экспресс-анализ данных, построив графики
```{r warning=FALSE}
boxplot(Hwt~Sex,data=cats,main="Распределение веса кошачьих сердец\nв зависимости от пола")

```
Можно еще построить диаграмму рассеивания. Обратите внимание, чтобы русский язык отображался корректно следует использовать кодировку UTF-8. Также зададим побольше размер для картинки.

```{r pic1,fig.width=10,fig.height=10, warning=FALSE}
plot(Hwt~Bwt,data=cats,col=Sex,pch=19,main="Вес кошачьих сердец в зависимости от массы тушки")

```
Русские заголовки имеют привычку создавать разные *warnings*, которые мы подавили добавив выше сответсвующую опцию.

Построение модели
-----------------

Дальше можно попробовать построить регрессионуую модель

```{r}
summary(mdl1<-lm(Hwt~Bwt+Sex,data=cats))
```
С помощью дополнительного пакета **xtabs** можно вывести красиво таблицу.

```{r results='asis'}
library(xtable)
print(xtable(mdl1),type = "html")
```

Вернемся к нашему примеру. В данном случае пол скотинки нам совершенно не значим, можно его убрать

```{r}
summary(mdl2<-lm(Hwt~Bwt,data=cats))
```
Ну, теперь самое время построить диагностические графики

```{r fig.width=10,fig.height=10}
par(mfrow=c(2,2))
plot(mdl2)
par(mfrow=c(1,1))
```

Вот так с шутками и прибаутками была простоена регрессия.

$$ y=`r mdl2$coeff[2] ` x + `r mdl2$coeff[1]`+\epsilon $$

(это, типа, формула)

Заключение
------------
По-моему, очень удобный инструмент.

Дальше, жмем кнопку "Knit HTML" и получаем файл html, картинки, а также markup файл, в котором все команды R уже выполнены.

Вот тут нам следует обратить внимание на другую чудесную программу с названием pandoc. Установив оную и запустив из командной строки в нужном каталоге можно получить на выходе Word документ.
Запуск выглядит так

pandoc -o MarkdownDemo.docx MarkdownDemo.md

Как видно в приведенном во вложении примере, пакет xtable в данном случае бесполезен, поскольку генерит Html. Понятно, что правильнее и интереснее все делать через LaTeX, но, не разводя лишнего холивара, смею заметить, что не все еще им прониклись. Поэтому путь будет Markdown.