Портал аспирантов - GNU R: Вопросы и ответы

Портал аспирантов (http://www.aspirantura.spb.ru/forum/index.php)

- Software (программное обеспечение) (http://www.aspirantura.spb.ru/forum/forumdisplay.php?f=107)

- - GNU R: Вопросы и ответы (http://www.aspirantura.spb.ru/forum/showthread.php?t=10501)

gnu r регрессионный анализ. Ответы на вопросы.

Как выяснилось в переписке, вопрос касался также 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

Работа с SQLite в в GNU R

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

Вот тут есть статья
Перенос табличных данных в 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. Получается быстрее и сразу виден результат.

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

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.

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

Цитата:

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

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

Цитата:

Сообщение от will (Сообщение 315032)

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

Код:

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

Цитата:

Сообщение от will (Сообщение 315032)

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

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

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

GNU R - краткий обзор литературы на английском языке

R Cookbook (O'Reilly Cookbooks)
http://ecx.images-amazon.com/images/...SH20_OU01_.jpg

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

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

http://ecx.images-amazon.com/images/...SH20_OU01_.jpg

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

Оценка: 5 из 5

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

http://ecx.images-amazon.com/images/...SH20_OU01_.jpg

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

Оценка: 3 из 5

R Graphics Cookbook [Paperback]

http://ecx.images-amazon.com/images/...SH20_OU01_.jpg

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

Оценка: 5 из 5

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

http://ecx.images-amazon.com/images/...SH20_OU01_.jpg

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

Оценка: 5 из 5

GNU R - лица Чернова

Задал мне тут вопрос Дмитрий В., говорит, Хогфазер, а есть ли в 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"

http://aspirantura.spb.ru/forum/pict...pictureid=1092

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

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

http://aspirantura.spb.ru/forum/pict...pictureid=1093

Добавлено через 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"))

http://aspirantura.spb.ru/forum/pict...pictureid=1094

R, Markdown, pandoc и другие

В процессе обучения курсам мастерства работы с 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.