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

[Hogfather]

3^http://aspirantura.spb.ru/forum/picture.php?albumid=36&pictureid=979 | Внимание![BR]Автор темы совместно с Администрацией портала просит писать в эту тему только относящееся к GNU R. Благодарности, разговоры о погоде будут безжалостно удаляться. Все "чмоки", пожалуйста, во флейме. Вопросы по R прошу задавать в личных сообщениях, чтобы тема была удобна для восприятия. Надеюсь на понимание.


Итак, поскольку я показал как пользоваться R для решения задачи подбора распределения (http://aspirantura.spb.ru/forum/showpost.php?p=289709&postcount=46), а добро не должно остаться безнаказанным, в личные сообщения посыпались вопросы. Чтобы два раза не вставать, ответы публикую сюда.

Рассмотрено в других темах (пример кода)

1. Подгонка распределения. fitdistr из MASS (http://aspirantura.spb.ru/forum/showpost.php?p=289709&postcount=46)
2. Подгонка распределения. fitdist из fitdistrplus (http://aspirantura.spb.ru/forum/showpost.php?p=290182&postcount=56). Лучше!
3. Нелинейный метод наименьших квадратов (http://aspirantura.spb.ru/forum/showpost.php?p=290073&postcount=50)
4. Гамма-распределение. Вывод графика. (http://aspirantura.spb.ru/forum/showthread.php?t=10388)
5. Диаграмма "Ящик с усами" для нескольких векторов (http://aspirantura.spb.ru/forum/showpost.php?p=290226&postcount=58)
6. Упрощение себе жизни с помощью написания функции в R (http://aspirantura.spb.ru/forum/showpost.php?p=290374&postcount=60)
7. Подгонка МНК (http://aspirantura.spb.ru/forum/showpost.php?p=290854&postcount=5)

Q: Где взять R?

A: Лучше всего в специально предназначенном для этого месте, называется CRAN. Русское зеркало находится вот тут (http://cran.gis-lab.info/). Выбираете операционную систему, далее скачиваете и устанавливаете.

Q:C чего начать?

A: Почитать статьи вот тут (http://herba.msu.ru/shipunov/software/r/r-ru.htm)

Q: Как попробовать то, что показано в примерах?

A: Очень просто. Запускаем R, копируем команду из примера в консоль интерпретатора. Команды начинаются с ">", сам символ ">" приглашения интерпретатора R не копируем.
Например, скопируйте все строки из примера ниже (я оставил только команды из вышеприведенной заметки, а символ ">" стёр) и запустите в R (Не забудьте в конце нажать Enter, чтобы выполнилась последняя команда).

# Тестовый пример. Скопировать все строки и вставить в окно интерпретатора R
LT<-c(rep(1,9),rep(2,267),rep(3,2843),rep(4,5450),rep( 5,6564),rep(6,7044),rep(7,7518),rep(8,7071),rep(9, 5620),rep(10,4016),rep(11,2545),rep(12,1494),rep(1 3,854),rep(14,416),rep(15,214),rep(16,122),rep(17, 53),rep(18,16),rep(19,7),rep(20,2),21,22)
summary(LT)
length(LT)
# строим график
old.par <- par(mfrow=c(2,2))
hist(LT,main="Распределение букв",ylab="Число наблюдений",xlab="Число букв")
hist(LT,freq=F,ylab="Вероятность",xlab="Число букв",main="Распределение букв")
plot(ecdf(LT),verticals=T,main="График функции распределения")
boxplot(LT,main="Диаграмма Ящик-с-Усами",xlab="Число букв",horizontal=T)
par(old.par)


Вот что получится в результате:
http://aspirantura.spb.ru/forum/picture.php?albumid=36&pictureid=973

Q: Замечательно, получилось. А как теперь скопировать картинку?

A: В MS Windows -- правая кнопка мыши на рисунке и выбрать нужный Вам вариант. Более продвинутые способы рассмотрим позже.

Q: А как бы мне открыто новую картинку на экране, да не поверх предыдущей?

A: Воспользуйтесь командой windows()
>windows()
Можно задать размеры окна в дюймах и еще кучу всего. Подробнее в справке.
>?windows


Q: У меня не подключилась библиотека MASS.

A: Понимаю, нужно её установить. В меню "Пакеты" выберите "Установить пакеты", дальше разберетесь сами.

Q: Я ввел неверную команду или хочу повторить предыдущую, опять заново набивать надо?

A: Нет. Достаточно нажать стрелку вверх. Это позволит "путешествовать" по истории команд, а если проскочили лишнее, то стрелка вниз. Далее можно редактировать строку. Потом нажмите Enter и она выполнится. Либо наберите команду history(), в отдельном окне покажется список последних введенных команд (25 штук). Если нужно больше, то поставьте соответствующую цифру в скобках, например history(100). Из открывшегося окна можно копировать нужные команды.

[Hogfather]

Q; Как корректно вывести графики с русскими буквами в PDF (SVG)?

A: Путем долгих экспериментов я пришел к выводу, что для меня лучше всего использовать библиотеку cairoDevice, которая использует библиотеку GTK+ .

Пример (использованы данные из темы про параметрические методы в филологии)


> library(cairoDevice)
> ##### Формируем PDF
> Cairo_pdf("test.pdf", width = 7, height = 7, pointsize = 10)
> boxplot(LT,main="Диаграмма ящик-с усами")
> dev.off()
windows
2
> ##### Формируем SVG
> Cairo_svg("test.svg", width = 7, height = 7, pointsize = 10)
> boxplot(LT,main="Диаграмма ящик-с усами")
> dev.off()
windows
2

> ##### Формируем PNG
> Cairo_png("test.png", width = 7, height = 7, pointsize = 10)
> boxplot(LT,main="Диаграмма ящик-с усами")
> dev.off()
windows
2


Результаты во вложении. Недавно вышла новая библиотека Cairo 2012 года, проверил, тоже работает.

Q: Можно ли редактировать результирующий график?
A: Да. Достаточно вывести его в SVG и воспользоваться графическим редактором, который понимает это формат. Например, Inkscape (http://sourceforge.net/projects/inkscape/).

[Hogfather]

Q: Как удобнее редактировать скрипт в R.
A: Меню Файл->Новый скрипт. В новом окне пишите команды подряд.
Любую команду можно запустить в окне интерпретатора, нажав Ctrl-R.
Если нужно запустить весь скрипт (отдельные строки), то выделите нужные строки и нажмите Ctrl-R. Выделение всего скрипта -- Ctrl-A.

Q: Как отредактировать данные.
A: Например, воспользоваться командой edit или fix. Команда edit не меняет исходные данные, а создает копию. Поэтому её следует использовать так
>LT1<-edit(LT)
Если нужно внести правки, без сохранения исходных данных, то тогда всё проще.
>fix(LT)

Q: Как загрузить данные в R.
Допустим, у вас имеются данные, которые нужно загрузить в R.
Существует несколько способов это сделать. Самые полезные.
1) Загрузить из csv файла
>Data2012<-read.csv("D://R//data2012.csv", header=TRUE, sep = ";", quote="\"", dec=".")
По команде читается файл по адресу D:/R/data2012.csv (Обратите внимание на двойные слеш-символы), первая строка файла содержит заголовки (английский, без пробелов), разделитель записей символ ";", символ кавычки ["], десятичная точка -- ".".

2) Получить из Clipboard
> xx<-read.table("clipboard",header=T)

http://aspirantura.spb.ru/forum/picture.php?albumid=36&pictureid=978


3) Получить прямо из Интернет
Для этого нужна библиотека XML. Вот таким нехитрым способом скачиваем 5 таблицу из статьи в Википедии.

> library(XML)
> url<-"http://en.wikipedia.org/wiki/World_population"
> (tbls<-readHTMLTable(url,which=5))
Rank Country / Territory Population Date
1 1 В*China[note 2] 1,354,210,000 November 12, 2012
2 2 В*India 1,210,193,422 March 2011
3 3 В*United States 314,750,000 November 12, 2012
4 4 В*Indonesia 238,400,000 May 2010
5 5 В*Brazil 197,268,000 November 12, 2012
6 6 В*Pakistan 181,201,000 November 12, 2012
7 7 В*Nigeria 170,123,740 July 2012
8 8 В*Bangladesh 161,083,804 July 2012
9 9 В*Russia 141,927,297 January 1, 2010
10 10 В*Japan 127,610,000 May 1, 2012
В*% of world\npopulation Source
1 19.2% [73]
2 17% [74]
3 4.46% [75]
4 3.32% [76]
5 2.8% [77]
6 2.57% [78]
7 2.41% [79]
8 2.28% [80]
9 2.013% [81]
10 1.81% [82]
>



4) Ввести вручную

> mydata <- data.frame(age=numeric(0), gender=character(0), weight=numeric(0))
> mydata <- edit(mydata)


5) Получить из базы данных
Подробнее см. тут (http://www.statmethods.net/input/dbinterface.html)

[SergeyMastitsky]

Коллеги, приглашаю посетить мой блог "R: Анализ и визуализация данных": r-analytics.blogspot.com
Надеюсь, он окажется кому-то полезным в поисках ответов на вопросы, касающиеся R.

[Hogfather]

приглашаю посетить мой блог "R: Анализ и визуализация данных": r-analytics.blogspot.com
Подтверждаю. Блог (http://r-analytics.blogspot.ru) полезный, кошерный, рекомендуется к посещению и вдумчивому прочтению.

[will]

Вопрос
Есть ли в R
1)Фурье-анализ
2)подбор нелинейного уравнения (из фиксированного перечня+пользовательск.) для данного набора 2-х или 3-х мерных данных , а потом отобразить графически исходные данные и аппроксимирующую кривую \поверхность + повращать (для трехмерной графики)?

[Hogfather]

Есть ли в R
1)Фурье-анализ
Конечно, входит. Для быстрого преобразования Фурье даже пакетов подключать дополнительно не надо (пакет stats)

fft {stats} R Documentation
Fast Discrete Fourier Transform

Description

Performs the Fast Fourier Transform of an array.

Usage

fft(z, inverse = FALSE)
mvfft(z, inverse = FALSE)
Arguments

z
a real or complex array containing the values to be transformed.

inverse
if TRUE, the unnormalized inverse transform is computed (the inverse has a + in the exponent of e, but here, we do not divide by 1/length(x)).

Value

When z is a vector, the value computed and returned by fft is the unnormalized univariate Fourier transform of the sequence of values in z. When z contains an array, fft computes and returns the multivariate (spatial) transform. If inverse is TRUE, the (unnormalized) inverse Fourier transform is returned, i.e., if y <- fft(z), then z is fft(y, inverse = TRUE) / length(y).

By contrast, mvfft takes a real or complex matrix as argument, and returns a similar shaped matrix, but with each column replaced by its discrete Fourier transform. This is useful for analyzing vector-valued series.

The FFT is fastest when the length of the series being transformed is highly composite (i.e., has many factors). If this is not the case, the transform may take a long time to compute and will use a large amount of memory.

Source

Uses C translation of Fortran code in Singleton (1979).

References

Becker, R. A., Chambers, J. M. and Wilks, A. R. (1988) The New S Language. Wadsworth & Brooks/Cole.

Singleton, R. C. (1979) Mixed Radix Fast Fourier Transforms, in Programs for Digital Signal Processing, IEEE Digital Signal Processing Committee eds. IEEE Press.


> x <- 1:4
> fft(x)
[1] 10+0i -2+2i -2+0i -2-2i
> fft(fft(x), inverse = TRUE)/length(x)
[1] 1+0i 2+0i 3+0i 4+0i

2)подбор нелинейного уравнения (из фиксированного перечня+пользовательск.) для данного набора 2-х или 3-х мерных данных , а потом отобразить графически исходные данные и аппроксимирующую кривую \поверхность + повращать (для трехмерной графики)?
Естественно есть, причем разными методами. Я уже рассматривал (http://aspirantura.spb.ru/forum/showpost.php?p=290073&postcount=50) нелинейный МНК. Трехмерные графики тоже есть (http://www.personal.psu.edu/ljk20/3dinR.pdf). Вот тут (http://stackoverflow.com/questions/6720526/plotting-3d-data-in-r) есть разбор пакетов.

Классический эффектный пример нормального распределения в двух измерениях [1] (http://www.ejwagenmakers.com/misc/Plotting_3d_in_R.pdf).


mu1<-0 # setting the expected value of x1
mu2<-0 # setting the expected value of x2
s11<-10 # setting the variance of x1
s12<-15 # setting the covariance between x1 and x2
s22<-10 # setting the variance of x2
rho<-0.5 # setting the correlation coefficient between x1 and x2
x1<-seq(-10,10,length=41) # generating the vector series x1
x2<-x1 # copying x1 to x2
#
f<-function(x1,x2)
{
term1<-1/(2*pi*sqrt(s11*s22*(1-rho^2)))
term2<--1/(2*(1-rho^2))
term3<-(x1-mu1)^2/s11
term4<-(x2-mu2)^2/s22
term5<--2*rho*((x1-mu1)*(x2-mu2))/(sqrt(s11)*sqrt(s22))
term1*exp(term2*(term3+term4-term5))
} # setting up the function of the multivariate normal density
#
z<-outer(x1,x2,f) # calculating the density values
#
persp(x1, x2, z,
main="Two dimensional Normal Distribution",
sub=expression(italic(f)~(bold(x))==frac(1,2~pi~sq rt(sigma[11]~
sigma[22]~(1-rho^2)))~phantom(0)^bold(.)~exp~bgroup("{",
list(-frac(1,2(1-rho^2)),
bgroup("[", frac((x[1]~-~mu[1])^2, sigma[11])~-~2~rho~frac(x[1]~-~mu[1],
sqrt(sigma[11]))~ frac(x[2]~-~mu[2],sqrt(sigma[22]))~+~
frac((x[2]~-~mu[2])^2, sigma[22]),"]")),"}")),
col="lightgreen",
theta=30, phi=20,
r=50,
d=0.1,
expand=0.5,
ltheta=90, lphi=180,
shade=0.75,
ticktype="detailed",
nticks=5) # produces the 3-D plot
#
mtext(expression(list(mu[1]==0,mu[2]==0,sigma[11]==10,sigma[22]==10,sigma[12
]==15,rho==0.5)), side=3) # adding a text line to the graph

http://aspirantura.spb.ru/forum/picture.php?albumid=36&pictureid=991

[Hogfather]

Рассмотренная ранее диаграмма Каллена и Фрея (Cullen and Frey graph) всем хороша, но она, во-первых, на языке эвентуального противника, во-вторых, нет возможности менять заголовок. Такое безобразие долго продолжаться не могло и, в процессе работы над статьёй, пришлось подправить функцию и русифицировать её вывод. Для этого идём на CRAN (http://cran.r-project.org/web/packages/fitdistrplus/index.html), скачиваем "Package source", распаковываем, в папке R находим нужную функцию и редактируем.

Результат во вложении. Пользуйтесь на здоровье.

Если работаете в стандартной оболочке R, то достаточно распаковать архив, в меню "Файл" выбрать "Загрузить код R..." и дальше спокойно пользоваться оной функцией, например так, как показано ниже. Обратите внимание, параметр "main" меняет заголовок диаграммы. Этого нет в стандартном функционале descdistr.


> descdist(LT,boot=100,main="Диграмма Каллена и Фрея\n(Редакция Hogfather'а для Портала Аспирантов)")
summary statistics
------
min: 1 max: 22
median: 7
mean: 7.199839
estimated sd: 2.628829
estimated skewness: 0.519882
estimated kurtosis: 3.143716
>


http://aspirantura.spb.ru/forum/picture.php?albumid=36&pictureid=992

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

[Hogfather]

Прочитав про RStudio (http://r-analytics.blogspot.ru/p/rstudio.html) решил её обжить и делюсь решением первых проблем.

Проблемы с кодировкой
Кодировка внутри файлов воспринимается весьма неплохо, если скрипт прочитался неверно, его можно прочитать в другой кодировке с помощью команды Файл->Reopen with Encoding...

Но! Если вы на компьютере имеете логин на русском языке, и в папке Users (Пользователи) у Вас папка на русском языке, соответствующая логину, то в Win7 возникают проблемы. Главная, RStudio категорически не видит установленные пакеты, поскольку они устанавливаются в C:/Users/Логин/Documents/R/... Вместо вот этого "Логин", который у Вас вполне может быть "Милый суслик", программа пойдет в папку C:/Users/Логин/Documents/R/ (Логин = попытка прочитать в кодировке Windows-1251 UTF-8 слово "Логин")
Поскольку ни считать её, ни создать она не может, у программы возникает печалька, и она не видит установленных пакетов, а также не может поставить новые. Вот такая ерунда.

Выход есть! Нужно сделать ссылку типа "связь каталогов" с таким идиотским именем в папке Users
Для этого (как вариант).
1. Скачиваем и устанавливаем Far (http://www.farmanager.com/). Хорошему грустному человеку он завсегда пригодится.
2. Заходим в каталог C:/Users/
3. Нажимаем Alt-F6 и в меню выбираем тип ссылки "Связь каталогов" или "directory junction", зависит от языка в Far.
4. На противоположной панели появляется ссылка.
5. Переименовываем её в кракозябровое имя, которое мы подглядим в ошибке RStudio, выделим, скопируем, нажмем в Far F6 и вставим его в нужное место окна, дополнив строку пути нашими кракозябрами
6. Перезапускаем RStudio и радуемся жизни.

[Hogfather]

Наши постоянные читатели спрашивают

Будет замечательно, если Вы рассмотрите данную задачу как пример.

Задача: получить зависимость z(x,y)
Исходные данные:

x <- c(75, 75, 75, 75, 100, 125, 125, 150, 150, 175, 175, 200, 200, 225, 250, 250, 250, 275, 300, 300, 300, 300)
y <- c(375, 625, 1000, 1500, 500, 375, 1125, 750, 1250, 825, 1500, 375, 1000, 1125, 750, 1250, 1500, 1375, 375, 825, 1500, 1250)
z <- c(1.452, 3.705, 18.857, 50.434, 1.956, 2.272, 40.306, 3.368, 47.686, 4.464, 55.63, 2.611, 9.077, 15.304, 7.559, 27.257, 60.174, 35.965, 3.355, 22.491, 38.036, 54.237)


Ну, что ж. Начнем. Загоняем данные в R и делаем Data Frame (таблицу)

> x <- c(75, 75, 75, 75, 100, 125, 125, 150, 150, 175, 175, 200, 200, 225, 250, 250, 250, 275, 300, 300, 300, 300)
> y <- c(375, 625, 1000, 1500, 500, 375, 1125, 750, 1250, 825, 1500, 375, 1000, 1125, 750, 1250, 1500, 1375, 375, 825, 1500, 1250)
> z <- c(1.452, 3.705, 18.857, 50.434, 1.956, 2.272, 40.306, 3.368, 47.686, 4.464, 55.63, 2.611, 9.077, 15.304, 7.559, 27.257, 60.174, 35.965, 3.355, 22.491, 38.036, 54.237)
> myData<-data.frame(x,y,z)

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

> library(lattice)
> splom(~myData[,c("x","y","z")],cex.labels=1.2,varnames=c("x","y","z"))

http://aspirantura.spb.ru/forum/picture.php?albumid=36&pictureid=999

Кого-как, а меня не впечатлило. Икс у нас явно не при делах, а Игрек на линейную модель никак не тянет. Напрашивается экспонента. Проверим эту нехитрую мысль.
Введем фиктивную переменную lnZ равную натуральному логарифму от Z. В R это делается так.

> myData$lnZ<-log(myData$z)
> splom(~myData[,c("x","y","lnZ")],cex.labels=1.2,varnames=c("x","y","ln(z)"),col="red")

http://aspirantura.spb.ru/forum/picture.php?albumid=36&pictureid=1001
Поясняю, мы от модель вида http://www.texify.com/img/%5Cnormalsize%5C%21z%3D%5Cbeta_%7B1.1%7D%20x%2B%5C beta_%7B1.2%7D%20y%20%2B%5Cbeta_0%20%2B%20%5Cepsil on.gif плавно перешли к http://www.texify.com/img/%5Cnormalsize%5C%21%5Cnormalsize%5C%21z%3De%5E%7B% 28%5Cbeta_%7B1.1%7D%20x%2B%5Cbeta_%7B1.2%7D%20y%20 %2B%5Cbeta_0%20%2B%20%5Cepsilon%29%7D.gif. В данном случае эпсилон у нас изображает остатки модели, которые, в идеальном варианте, должны иметь нормальное распределение.

Третий вариант мы можем получить логарифмированием всего и вся. В таком случае имеем http://www.texify.com/img/%5Cnormalsize%5C%21z%3DAx%5E%7B%5Cbeta_%7B1.1%7D%7 D%20y%5E%7B%5Cbeta_%7B1.2%7D%7D%5Ctimes%5Cmathscr% 7BE%7D.gif

> splom(~log(myData[,c("x","y","z")]),cex.labels=1.2,varnames=c("ln(x)","ln(y)","ln(z)"),col="dark green")

http://aspirantura.spb.ru/forum/picture.php?albumid=36&pictureid=1002

После такой работы уже можно и к моделированию переходить. Чисто чтоб поржать сделаем все три модели.

> summary(lm1<-lm(z~x+y))

Call:
lm(formula = z ~ x + y)

Residuals:
Min 1Q Median 3Q Max
-15.717 -8.306 1.616 6.760 17.840

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -21.233747 6.918102 -3.069 0.00631 **
x 0.001910 0.028261 0.068 0.94683
y 0.045646 0.005698 8.010 1.64e-07 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 10.24 on 19 degrees of freedom
Multiple R-squared: 0.7846, Adjusted R-squared: 0.762
F-statistic: 34.61 on 2 and 19 DF, p-value: 4.624e-07

> summary(lm2<-lm(log(z)~x+y))

Call:
lm(formula = log(z) ~ x + y)

Residuals:
Min 1Q Median 3Q Max
-0.61700 -0.30351 -0.03218 0.19850 0.81873

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.5939071 0.3067107 -1.936 0.0678 .
x 0.0016541 0.0012529 1.320 0.2025
y 0.0029021 0.0002526 11.488 5.4e-10 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.4538 on 19 degrees of freedom
Multiple R-squared: 0.889, Adjusted R-squared: 0.8773
F-statistic: 76.08 on 2 and 19 DF, p-value: 8.53e-10

> summary(lm3<-lm(log(z)~log(x)+log(y)))

Call:
lm(formula = log(z) ~ log(x) + log(y))

Residuals:
Min 1Q Median 3Q Max
-0.9323 -0.2865 0.1118 0.3488 0.6532

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -14.5313 1.6847 -8.626 5.38e-08 ***
log(x) 0.2550 0.2275 1.121 0.276
log(y) 2.3263 0.2280 10.203 3.81e-09 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.5058 on 19 degrees of freedom
Multiple R-squared: 0.8621, Adjusted R-squared: 0.8476
F-statistic: 59.41 on 2 and 19 DF, p-value: 6.682e-09

Смотрим колонку Pr(>|t|) и количество звездочек сбоку. Начинаем модифицировать модель "a posteriori". Собственно говоря, Икс нам совершенно не нужен во всех его проявлениях. Непонятно зачем он вообще тут оказался. Возьмем вторую модель, она посимпатичнее и выкинем лишнее.
> summary(lm4<-lm(log(z)~y))

Call:
lm(formula = log(z) ~ y)

Residuals:
Min 1Q Median 3Q Max
-0.66190 -0.35391 -0.08276 0.28610 1.01267

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.3660607 0.2582089 -1.418 0.172
y 0.0029897 0.0002483 12.043 1.28e-10 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.4622 on 20 degrees of freedom
Multiple R-squared: 0.8788, Adjusted R-squared: 0.8727
F-statistic: 145 on 1 and 20 DF, p-value: 1.277e-10

Смотрим дальше. Свободный член (не радуйтесь девушки, это не то, что вы подумали) "Intercept" нам тоже не нужен. Да и пёс с ним.
> summary(lm5<-lm(log(z)~-1+y))

Call:
lm(formula = log(z) ~ -1 + y)

Residuals:
Min 1Q Median 3Q Max
-0.78397 -0.35742 -0.05766 0.18365 0.91500

Coefficients:
Estimate Std. Error t value Pr(>|t|)
y 2.664e-03 9.699e-05 27.47 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.4732 on 21 degrees of freedom
Multiple R-squared: 0.9729, Adjusted R-squared: 0.9716
F-statistic: 754.6 on 1 and 21 DF, p-value: < 2.2e-16

Обратите внимание какой у нас славный получился Эр квадрат. Дай бог всякому.
Моделька имеет вид http://www.texify.com/img/%5Cnormalsize%5C%21%5Cnormalsize%5C%21z%3De%5E%7B2 .664%5Ctimes%2010%5E%7B-3%7D%20y%7D%20%20%20%5Ctimes%20%5Cmathscr%7BE%7D.g if
Едем дальше. Впереди дисперсионный анализ.
> anova(lm5)
Analysis of Variance Table

Response: log(z)
Df Sum Sq Mean Sq F value Pr(>F)
y 1 168.946 168.946 754.61 < 2.2e-16 ***
Residuals 21 4.702 0.224
И графическое отображение остатков модели.
> par(mfrow=c(2,2))
> plot(lm5)
> par(mfrow=c(1,1))
Ну, лично мне нравится.
http://aspirantura.spb.ru/forum/picture.php?albumid=111&pictureid=1003

Итак, у нас есть модель и остатки. Если все хорошо, они должны иметь нормальное распределение с матожиданием 0. Судя по картинке (график квантилей) мы должны получить неплохой результат в тесте не нормальность распределения остатков модели.
Делаем тест Шапиро-Уилка. Смотрим и балдеем.
> shapiro.test(lm5$residuals)

Shapiro-Wilk normality test

data: lm5$residuals
W = 0.9628, p-value = 0.5488

Раз все так здорово, находим стандартное отклонение.

> sd(lm5$residuals)
[1] 0.4700086
> summary(lm5$residuals)
Min. 1st Qu. Median Mean 3rd Qu. Max.
-0.78400 -0.35740 -0.05766 -0.05331 0.18360 0.91500


Sapienti sat.

Понятно, что в форумной заметке нельзя рассказать весь курс, но основные идеи я, надеюсь, рассказал.

P.S. Вот тут спрашивают, а как нарисовать 3D-график этого безобразия. Отвечаю. Просто. Например, так.
> library(scatterplot3d)
> scatterplot3d(x , y, z,highlight.3d=T,main="Вот такой график")

http://aspirantura.spb.ru/forum/picture.php?albumid=111&pictureid=1004

[Hogfather]

Попробуем резюмировать.
Для подгонки модели МНК, используем код

> fit <- lm(lnZ ~ x + y, data=myData)
> summary(fit)

Call:
lm(formula = lnZ ~ x + y, data = myData)

Residuals:
Min 1Q Median 3Q Max
-0.61700 -0.30351 -0.03218 0.19850 0.81873

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.5939071 0.3067107 -1.936 0.0678 .
x 0.0016541 0.0012529 1.320 0.2025
y 0.0029021 0.0002526 11.488 5.4e-10 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.4538 on 19 degrees of freedom
Multiple R-squared: 0.889, Adjusted R-squared: 0.8773
F-statistic: 76.08 on 2 and 19 DF, p-value: 8.53e-10



Также можно вызвать следующие функции


> coefficients(fit) # коэффициенты модели
(Intercept) x y
-0.593907127 0.001654092 0.002902143
> confint(fit, level=0.95) # доверительные интервалы для параметров модели
2.5 % 97.5 %
(Intercept) -1.2358599183 0.048045664
x -0.0009683542 0.004276538
y 0.0023733723 0.003430914
> fitted(fit) # вычисленные значения
1 2 3 4 5 6 7 8
0.6184535 1.3439892 2.4322929 3.8833645 1.0225736 0.7011581 2.8777654 1.8308140
9 10 11 12 13 14 15 16
3.2818856 2.0898271 4.0487737 0.8252150 2.6390544 3.0431746 1.9962232 3.4472948
17 18 19 20 21 22
4.1728306 3.8514150 0.9906242 2.2965886 4.2555352 3.5299994
> residuals(fit) # остатки модели
1 2 3 4 5 6 7
-0.24551154 -0.03430598 0.50459127 0.03730105 -0.35167208 0.11950244 0.81873492
8 9 10 11 12 13 14
-0.61649494 0.58275224 -0.59378185 -0.03005108 0.13451833 -0.43331068 -0.31506039
15 16 17 18 19 20 21
0.02651567 -0.14198445 -0.07559024 -0.26886877 0.21982761 0.81652665 -0.61700212
22
0.46336392
> anova(fit) # дисперсионный анализ
Analysis of Variance Table

Response: lnZ
Df Sum Sq Mean Sq F value Pr(>F)
x 1 4.1589 4.1589 20.192 0.0002488 ***
y 1 27.1795 27.1795 131.963 5.395e-10 ***
Residuals 19 3.9133 0.2060
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> vcov(fit) # ковариационная матрица для параметров модели
(Intercept) x y
(Intercept) 9.407143e-02 -2.162456e-04 -4.568268e-05
x -2.162456e-04 1.569874e-06 -8.310037e-08
y -4.568268e-05 -8.310037e-08 6.382438e-08
> influence(fit) # анализ регрессионной модели
$hat
1 2 3 4 5 6 7
0.19665640 0.14809675 0.14788589 0.28317818 0.13829845 0.15275498 0.09302145
8 9 10 11 12 13 14
0.06408668 0.09164807 0.05113364 0.14270474 0.15836003 0.04654710 0.05902938
15 16 17 18 19 20 21
0.09846547 0.08567971 0.13738973 0.12649809 0.29922020 0.15799923 0.18148452
22
0.13986129

$coefficients
(Intercept) x y
1 -0.090100375 1.924032e-04 4.151880e-05
2 -0.009639360 2.941447e-05 2.351114e-06
3 0.092492681 -5.221326e-04 3.423996e-05
4 0.002356970 -5.638034e-05 1.107148e-05
5 -0.098292644 1.997504e-04 4.375226e-05
6 0.034179275 -3.504510e-05 -2.200750e-05
7 0.068581503 -4.974500e-04 6.895179e-05
8 -0.087542609 1.378093e-04 3.287562e-05
9 0.014114197 -2.636425e-04 6.738312e-05
10 -0.056330941 3.061248e-05 2.300041e-05
11 0.002092611 1.126132e-05 -6.043748e-06
12 0.026144651 5.165601e-05 -2.977428e-05
13 -0.011340740 -3.227886e-05 -3.357501e-06
14 0.009716284 -7.069650e-05 -1.206585e-05
15 0.000820820 1.626464e-05 -2.654591e-06
16 0.012887806 -5.454754e-05 -1.004487e-05
17 0.012131630 -2.194206e-05 -1.245702e-05
18 0.042158631 -1.512501e-04 -2.872849e-05
19 0.018378098 3.404810e-04 -7.109334e-05
20 -0.039976009 8.764986e-04 -8.455088e-05
21 0.143930409 -4.760290e-04 -9.195059e-05
22 -0.072988620 3.945329e-04 2.397855e-05

$sigma
1 2 3 4 5 6 7 8
0.4617764 0.4661857 0.4481136 0.4661523 0.4576378 0.4652627 0.4199357 0.4414127
9 10 11 12 13 14 15 16
0.4434361 0.4435794 0.4662052 0.4649853 0.4543848 0.4599405 0.4662215 0.4649526
17 18 19 20 21 22
0.4658732 0.4613113 0.4621416 0.4164320 0.4376836 0.4511520

$wt.res
1 2 3 4 5 6 7
-0.24551154 -0.03430598 0.50459127 0.03730105 -0.35167208 0.11950244 0.81873492
8 9 10 11 12 13 14
-0.61649494 0.58275224 -0.59378185 -0.03005108 0.13451833 -0.43331068 -0.31506039
15 16 17 18 19 20 21
0.02651567 -0.14198445 -0.07559024 -0.26886877 0.21982761 0.81652665 -0.61700212
22
0.46336392



Сравнение двух моделей

> fit1 <- lm(lnZ ~ y, data=myData)
> anova(fit,fit1)
Analysis of Variance Table

Model 1: lnZ ~ x + y
Model 2: lnZ ~ y
Res.Df RSS Df Sum of Sq F Pr(>F)
1 19 3.9133
2 20 4.2723 -1 -0.35896 1.7428 0.2025


Остановлюсь на этом месте поподробнее. Когда мы добавляем или убавляем регрессоры, мы должны быть уверены, что принципиально это на дисперсию остатков не влияет. В данном случае используется нулевая гипотеза, что коэффициенты регрессоров равны между собой. Примерно то же самое мы делаем, когда проверяем саму модель на значимость(см. ниже), но в таком случаем сравниваем её с уравнением, где все регрессоры имеют коэффициенты 0. Аналогично описанному ниже, если Pr(>F) превышает 0.05, то гипотеза принимается. Т.е. в приведенном выше примере Pr(>F)=0.2025, отсюда разница между Model 1: lnZ ~ x + y и Model 2: lnZ ~ y является статистически незначимой.

Автоматическая подгонка модели посредством исключения незначимых членов.


> library(MASS)
> step <- stepAIC(fit, direction="both")
Start: AIC=-31.99
lnZ ~ x + y

Df Sum of Sq RSS AIC
- x 1 0.359 4.2723 -32.056
<none> 3.9133 -31.987
- y 1 27.180 31.0928 11.611

Step: AIC=-32.06
lnZ ~ y

Df Sum of Sq RSS AIC
<none> 4.272 -32.056
+ x 1 0.359 3.913 -31.987
- y 1 30.979 35.252 12.372
> step$anova # отображаем результаты
Stepwise Model Path
Analysis of Deviance Table

Initial Model:
lnZ ~ x + y

Final Model:
lnZ ~ y


Step Df Deviance Resid. Df Resid. Dev AIC
1 19 3.913305 -31.98653
2 - x 1 0.3589588 20 4.272264 -32.05577


Прочую красоту смотреть вот здесь (http://www.statmethods.net/stats/regression.html). А подробно проверка регрессионной модели описана вот тут (http://www.statmethods.net/stats/rdiagnostics.html).

[Hogfather]

Хотя меня особо не спрашивают, но имея опыт работы со студентами (http://aspirantura.spb.ru/forum/showthread.php?t=10583) постараюсь кратко и без лишнего занудства изложить, что же Пытливый Ум может увидеть в сборище букв на языке эвентуального противника, когда строит модель.

Рассмотрим пример, который был приведен выше.

> summary(fit)

Call:
lm(formula = lnZ ~ x + y, data = myData)

Residuals:
Min 1Q Median 3Q Max
-0.61700 -0.30351 -0.03218 0.19850 0.81873

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.5939071 0.3067107 -1.936 0.0678 .
x 0.0016541 0.0012529 1.320 0.2025
y 0.0029021 0.0002526 11.488 5.4e-10 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.4538 on 19 degrees of freedom
Multiple R-squared: 0.889, Adjusted R-squared: 0.8773
F-statistic: 76.08 on 2 and 19 DF, p-value: 8.53e-10


Итак, первое на что мы смотрим, это не на разные Эр квадрат и прочие коэффициенты, нас интересует только вот эта строчка:
F-statistic: 76.08 on 2 and 19 DF, p-value: 8.53e-10

В данном случае нулевая гипотеза выглядит так: H0: β1=β2=...βn=0, т.е. все коэффициенты кроме свободного члена равны 0.

Это мы проводим проверку значимости модели. Если p-value у нас превышает 0.05, то дальнейшая работа с моделью просто бессмысленна, поскольку модель незначимая.



> lm0<-lm(y~x,data=test)
> summary(lm0)

Call:
lm(formula = y ~ x, data = test)

Residuals:
Min 1Q Median 3Q Max
-3.3245 -1.4632 0.6214 1.3782 2.6417

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 3.4424 0.8937 3.852 0.00117 **
x 0.1102 0.0746 1.478 0.15681
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.924 on 18 degrees of freedom
Multiple R-squared: 0.1082, Adjusted R-squared: 0.05862
F-statistic: 2.183 on 1 and 18 DF, p-value: 0.1568


http://aspirantura.spb.ru/forum/picture.php?albumid=111&pictureid=1016



Второе, на что мы смотрим, это на остатки модели (Residuals).

Residuals:
Min 1Q Median 3Q Max
-0.61700 -0.30351 -0.03218 0.19850 0.81873


Алгоритм МНК подразумевает, что матожидание остатков должно быть равно 0 и мы надеемся, что остатки имеют нормальное распределение. Напоминаю, что у нормального распределения матожидание и медиана совпадают. В данном случае -- не очень. Поскольку мы видим отрицательную медиану (-0.03218), то распределение остатков незначительно смещено влево. Квартили, по идее, должны быть симметричны относительно 0, здесь же они значительно смещены влево. Максимум и минимум позволяют оценить выбросы (http://ru.wikipedia.org/wiki/%D0%92%D1%8B%D0%B1%D1%80%D0%BE%D1%81_(%D1%81%D1%82 %D0%B0%D1%82%D0%B8%D1%81%D1%82%D0%B8%D0%BA%D0%B0)) (outliers) модели.

Третье, на что мы смотрим -- коэффициенты.

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.5939071 0.3067107 -1.936 0.0678 .
x 0.0016541 0.0012529 1.320 0.2025
y 0.0029021 0.0002526 11.488 5.4e-10 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


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

Дальше, Residual standard error: 1.924 on 18 degrees of freedom говорит нам о средней квадратической ошибке остатков (σ).

Multiple R-squared: 0.889, Adjusted R-squared: 0.8773 говорит о коэффициентах детерминации и качестве модели. Чем они ближе к единице, тем лучше.

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


> test<-data.frame(x<-1:20)
> test$y<-sqrt(test$x)
> lm0<-lm(y~x,data=test)
> summary(lm0)

Call:
lm(formula = y ~ x, data = test)

Residuals:
Min 1Q Median 3Q Max
-0.46467 -0.08467 0.04310 0.13672 0.17227

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.294290 0.081968 15.79 5.44e-12 ***
x 0.170382 0.006843 24.90 2.13e-15 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.1765 on 18 degrees of freedom
Multiple R-squared: 0.9718, Adjusted R-squared: 0.9702
F-statistic: 620 on 1 and 18 DF, p-value: 2.127e-15

Обратите внимание на R2
А теперь смотрите на график

http://aspirantura.spb.ru/forum/picture.php?albumid=111&pictureid=1017

Что должно быть на графике:
1. Точки на графике Residuals vs Fitted должны быть случайно разбросаны, без явно выраженного тренда. Тут этого нет.
2. Точки на графике Normal Q-Q должны более-менее лежать на линии, показывая что остатки имеют нормальное распредение
3. Точки на двух остальных графиках должны быть в группах не слишком далеко от центра.

Т.е. модель у нас неадекватна, не смотря на милый Эр Квадрат. Такие дела.
Для сравнения посмотрите график на предыдущей странице (http://aspirantura.spb.ru/forum/album.php?albumid=111&pictureid=1003)



Добавлено через 4 часа 3 минуты

Использование метода Бокса-Кокса для преобразования регрессионной модели

В предыдущем примере мы рассмотрели случай, когда модель выбрана неадекватной. Судя по параболе, которую мы наблюдаем, здесь необходимо какое-то степенное преобразование.
Попробуем использовать преобразование Бокса-Кокса, которое входит в пакет MASS. Т.е. найдем такую λ, чтобы yλ сделал нам аж боже мой.

Исходные данные у нас те же и та же модель.

> test<-data.frame(x<-1:20)
> test$y<-sqrt(test$x)
> lm0<-lm(y~x,data=test)


А дальше начинается колдунство. В данном случае lambda взята в интервале 0 до 4 с шагом 0.01 ( по умолчанию -2 до 2 с шагом 0.1)

> library(MASS)
> bc<-boxcox(lm0,lambda = seq(0,4,0.01))


Получаем следующий график зависимости логарифмической функции правдоподобия от лямбды.
http://aspirantura.spb.ru/forum/picture.php?albumid=111&pictureid=1018

К сожалению, конкретного результата, который нужно подставить у нас нигде нет. Но мы можем его получить. Функция which.max(bc$y) вернёт нам номер элемента, у которого Игрек наибольший, а bc$x[which.max(bc$y)]) -- собственно нужную лямбду.

> (lambda<-bc$x[which.max(bc$y)])
[1] 2

Лямбда получилась равной 2. Ах, удивительно!
Делаем преобразование и строим новую модель.
> lm1<-lm(I(y^lambda)~x,data=test)
> summary(lm1)

Call:
lm(formula = I(y^lambda) ~ x, data = test)

Residuals:
Min 1Q Median 3Q Max
-7.049e-15 -1.023e-15 -1.436e-16 1.235e-15 7.585e-15

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.000e+00 1.324e-15 0.000e+00 1
x 1.000e+00 1.105e-16 9.047e+15 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.85e-15 on 18 degrees of freedom
Multiple R-squared: 1, Adjusted R-squared: 1
F-statistic: 8.185e+31 on 1 and 18 DF, p-value: < 2.2e-16

Нарисуем первый график из четырех диагностических графиков модели.


> plot(lm1,which=1)


http://aspirantura.spb.ru/forum/picture.php?albumid=111&pictureid=1019

Воот! Теперь совсем другое дело.

Следует отметить особые случаи в преобразовании Бокса-Кокса. Если λ=1, то преобразование не требуется. Если λ=0, то следует использовать логарифмическое преобразование log(y). Пример.


> test<-data.frame(x<-1:20)
> test$y<-exp(-test$x)
> lm0<-lm(y~x,data=test)
> summary(lm0)

Call:
lm(formula = y ~ x, data = test)

Residuals:
Min 1Q Median 3Q Max
-0.06529 -0.04258 -0.01349 0.02358 0.26464

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.111047 0.034464 3.222 0.00473 **
x -0.007805 0.002877 -2.713 0.01426 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.07419 on 18 degrees of freedom
Multiple R-squared: 0.2902, Adjusted R-squared: 0.2508
F-statistic: 7.359 on 1 and 18 DF, p-value: 0.01426

> bc<-boxcox(lm0,lambda = seq(-4,4,0.01))
> (lambda<-bc$x[which.max(bc$y)])
[1] 0
> lm1<-lm(log(y)~x,data=test)
> summary(lm1)

Call:
lm(formula = log(y) ~ x, data = test)

Residuals:
Min 1Q Median 3Q Max
-2.038e-15 -6.074e-16 -3.526e-16 -8.400e-18 6.208e-15

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 4.766e-15 7.696e-16 6.193e+00 7.6e-06 ***
x -1.000e+00 6.425e-17 -1.557e+16 < 2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.657e-15 on 18 degrees of freedom
Multiple R-squared: 1, Adjusted R-squared: 1
F-statistic: 2.423e+32 on 1 and 18 DF, p-value: < 2.2e-16



И для закрепления материала два графика. Слева -- плохой, негодный, а справа -- хороший, годный.



> par(mfrow=c(1,2))
> plot(lm0,which=1)
> plot(lm1,which=1)
> par(mfrow=c(1,1))



http://aspirantura.spb.ru/forum/picture.php?albumid=111&pictureid=1020

Добавлено через 18 часов 12 минут

Кореллограммы
В рассмотренной задаче (#10 (http://aspirantura.spb.ru/forum/showpost.php?p=294873&postcount=10)) рисовалась корреляционная матрица в виде соответствующих диаграмм рассеяния. Также можно получить корреляционную матрицу в виде коэффициентов корреляции в текстовом виде.


> # Инициализация данных
> x <- c(75, 75, 75, 75, 100, 125, 125, 150, 150, 175, 175, 200, 200, 225, 250, 250, 250, 275, 300, 300, 300, 300)
> y <- c(375, 625, 1000, 1500, 500, 375, 1125, 750, 1250, 825, 1500, 375, 1000, 1125, 750, 1250, 1500, 1375, 375, 825, 1500, 1250)
> z <- c(1.452, 3.705, 18.857, 50.434, 1.956, 2.272, 40.306, 3.368, 47.686, 4.464, 55.63, 2.611, 9.077, 15.304, 7.559, 27.257, 60.174, 35.965, 3.355, 22.491, 38.036, 54.237)
> myData<-data.frame(x,y,z)
> myData$lnZ<-log(myData$z)
>
> # Собственно, вывод корреляционной матрицы
> cor(myData[,c("x","y","z")])
x y z
x 1.0000000 0.2625290 0.2394832
y 0.2625290 1.0000000 0.8857714
z 0.2394832 0.8857714 1.0000000


Обратите внимание, мы взяли не 4 столбца таблицы, а только три. Столбец lnZ опущен. Хотя данная матрица неплохо читается, её можно сделать более удобной для восприятия с помощью функции symnum, заменяющей числовые значения на символы.
> symnum(cor(myData[,c("x","y","z")]))
x y z
x 1
y 1
z + 1
attr(,"legend")
[1] 0 ‘ ’ 0.3 ‘.’ 0.6 ‘,’ 0.8 ‘+’ 0.9 ‘*’ 0.95 ‘B’ 1

Для сравнения

> symnum(cor(myData[,c("x","y","lnZ")]))
x y l
x 1
y 1
lnZ . * 1
attr(,"legend")
[1] 0 ‘ ’ 0.3 ‘.’ 0.6 ‘,’ 0.8 ‘+’ 0.9 ‘*’ 0.95 ‘B’ 1


Однако, есть и другой способ сделать это, да и к тому же красиво, нарисовав коррелограмму. Такой рисунок несказанно обогатит статью и поразит всех читателей. Для этого нужна установленная библиотека corrgram.
Примеры вызовов функции приведены ниже.

Существует 4 возможных отображения корреляции:
panel.pts -- диаграмма разброса
panel.pie -- круговая диаграмма, меняет цвет и заполненность в зависимость от коэффициента корреляции (отрицательная -- красная, положительная -- синяя).
panel.ellipse -- рисует красоту, как показана на рисунке 3. В общем, границы и аппроксимацию.
panel.shade -- закрашивает прямоугольник в зависимости от коэффициента корреляции.

Их присваивают, соответственно lower.panel или upper.panel, в зависимости от того, где нужен рисунок: сверху или снизу.
Также можно выводить название переменных и минимум/максимум значений, как это показано ниже.


>
> library(corrgram)
>
>
> corrgram(myData[,c("x","y","z")], lower.panel=panel.pie,
+ upper.panel=panel.pts, text.panel=panel.txt,
+ diag.panel=panel.minmax, main="Кореллограмма 1")
>
http://aspirantura.spb.ru/forum/picture.php?albumid=111&pictureid=1021
> corrgram(myData[,c("x","y","lnZ")], lower.panel=panel.pts,
+ upper.panel=panel.pie, text.panel=panel.txt,
+ diag.panel=panel.minmax, main="Кореллограмма 2")
>

http://aspirantura.spb.ru/forum/picture.php?albumid=111&pictureid=1022

> corrgram(log(myData[,c("x","y","z")]), order=TRUE, lower.panel=panel.shade,
+ upper.panel=panel.ellipse, text.panel=panel.txt,
+ main="Кореллограмма 3\n Логарифмы величин")

http://aspirantura.spb.ru/forum/picture.php?albumid=111&pictureid=1023

Понятно, что на трех переменных это не сильно показательно, но если их будет десяток, то картинка будет гораздо симпатичнее.

[Hogfather]

Чем хорош R, что в нём можно много всякого забавного сделать. например, захотелось поиграть в регрессионные деревья (http://ru.wikipedia.org/wiki/%D0%94%D0%B5%D1%80%D0%B5%D0%B2%D0%BE_%D0%BF%D1%80% D0%B8%D0%BD%D1%8F%D1%82%D0%B8%D1%8F_%D1%80%D0%B5%D 1%88%D0%B5%D0%BD%D0%B8%D0%B9) -- никто слова плохого не скажет.

Рассмотрим ещё раз нашу задачу.



> library(rpart)
> (fit <- rpart(lnZ ~ x + y, data = myData))
n= 22

node), split, n, deviance, yval
* denotes terminal node

1) root 22 35.251650 2.508130
2) y< 912.5 10 5.460743 1.319059 *
3) y>=912.5 12 3.869613 3.499022 *
> summary(fit)
Call:
rpart(formula = lnZ ~ x + y, data = myData)
n= 22

CP nsplit rel error xerror xstd
1 0.7353214 0 1.0000000 1.0863277 0.1579386
2 0.0100000 1 0.2646786 0.7911578 0.1943923

Variable importance
y x
91 9

Node number 1: 22 observations, complexity param=0.7353214
mean=2.50813, MSE=1.602348
left son=2 (10 obs) right son=3 (12 obs)
Primary splits:
y < 912.5 to the left, improve=0.7353214, (0 missing)
x < 212.5 to the left, improve=0.1397332, (0 missing)
Surrogate splits:
x < 112.5 to the left, agree=0.591, adj=0.1, (0 split)

Node number 2: 10 observations
mean=1.319059, MSE=0.5460743

Node number 3: 12 observations
mean=3.499022, MSE=0.3224678

> par(xpd = NA) # сильное колдунство, чтобы текст не обрезался
> plot(fit)
> text(fit, use.n = TRUE)


http://aspirantura.spb.ru/forum/picture.php?albumid=111&pictureid=1024

Результат, как видите следующий. Если y<912.5, то ln(z) равен 1.319, а если y>=912.5, то 3.499. Смысла для данной задачи никакого, разве что лишний раз убедились, что Икс у нас не при делах.



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

В качестве исходных данных берется таблица данных по кифозу у детей из 81 строк and 4 столбцов, описывающая детей у которых были корректирующие операции на позвоночнике

Таблица содержит следующие колонки:

Kyphosis - качественные данные (factor) с вариантами "absent" и "present", показывающая был ли кифоз после операции.
Age- возраст в месяцах
Number - число затронутых позвонков
Start - номер верхнего позвонка, на котором происходила операция

> # Первые строки таблицы
> head(kyphosis)
Kyphosis Age Number Start
1 absent 71 3 5
2 absent 158 3 14
3 present 128 4 5
4 absent 2 5 1
5 absent 1 4 15
6 absent 1 2 16
> (fit <- rpart(Kyphosis ~ Age + Number + Start, data = kyphosis))
n= 81

node), split, n, loss, yval, (yprob)
* denotes terminal node

1) root 81 17 absent (0.79012346 0.20987654)
2) Start>=8.5 62 6 absent (0.90322581 0.09677419)
4) Start>=14.5 29 0 absent (1.00000000 0.00000000) *
5) Start< 14.5 33 6 absent (0.81818182 0.18181818)
10) Age< 55 12 0 absent (1.00000000 0.00000000) *
11) Age>=55 21 6 absent (0.71428571 0.28571429)
22) Age>=111 14 2 absent (0.85714286 0.14285714) *
23) Age< 111 7 3 present (0.42857143 0.57142857) *
3) Start< 8.5 19 8 present (0.42105263 0.57894737) *
> (fit2 <- rpart(Kyphosis ~ Age + Number + Start, data = kyphosis,
+ parms = list(prior = c(.65,.35), split = "information")))
n= 81

node), split, n, loss, yval, (yprob)
* denotes terminal node

1) root 81 28.350000 absent (0.65000000 0.35000000)
2) Start>=12.5 46 3.335294 absent (0.91563089 0.08436911) *
3) Start< 12.5 35 16.453120 present (0.39676840 0.60323160)
6) Age< 34.5 10 1.667647 absent (0.81616742 0.18383258) *
7) Age>=34.5 25 9.049219 present (0.27932897 0.72067103) *
> par(mfrow = c(1,2), xpd = NA) # otherwise on some devices the text is clipped
> plot(fit)
> text(fit, use.n = TRUE)
> plot(fit2)
> text(fit2, use.n = TRUE)
> par(mfrow = c(1,1))
> ############################################
> summary(fit)
Call:
rpart(formula = Kyphosis ~ Age + Number + Start, data = kyphosis)
n= 81

CP nsplit rel error xerror xstd
1 0.17647059 0 1.0000000 1.0000000 0.2155872
2 0.01960784 1 0.8235294 0.9411765 0.2107780
3 0.01000000 4 0.7647059 1.0000000 0.2155872

Variable importance
Start Age Number
64 24 12

Node number 1: 81 observations, complexity param=0.1764706
predicted class=absent expected loss=0.2098765 P(node) =1
class counts: 64 17
probabilities: 0.790 0.210
left son=2 (62 obs) right son=3 (19 obs)
Primary splits:
Start < 8.5 to the right, improve=6.762330, (0 missing)
Number < 5.5 to the left, improve=2.866795, (0 missing)
Age < 39.5 to the left, improve=2.250212, (0 missing)
Surrogate splits:
Number < 6.5 to the left, agree=0.802, adj=0.158, (0 split)

Node number 2: 62 observations, complexity param=0.01960784
predicted class=absent expected loss=0.09677419 P(node) =0.7654321
class counts: 56 6
probabilities: 0.903 0.097
left son=4 (29 obs) right son=5 (33 obs)
Primary splits:
Start < 14.5 to the right, improve=1.0205280, (0 missing)
Age < 55 to the left, improve=0.6848635, (0 missing)
Number < 4.5 to the left, improve=0.2975332, (0 missing)
Surrogate splits:
Number < 3.5 to the left, agree=0.645, adj=0.241, (0 split)
Age < 16 to the left, agree=0.597, adj=0.138, (0 split)

Node number 3: 19 observations
predicted class=present expected loss=0.4210526 P(node) =0.2345679
class counts: 8 11
probabilities: 0.421 0.579

Node number 4: 29 observations
predicted class=absent expected loss=0 P(node) =0.3580247
class counts: 29 0
probabilities: 1.000 0.000

Node number 5: 33 observations, complexity param=0.01960784
predicted class=absent expected loss=0.1818182 P(node) =0.4074074
class counts: 27 6
probabilities: 0.818 0.182
left son=10 (12 obs) right son=11 (21 obs)
Primary splits:
Age < 55 to the left, improve=1.2467530, (0 missing)
Start < 12.5 to the right, improve=0.2887701, (0 missing)
Number < 3.5 to the right, improve=0.1753247, (0 missing)
Surrogate splits:
Start < 9.5 to the left, agree=0.758, adj=0.333, (0 split)
Number < 5.5 to the right, agree=0.697, adj=0.167, (0 split)

Node number 10: 12 observations
predicted class=absent expected loss=0 P(node) =0.1481481
class counts: 12 0
probabilities: 1.000 0.000

Node number 11: 21 observations, complexity param=0.01960784
predicted class=absent expected loss=0.2857143 P(node) =0.2592593
class counts: 15 6
probabilities: 0.714 0.286
left son=22 (14 obs) right son=23 (7 obs)
Primary splits:
Age < 111 to the right, improve=1.71428600, (0 missing)
Start < 12.5 to the right, improve=0.79365080, (0 missing)
Number < 4.5 to the left, improve=0.07142857, (0 missing)

Node number 22: 14 observations
predicted class=absent expected loss=0.1428571 P(node) =0.1728395
class counts: 12 2
probabilities: 0.857 0.143

Node number 23: 7 observations
predicted class=present expected loss=0.4285714 P(node) =0.08641975
class counts: 3 4
probabilities: 0.429 0.571



http://aspirantura.spb.ru/forum/picture.php?albumid=111&pictureid=1025

[Olafson]

Hogfather,

мне бывает нужен МНК на очень подробных сетках (1000\times 1000, например). Потянет это оболочка? (без ограничения производительности железа)

[Hogfather]

мне бывает нужен МНК на очень подробных сетках (1000\times 1000, например). Потянет это оболочка? (без ограничения производительности железа)
Специально для Вас прогнал простой тест. Помучался с подобром матрицы, поскольку сингулярные матрицы решаются на раз.
Обратите внимание, я специально выводил время. 13 секунд на модель. Дольше таблицу инициализировал...
Процесcор старый, Intel Core Duo T7250 2 ГГц. Озу 3 Гб



> # Инициализируем матрицу 1000х1000 для x1..x1000, с матожидаением i
> # y присваеваем вектор в 1000 элементов
> Sys.time()
[1] "2012-12-11 12:01:09 MSK"
> MyData<-data.frame(y=rep(0,1000))
> for(i in 1:1000) {
+ MyData[paste("x",i,sep="")]<-(1:1000)^(1/i)
+ MyData["y"]<-MyData["y"]+MyData[paste("x",i,sep="")]*i+rnorm(10)
+ }
>
> # Понеслась
> Sys.time()
[1] "2012-12-11 12:01:49 MSK"
> test<-lm(y~.,data=MyData)
> Sys.time()
[1] "2012-12-11 12:02:03 MSK"
> summary(test)

Call:
lm(formula = y ~ ., data = MyData)

Residuals:
Min 1Q Median 3Q Max
-55.804 -16.049 4.162 14.695 61.563

Coefficients: (993 not defined because of singularities)
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.290e+06 2.890e+06 0.446 0.655
x1 -2.337e+00 1.785e+01 -0.131 0.896
x2 1.271e+04 5.447e+04 0.233 0.815
x3 -6.552e+05 2.619e+06 -0.250 0.803
x4 5.647e+06 2.185e+07 0.258 0.796
x5 -1.424e+07 5.407e+07 -0.263 0.792
x6 1.048e+07 3.930e+07 0.267 0.790
x7 NA NA NA NA
x8 NA NA NA NA
x9 NA NA NA NA
x10 NA NA NA NA
x11 NA NA NA NA
x12 NA NA NA NA
x13 NA NA NA NA
x14 NA NA NA NA
x15 NA NA NA NA
x16 NA NA NA NA
x17 NA NA NA NA
x18 NA NA NA NA
x19 NA NA NA NA
x20 NA NA NA NA
x21 -2.029e+06 7.406e+06 -0.274 0.784
x22 NA NA NA NA
x23 NA NA NA NA
x24 NA NA NA NA
x25 NA NA NA NA
x26 NA NA NA NA
x27 NA NA NA NA
x28 NA NA NA NA
x29 NA NA NA NA
.....
x993 NA NA NA NA
x994 NA NA NA NA
x995 NA NA NA NA
x996 NA NA NA NA
x997 NA NA NA NA
x998 NA NA NA NA
x999 NA NA NA NA
x1000 NA NA NA NA

Residual standard error: 30.74 on 992 degrees of freedom
Multiple R-squared: 0.9994, Adjusted R-squared: 0.9994
F-statistic: 2.538e+05 on 7 and 992 DF, p-value: < 2.2e-16

[Olafson]

Hogfather, спасибо!

Синтаксис не точно понимаю, но я понял, что R справился. Как он отдает результаты?

[Hogfather]

Как он отдает результаты?
Можно конкретизировать вопрос: а как хотите?

[Olafson]

а как хотите?

Скорее всего, он визуализирует и препарирует любые агрегаты из данных сам.
Но может понадобиться что-нибудь сугубо универсальное, типа текстового файла, для демонстрации результатов через всякие Excel'и с Maple'ами. (ведь R, как я понял, -- GNU, а те далеко нет)

[Hogfather]

Но может понадобиться что-нибудь сугубо универсальное, типа текстового файла
Функция write.table() позволяет сохранить нужные данные в формате csv. А дальше уже дело техники.

Дублирую табличку из темы про GRETL.
Откуда берем данные|Формат обмена|Куда перегоняем|Примечание
Excel|xls, xlsx или csv|Gretl| Читается прямо в формате Excel
Excel|csv|GNU R |Можно и напрямую, но с доп. библиотекой, а также через clipboard (см. read.table("clipboard",header=T) (http://aspirantura.spb.ru/forum/showpost.php?p=290293&postcount=3) )
Gretl|csv|Excel| Используйте правильный разделитель данных. Лучше ";". Например, сохраняем таблицу MyData: write.table(MyData,"myData.csv", sep = ";", dec = "."). Для небольшого объема данных можно через clipboard, заменив имя файла на "clipboard" и использовав sep = "\t" (табулятор).
Gretl|GNU R (.RData), csv|GNU R| Лучше использовать внутренний формат .RData, потом в R воспользоваться командой load("имя файла")
GNU R|csv,xls|Excel|Используйте правильный разделитель данных. Лучше "\t". Для экспорта в Excel можно использовать библиотеку xlsReadWrite
GNU R|csv|gretl|В случае, если запущена сессия из gretl, можно также воспользоваться командой gretl.export(). Подробности см. в руководстве пользователя раздел "Gretl and R"

[Olafson]

Hogfather, merci за консультацию.

[Hogfather]

Как выяснилось в переписке, вопрос (http://aspirantura.spb.ru/forum/showthread.php?p=301094#post301094) касался также 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 (http://ru.wikipedia.org/wiki/RESET-%D1%82%D0%B5%D1%81%D1%82_%D0%A0%D0%B0%D0%BC%D1%81% D0%B5%D1%8F).
Естественно, что он есть в 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 (http://alexandre-putt.livejournal.com/9837.html)

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

Во-первых, пакет называется "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]

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

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

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

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

P.S. На большинство вопросов есть ответ вот тут http://cran.rstudio.com/bin/windows/base/rw-FAQ.html#What_0027s-the-best-way-to-upgrade_003f
Достаточно быстро проходит операция и все пакеты переносятся. Проверял.

[Hogfather]

R Cookbook (O'Reilly Cookbooks)
http://ecx.images-amazon.com/images/I/51ODN7aoMqL._BO2,204,203,200_PIsitb-sticker-arrow-click,TopRight,35,-76_AA300_SH20_OU01_.jpg (http://www.amazon.com/gp/reader/0596809158/ref=sib_dp_pt#reader-link)

Книга из серии 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/I/31IUhlByl1L._BO2,204,203,200_PIsitb-sticker-arrow-click-small,TopRight,12,-30_AA300_SH20_OU01_.jpg (http://www.amazon.com/gp/reader/0521762936/ref=sib_dp_pt#reader-link)

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

Оценка: 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/I/51Uj0FIObeL._BO2,204,203,200_PIsitb-sticker-arrow-click,TopRight,35,-76_AA300_SH20_OU01_.jpg (http://www.amazon.com/gp/reader/1439810184/ref=sib_dp_pt#reader-link)

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

Оценка: 3 из 5

[Hogfather]

R Graphics Cookbook [Paperback]

http://ecx.images-amazon.com/images/I/51gjgnl23VL._BO2,204,203,200_PIsitb-sticker-arrow-click,TopRight,35,-76_AA300_SH20_OU01_.jpg (http://www.amazon.com/gp/reader/1449316956/ref=sib_dp_pt#reader-link)

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

Оценка: 5 из 5


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

http://ecx.images-amazon.com/images/I/51N1ztEPh%2BL._BO2,204,203,200_PIsitb-sticker-arrow-click,TopRight,35,-76_AA300_SH20_OU01_.jpg (http://www.amazon.com/gp/reader/0123814855/ref=sib_dp_pt#reader-link)

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

Оценка: 5 из 5

[Hogfather]

Задал мне тут вопрос Дмитрий В., говорит, Хогфазер, а есть ли в R лица Чернова (http://nordisk.pp.ru/44/), зело они мне на конференции понравились, хочу научному руководителю рожу гнусную состроить, но на научной основе?

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

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

> 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/picture.php?albumid=111&pictureid=1092

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

http://aspirantura.spb.ru/forum/picture.php?albumid=111&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/picture.php?albumid=111&pictureid=1094

[Hogfather]

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

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

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

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

Обзор возможностей 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 (http://johnmacfarlane.net/pandoc/). Установив оную и запустив из командной строки в нужном каталоге можно получить на выходе Word документ.
Запуск выглядит так

pandoc -o MarkdownDemo.docx MarkdownDemo.md


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

[Hogfather]

1. Трансформация данных

Возник вопрос: как преобразовать данные в R на манер сводных таблиц Excel. Вопрос был решен, но в процессе открыты для себя два полезных пакета, о которых хотелось бы рассказать.

Пакет reshape2 - позволяет агрегировать данные на манер Excel. На сайте http://www.r-statistics.com/tag/aggregate/ есть неплохой разбор его возможностей. Основная идея должна быть понятна вот из этой схемы (кликабельна).

http://www.r-statistics.com/wp-content/uploads/2012/01/reshaping-data-using-melt-and-cast-300x214.png (http://www.r-statistics.com/wp-content/uploads/2012/01/reshaping-data-using-melt-and-cast.png)

Второй пакет -- sqldf. Позволяет писать прямо в коде SQL запросы к данным в синтаксисе SQLite.

Привожу простой код, который решает одну и ту же задачу с использованием двух этих пакетов. Имеются данные по весу кошачьих сердец, кошачьих тушек и пола (Пакет MASS данные cats), Попробуем найти для каждого пола число измерений и средний вес тушки.

Напишем тестовый пример, который бы в цикле 1000 раз пытался решить эту задачу и посчитаем затраченное время.

Возможное решение.

> library(MASS)
> library(sqldf)
> library(reshape2)
>
> head(cats)
Sex Bwt Hwt
1 F 2.0 7.0
2 F 2.0 7.4
3 F 2.0 9.5
4 F 2.1 7.2
5 F 2.1 7.3
6 F 2.1 7.6
>
> # Первый тест тест. Библиотека sqldf
> x<-Sys.time()
> for(i in 1:1000) z<-sqldf("select Sex,count(*) as cnt,avg(Bwt) as Avg_Bwt from cats group by Sex")
> y<-Sys.time()
> y-x
Time difference of 19.61596 secs
> z
Sex cnt Avg_Bwt
1 F 47 2.359574
2 M 97 2.900000
> # Второй тест. Библиотека reshape 2
> id<-1:nrow(cats)
> mcats<-melt(cbind(id,cats),id=c("id","Sex"))
> x<-Sys.time()
> for(i in 1:1000) {
+ df1<-dcast(subset(mcats,variable == "Bwt"),formula=Sex~variable,length)
+ df2<-dcast(subset(mcats,variable == "Bwt"),formula=Sex~variable,mean)
+ merge(df1,df2,intersect="Sex")
+ }
> y<-Sys.time()
> y-x
Time difference of 12.49125 secs
> z
Sex cnt Avg_Bwt
1 F 47 2.359574
2 M 97 2.900000


Как видно, второй вариант работает гораздо быстрее, но первый проще писать. Так что однозначного выбора нет.

Во втором примере также показано объединение двух таблиц по ключевому полю с использованием команды merge, поскольку команда cast не позволяет агрегировать сразу по двум функциям (может и позволяет, но я не умею).
Такие дела (с).


2. Графы в GNU R

Рассматривая интересные библиотеки в R можно упомянуть о возможности строить графы с помощью библиотеки igraph.

Простой пример кода. Пытаемся построить граф связей на портале аспирантов (фрагмент).

library(graph)
# Заполнем фрейм парами данных
left<-rep("Hogfather",16)
right<-c("Alextiger","caty-zharr","Dukar","Ink","IvanSpbRu","Martusya","osmos","saovu","Seta","Uzanka","Vica3","Димитриадис","Дмитрий В.","Домохозяйка","море","Степан Капуста")
left1<-rep("Степан Капуста",2)
right1<-c("Hogfather","Ink")
left<-c(left,left1)
right<-c(right,right1)
left1<-rep("море",10)
right1<-c("agasfer","Alextiger","bugo","fazotron","Hogfather","Ink","IvanSpbRu","Maksimus","Martusya","osmos")
left<-c(left,left1)
right<-c(right,right1)
myData<-data.frame(left,right)

# Формируем и строим граф
g<-graph.data.frame(myData,directed=F)
plot(g,vertex.label.cex=0.8,vertex.label.dist=1,ve rtex.size=10)


Результат

http://www.aspirantura.spb.ru/forum/picture.php?albumid=111&pictureid=1100

Нужно или нет, решать вам, но если надо рисовать серьезные графы, то рекомендую сперва глянуть в сторону graphviz (http://www.graphviz.org). Очень кошерная и мощная штука.

[Hogfather]

Вышла книга. Шитиков В.К., Розенберг Г.С. Рандомизация и бутстреп: статистический анализ в биологии и экологии с использованием R. - Тольятти: «Кассандра», 2013. - 289 с. (http://www.ievbras.ru/ecostat/Kiril/Article/A32/Starb.pdf)

Ссылка ведет на Интернет-версию на сайте авторов.

[Hogfather]

Добрый день!
Опишу в кратце иеющиеся данные. В моём распоряжении есть пара тысяч текстовых файлов, в каждом из которых десяток тысяч строк вида: 1.2.2013 12:12 1234 1234 234 23423 (т.е. дата, время, четыре числа). В файлах содержится статистическая информация о работе некого сложного технического устройства. Один файл - это статистика работы одного устройства за период в несколько лет с интервалом в несколько минут.
Сейчас я провёл анализ данных десятка файлов. Определил несколько интересных закономерностей, для описания которых создал мат.модель. Однако в ручную обработать более 2 000 файлов - процесс долгий, т.е. нужна автоматизация. Плюс есть множество идей по выявлению других закономерностей Анализ нужно проводить как по каждому отдельному файлу, т.к. и по группе файлов.

Не рассматривая собственно вопрос анализа, продемонстрирую как можно закачать данные в пакет R. Для тестового примера был создан каталог Data_R, в которые положены 2 одинаковых файла test1.txt и test2.txt, вида:


1.2.2013 12:12 1234 1234 234 23423
1.2.2013 12:12 1234 1234 234 23423
1.2.2013 12:12 1234 1234 234 23423
1.2.2013 12:12 1234 1234 234 23423
1.2.2013 12:12 1234 1234 234 23423
1.2.2013 12:12 1234 1234 234 23423


Простенький скрипт, который создает список всех txt файлов в данном каталоге, а потом импортирует их во внутреннюю таблицу R, добавляя колонку "Источник" прилагаю.


FileList<-dir("Data_R",pattern = "*.txt", full.names = TRUE, ignore.case = TRUE, include.dirs =TRUE)
MyResult<-data.frame()
for(FileName in FileList)
{
TempTable<-read.table(FileName,sep=" ")
TempTable$FileName<-FileName
MyResult<-rbind(MyResult,TempTable)
}


Для просмотре результатов работы вводит команду MyResult, которая, собственно, отобразит все значения этой таблицы


> MyResult
V1 V2 V3 V4 V5 V6 FileName
1 1.2.2013 12:12 1234 1234 234 23423 Data_R/test1.txt
2 1.2.2013 12:12 1234 1234 234 23423 Data_R/test1.txt
3 1.2.2013 12:12 1234 1234 234 23423 Data_R/test1.txt
4 1.2.2013 12:12 1234 1234 234 23423 Data_R/test1.txt
5 1.2.2013 12:12 1234 1234 234 23423 Data_R/test1.txt
6 1.2.2013 12:12 1234 1234 234 23423 Data_R/test1.txt
7 1.2.2013 12:12 1234 1234 234 23423 Data_R/test2.txt
8 1.2.2013 12:12 1234 1234 234 23423 Data_R/test2.txt
9 1.2.2013 12:12 1234 1234 234 23423 Data_R/test2.txt
10 1.2.2013 12:12 1234 1234 234 23423 Data_R/test2.txt
11 1.2.2013 12:12 1234 1234 234 23423 Data_R/test2.txt
12 1.2.2013 12:12 1234 1234 234 23423 Data_R/test2.txt


Как видим, все работает. Дальше можно уже делать срезы какие душе угодно.

Добавлено через 36 минут
С другой стороны, данных много, поэтому возможно удобнее будет загнать в базу данных SQLite (см. выше о работе в R (http://aspirantura.spb.ru/forum/showthread.php?t=10501&page=3#22)) и использовать срезы оттуда.

Модифицируем чуть-чуть вышеприведенный код

require(RSQLite)
drv <- dbDriver("SQLite")
fname<-"mytest.sqlite"
conn <- dbConnect(drv, dbname = fname)
if(dbExistsTable(conn, "RESULTS")) {
dbRemoveTable(conn, "RESULTS")
}
FileList<-dir("Data_R",pattern = "*.txt", full.names = TRUE, ignore.case = TRUE, include.dirs =TRUE)
for(FileName in FileList)
{
TempTable<-read.table(FileName,sep=" ")
TempTable$FileName<-FileName
if(dbExistsTable(conn, "RESULTS")) {
dbWriteTable(conn, "RESULTS", TempTable, append = T)
}
else
{
dbWriteTable(conn, "RESULTS", TempTable)
}
}
dbDisconnect(conn)
dbUnloadDriver(drv)


И посмотрим результат

http://aspirantura.spb.ru/forum/picture.php?albumid=111&pictureid=1162

Добавлено через 49 минут
В какой программе можно построить такие закрашенные полосы?
Пара кривых строится по точкам x-y, и область между ними закрашивается.

Kayra, Пример кода.


set.seed(666)
age <- 1:10
y.low <- rnorm(length(age), 150, 25) + 10*age
y.high <- rnorm(length(age), 250, 25) + 10*age
plot(age,y.high,type = 'n', ylim = c(100, 400), ylab = 'Y Range', xlab = 'Age (years)')
lines(age, y.low )
lines(age, y.high)
polygon(c(age, rev(age)), c(y.high, rev(y.low)), col = "lightblue", border = NA)

# Второй график, чтобы поржать

y.low <- rnorm(length(age), 150, 25) + 10*age
y.high <- rnorm(length(age), 250, 25) + 10*age
lines(age, y.low, col = "green" )
lines(age, y.high, col = "green")
polygon(c(age, rev(age)), c(y.high, rev(y.low)), col = "green", border = NA)



http://aspirantura.spb.ru/forum/picture.php?albumid=111&pictureid=1163

[Uzanka]

Hogfather,
а в R можно оценивать GARCH модели с разными распределениями ошибок?

а оценивать Stochastic volatility models с разными распределениями?

Или самим код надо писать?

ЗЫ. А что-то типа фильтра Калмана там есть?

[Hogfather]

а в R можно оценивать GARCH модели с разными распределениями ошибок?
Для примера потребуются две библиотеки fGarch и evir
Пример из книжки, упомянутой ниже


library(fGarch)
data(bmw,package="evir")

bmw.garch_norm = garchFit(~arma(1,0)+garch(1,1),data=bmw,cond.dist="norm")
options(digits=3)
summary(bmw.garch_norm)
options(digits=10)


x = bmw.garch_norm@residuals / bmw.garch_norm@sigma.t
n=length(bmw)
grid = (1:n)/(n+1)
fitdistr(x,"t")

par(mfrow=c(1,2))
qqnorm(x,datax=T,ylab= "Standardized residual quantiles",
main="(a) normal plot",
xlab="normal quantiles")
qqline(x,datax=T)
qqplot(sort(x), qt(grid,df=4),
main="(b) t plot, df=4",xlab= "Standardized residual quantiles",
ylab="t-quantiles")
abline( lm( qt(c(.25,.75),df=4)~quantile(x,c(.25,.75)) ) )

bmw.garch_t = garchFit(~arma(1,1)+garch(1,1),cond.dist="std",data=bmw)
options(digits=4)
summary(bmw.garch_t)


Справочник по финансовым пакетам: http://cran.r-project.org/web/views/Finance.html

Рекомендую почитать книжку Statistics and Data Analysis for Financial Engineering (http://www.springer.com/statistics/business,+economics+%26+finance/book/978-1-4419-7786-1). Я купил и мои волосы теперь чистые и шелковистые.

а оценивать Stochastic volatility models с разными распределениями?
http://cran.r-project.org/web/packages/stochvol/index.html

ЗЫ. А что-то типа фильтра Калмана там есть?
Пакет
http://cran.r-project.org/web/packages/FKF/index.html
Статья по теме
http://www.jstatsoft.org/v39/i02/paper

[Uzanka]

Hogfather,
просто огромное спасибо!!!!!!! :jump: :jump: :jump: пошла изучать

[Hogfather]

Приобрел неплохую книжку по Data Mining в GNU R, рекомендую.

http://ecx.images-amazon.com/images/I/41S5%2Be3ahXL._BO2,204,203,200_PIsitb-sticker-arrow-click,TopRight,35,-76_AA300_SH20_OU01_.jpg (http://www.amazon.com/gp/product/1441998896/ref=oh_details_o00_s01_i00?ie=UTF8&psc=1&tag=s601000020-20)


Пакет rattle предназначен для поиска закономерностях в данных (Data Mining) с помощью регрессионных деревьев, кластерного анализа и метода опорных векторов. В книжке разбирается порядок действий от загрузки данных до интерпретации результатов.

library(rattle)
rattle()

http://aspirantura.spb.ru/forum/picture.php?albumid=111&pictureid=1364

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

UPD: Появилось сообщение на форуме поддержки. На настоящий момент решение это проблемы выглядит вот так [1] . Вложение со старой функцией удалено, чтобы не смущало. (https://groups.google.com/forum/#!topic/rattle-users/p8tlZ0Kyu2s)


The pmmlCanExport is an issue with the new release of the pmml package - it is being addressed at present - a quick fix:

> pmmlCanExport <- function(...) TRUE
> rattle()

Regards,

Graham Williams
http://togaware.com



Т.е. есть еще некоторые проблемки с продуктом,а так, в целом, здорово!

[Hogfather]

UPD: Под МакОс rattle категорически отказался работать. Проблемы с GTK+

[Hogfather]

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

Имеет тест на 45 вопросов по 4 варианта ответа в каждом. Только один вариант верный.
Тест разбит на 3 блока по 15 вопросов. За правильные ответы на вопросы из 1-го блока
начисляется 1 балл, из 2-го блока - 2 балла, из 3-го блока - 3 балла, и соответственно,
максимальный балл, который можно набрать = 90. Определить вероятность сдачи теста,
тыкая ответы наугад, при условии, что для успешной сдачи нужна набрать минимум 61.

А давайте её еще немного усложним и ко всему прочему построим графики функции вероятности и функции распределения?

[Paul Kellerman]

Задан тест на 45 вопросов по 4 варианта ответа в каждом. Только один вариант верный.
Тест разбит на 3 блока по 15 вопросов. За правильные ответы на вопросы из 1-го блока
начисляется 1 балл, из 2-го блока - 2 балла, из 3-го блока - 3 балла, и соответственно,
максимальный балл, который можно набрать = 90. Определить вероятность сдачи теста,
тыкая ответы наугад, при условии, что для успешной сдачи нужна набрать минимум 61.
А давайте её еще немного усложним и ко всему прочему построим
графики функции вероятности и функции распределения?


Пусть T - дискретная случайная величина, равная количеству набранных баллов (0...90).
Функция вероятности случайной величины T: f(t) = P(T=t)
Функция распределения случайной величины T: F(t) = P(T<=t)

Решение в математической среде Waterloo Maple 15.0.

http://4put.ru/pictures/max/772/2372608.jpg

[Hogfather]

Имеет тест на 45 вопросов по 4 варианта ответа в каждом. Только один вариант верный.
Тест разбит на 3 блока по 15 вопросов. За правильные ответы на вопросы из 1-го блока
начисляется 1 балл, из 2-го блока - 2 балла, из 3-го блока - 3 балла, и соответственно,
максимальный балл, который можно набрать = 90. Определить вероятность сдачи теста,
тыкая ответы наугад, при условии, что для успешной сдачи нужна набрать минимум 61.


Итак, вернемся к нашим баранам. Собственно, решение задачи основано на знании формулы Бернулли (http://ru.wikipedia.org/wiki/%D0%A4%D0%BE%D1%80%D0%BC%D1%83%D0%BB%D0%B0_%D0%91% D0%B5%D1%80%D0%BD%D1%83%D0%BB%D0%BB%D0%B8) и понимания биноминального распределения (http://ru.wikipedia.org/wiki/%D0%91%D0%B8%D0%BD%D0%BE%D0%BC%D0%B8%D0%B0%D0%BB%D 1%8C%D0%BD%D0%BE%D0%B5_%D1%80%D0%B0%D1%81%D0%BF%D1 %80%D0%B5%D0%B4%D0%B5%D0%BB%D0%B5%D0%BD%D0%B8%D0%B 5).

Функция вероятности для биноминального распределения в R имеет вид dbinom(x, size, prob, log = FALSE), остальные связанные функции:
pbinom(q, size, prob, lower.tail = TRUE, log.p = FALSE) - функция распределения
qbinom(p, size, prob, lower.tail = TRUE, log.p = FALSE) - квантили распределения
rbinom(n, size, prob) - генерация случайного числа, распределенного по данному закону
где, соответственно: size - число попыток, prob - вероятность удачного исхода

Поэтому задачу:
Вот к примеру есть тупо тест тупо на 45 вопросов тупо по 4 ответа
в каждом, и только один ответ верный. Для положительной оценки
нужно минимум на 31 вопрос ответить верно. Какова вероятность
сдать его, ничего не зная, тупо тыкая наугад ответы?

можно решить двумя способами, через сложение устраивающих нас исходов испытания или как разницу между 1 и вероятностью неустраивающих нас результатов, которая равна куммулятивной сумме вероятностей исходов от 0 до 30, что соответствует значению функции распределения pbinom(30, 45, 0.25)


> sum(dbinom(31:45, 45, 0.25))
[1] 7.527228e-10
> 1-pbinom(30, 45, 0.25)
[1] 7.527228e-10


В MS Excel это считается так:
=1-БИНОМ.РАСП(30;45;0,25;ИСТИНА)
и имеем результат
7,52723E-10

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

> # Инициируем переменные
> P<-0 # Переменная для подсчета куммулятивной суммы, как увидим дальше, она лишняя. Просто для явности подсчетов для решения задачи
> Z<-data.frame(X=0:90,y=rep(0,91)) # Здесь мы будем накапливать вероятности для исходо теста 0:90
> for(i in 0:15) for (j in 0:15) for(k in 0:15) { # Перебираем возможные варианты для каждой части теста
+ Z$y[i+2*j+3*k+1]<-Z$y[i+2*j+3*k+1]+dbinom(i,15,0.25)*dbinom(j,15,0.25)*dbinom(k,15,0 .25) # считаем вероятность куммулятивно для точки
+ if(i+2*j+3*k>=61) P<-P+dbinom(i,15,0.25)*dbinom(j,15,0.25)*dbinom(k,15, 0.25) # считаем для точки
+ }
> Z$cdf<-cumsum(Z$y)
> cat("Ответ на задачу составляет: ",P)
Ответ на задачу составляет: 1.73944e-08
> oldpar<-par(mfrow=c(2,1))
> plot(Z$y~Z$X,main="Функция вероятности f(x)",xlab="x",ylab="f(x)",pch=19,cex=0.8,col="blue")
> points(Z$y[62:91]~Z$X[62:91],cex=0.9,col="red",pch=19)
> text(70,0.02,expression(sum(f[i],i==61,90)==1.73944 %*% 10^{-8}))
> plot(Z$cdf~Z$X,main="Рапределения F(x)",xlab="x",ylab="f(x)",pch=19,cex=0.8,col="blue")
> abline(h=Z$cdf[61],v=60,col="darkgrey",lty=2)
> points(Z$cdf[61]~Z$X[61],cex=0.9,col="red",pch=19)
> par(oldpar)
> # Собственно, как и было сказано выше, вероятности также можно посчитать как:
> 1-Z$cdf[61] # 1-F(60); 61, а не 60, потому как нумерация в массиве начинается с 1
[1] 1.73944e-08
> sum(Z$y[62:91]) # sum(f(t),t = 61...90)
[1] 1.73944e-08


http://aspirantura.spb.ru/forum/picture.php?albumid=111&pictureid=1391

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

[_Tatyana_]

о так вот вы куда ушли с решением !
пик и есть нужный нам процент?

[Hogfather]

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

Итак, пик на верхнем графике соответствует 22 баллам со значением вероятности 6.315466e-02 или 6%. Т.е. случайно в 6% случаев можно ответить на 22 балла. Вообще, верхний график показывает вероятности набрать случайным образом баллы от 0 до 90. Поэтому, если нам надо посчитать вероятность не менее 61, мы складываем вероятности набрать 61,62, ... , 90 баллов (точки выделены красным) и получаем нужный нам ответ 1.73944e-08.

Нижний график показывает куммулятивную сумму вероятностей, т.е. в точке 60 (красная точка) он равен сумме вероятностей набрать баллы от 0 до 60, иными словами, значение в каждой точке показывает вероятность набрать не более заданного количества баллов. Для того, чтобы найти вероятность набрать более чем заданное количество баллов (в нашем случае 60, потому как 61>60 и по условиям задачи оно входит в вопрос), учитывая, что сумма всех вероятностей равна 1, имеем 1-F(60)=1.73944e-08

[_Tatyana_]

таким образом вероятнее всего будет набрать в случае случайного тыка баллов около 20, что не является проходным
а вопрос задачи был о 31 балле. распределение говорит что это три процента
правильно?

[Дмитрий В.]

а вопрос задачи был о 31 балле. распределение говорит что это три процента
правильно?
_Tatyana_, вроде бы там написано
1.73944e-08

[_Tatyana_]

_Tatyana_, вроде бы там написано
я это не видела. а график тогда врет что ли?

[Дмитрий В.]

я это не видела. а график тогда врет что ли?
_Tatyana_, почему врет? Просто вероятность методом тыка правильно ответить на 61 балл уже мала и понижается с увеличением количества правильных вопросов. Объяснение тут:
Вообще, верхний график показывает вероятности набрать случайным образом баллы от 0 до 90. Поэтому, если нам надо посчитать вероятность не менее 61, мы складываем вероятности набрать 61,62, ... , 90 баллов (точки выделены красным) и получаем нужный нам ответ 1.73944e-08.

[Hogfather]

таким образом вероятнее всего будет набрать в случае случайного тыка баллов около 20, что не является проходным
а вопрос задачи был о 31 балле. распределение говорит что это три процента
правильно?
1. Вопрос задачи несколько другой (http://aspirantura.spb.ru/forum/showpost.php?p=403606&postcount=39). О 61 балле. О 31 балле более простая задача и она разбирается вначале, там другие условия (1 балл за каждый ответ в любом случае). Графика для неё я не рисовал.
Поскольку люди путаются, я подправил свой постинг.

2. Вы совершенно правы, для 31 балла примерно 3% на графике.

[_Tatyana_]

1. Вопрос задачи несколько другой (http://aspirantura.spb.ru/forum/showpost.php?p=403606&postcount=39). О 61 балле. О 31 балле более простая задача и она разбирается вначале, там другие условия (1 балл за каждый ответ в любом случае). Графика для неё я не рисовал.
Поскольку люди путаются, я подправил свой постинг.

2. Вы совершенно правы, для 31 балла примерно 3% на графике.

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

[Hogfather]

Немного о книгах по R

http://ecx.images-amazon.com/images/I/4196EYiwYCL._BO2,204,203,200_PIsitb-sticker-arrow-click,TopRight,35,-76_AA278_PIkin4,BottomRight,-68,22_AA300_SH20_OU01_.jpg (http://www.amazon.com/100-Statistical-Tests-N-D-Lewis-ebook/dp/B00C143F14)

Купил себе сегодня справочник. Толковый, рекомендую. Поставил на Kindle и на смартфон. Всегда под рукой.
В книге рассматривается 100 статистических тестов (оглавление доступно по ссылке), также есть в начале предметный указатель по проверяемым гипотезам. У каждого теста описана нулевая гипотеза, требования к выборке, интерпретация результатов. Даются ссылки на последние работы в разных областях, в которых эти тесты упомянуты.
5 из 5

http://ecx.images-amazon.com/images/I/41kUn%2BwfNoL._SY344_PJlook-inside-v2,TopRight,1,0_SH20_BO1,204,203,200_.jpg (http://www.amazon.com/gp/product/1461427495/ref=oh_details_o01_s00_i00?ie=UTF8&psc=1&tag=s601000020-20)

Отличная книжка. Есть всё, что нужно и даже больше, её уже разок упоминал выше. 5 из 5

http://ecx.images-amazon.com/images/I/31-KyVMiwnL._BO2,204,203,200_PIsitb-sticker-arrow-click,TopRight,35,-76_AA300_SH20_OU01_.jpg (http://www.amazon.com/gp/product/0387886974/ref=oh_details_o01_s00_i01?ie=UTF8&psc=1&tag=s601000020-20)

Книжка послабее, есть некая сумбурность в подаче материала, но как основа очень неплохо. 5 из 5.

[Hogfather]

Немного об анализе текста в R. Нижеприведенный код разбирает текст (вектор mydata содержит название тем диссертаций) на слова, выделяет корни, строит частотную матрицу слов и выводит в файл картинку.

Версия с перекодировкой в UTF8. Работает под Windows и MacOS.
plot.wordcloud<-function(mydata) {

library(tm)
library(wordcloud)
library(RColorBrewer)
library(SnowballC)

strsplit_space_tokenizer <- function(x) unlist(strsplit(x, "[[:space:]]+"))

# Надо преобразовать
temp.x<-enc2utf8(tolower(paste(as.vector(mydata),collapse = " ")))
ds <- DataframeSource(data.frame(temp.x),encoding ="UTF8")

xkcd.corpus <- Corpus(ds,readerControl = list(reader = readPlain, language = "ru"))
xkcd.corpus <- tm_map(xkcd.corpus, removePunctuation)
xkcd.corpus<- tm_map(xkcd.corpus, removeWords, stopwords("russian"))
xkcd.corpus <- tm_map(xkcd.corpus, stemDocument, language = "russian")
xkcd.corpus <- tm_map(xkcd.corpus, stripWhitespace)

tdm <- TermDocumentMatrix(xkcd.corpus,control = list(tokenize=strsplit_space_tokenizer))
m <- as.matrix(tdm)
v <- sort(rowSums(m),decreasing=TRUE)
d <- data.frame(word = names(v),freq=v)
pal <- brewer.pal(9, "BuGn")
pal <- pal[-(1:2)]
pal2 <- brewer.pal(8,"Dark2")
png("wordcloud.png", width=1000,height=1000)
wordcloud(d$word,d$freq, scale=c(8,.2),min.freq=3,
max.words=Inf, random.order=FALSE, rot.per=.15, colors=pal2)
dev.off()
return(xkcd.corpus)
}


http://aspirantura.spb.ru/forum/picture.php?albumid=158&pictureid=1429

Во вложении находится архив, содержащий данные о всех диссертациях, которые находятся на сайте ВАК с 2011 года.
При желании, с этим файлом можно делать что угодно.

Загрузка в R

# Для MacOS и Linux
dissers<-read.csv2("dissers.csv",header=T,fileEncoding="cp1251",encoding="cp1251")

# Для windows
dissers<-read.csv2("dissers.csv",header=T)



А дальше строим облако для 08 специальности
plot.wordcloud(subset(dissers,substr(Nspec,1,2)=="08")$diser)

http://aspirantura.spb.ru/forum/picture.php?albumid=158&pictureid=1431

[Hogfather]

Подправил код, указал в явном виде русский в команде Corpus. Иначе сжирало буквы "ч" и "я"

Добавлено через 1 час 11 минут
А вот несколько иная задача. Нас интересует конкретное слово -- "оптимизация" .
Ну, неравнодушен я к нему. Будем отбирать специальности. Подвох в том, что синтаксический разбор тут делать нельзя, но можно воспользоваться несколькими недокументированными трюками.


x<-subset(dissers,grepl("оптимизац",tolower(diser)))$Nspec
z<-summary(x,maxsum=1000)
od<-par(mar=c(0,0,0,0))
pal2 <- brewer.pal(8,"Dark2")
wordcloud(names(z),z, scale=c(6,0.2),min.freq=2,
max.words=Inf, random.order=FALSE, rot.per=.15, colors=pal2)


http://aspirantura.spb.ru/forum/picture.php?albumid=158&pictureid=1435

Результат забавный, никогда не думал, что оптимизацией занимается медицина

[Bronepoezd]

Здравствуйте!

Подскажите, пожалуйста, из-за чего это может быть?

В файле Code.R скрипт Hogfather'а без изменений

plot.wordcloud<-function(mydata) {
library(tm)
library(wordcloud)
library(RColorBrewer)
library(SnowballC)

xkcd.corpus <- Corpus(DataframeSource(data.frame(mydata)),readerC ontrol = list(reader = readPlain, language = "ru"))
xkcd.corpus <- tm_map(xkcd.corpus, removePunctuation)
xkcd.corpus<- tm_map(xkcd.corpus, removeWords, stopwords("russian"))
xkcd.corpus <- tm_map(xkcd.corpus, stemDocument, language = "russian")
xkcd.corpus <- tm_map(xkcd.corpus, stripWhitespace)
xkcd.corpus <- tm_map(xkcd.corpus, tolower)
tdm <- TermDocumentMatrix(xkcd.corpus)
m <- as.matrix(tdm)
v <- sort(rowSums(m),decreasing=TRUE)
d <- data.frame(word = names(v),freq=v)
pal <- brewer.pal(9, "BuGn")
pal <- pal[-(1:2)]
pal2 <- brewer.pal(8,"Dark2")
png("wordcloud.png", width=1000,height=1000)
wordcloud(d$word,d$freq, scale=c(8,.2),min.freq=3,
max.words=Inf, random.order=FALSE, rot.per=.15, colors=pal2)
dev.off()
}

Результат запуска в R ниже



Введите 'demo()' для запуска демонстрационных программ, 'help()' -- для
получения справки, 'help.start()' -- для доступа к справке через браузер.
Введите 'q()', чтобы выйти из R.

> source("G:\\Program_Files\\R-3.0.2\\dissers\\code.R")
> dissers<-read.csv2("G:\\Program_Files\\R-3.0.2\\dissers\\dissers.csv",header=T,fileEncoding="cp1251",encoding="cp1251")
Предупреждение
In read.table(file = file, header = header, sep = sep, quote = quote, :
incomplete final line found by readTableHeader on 'G:\Program_Files\R-3.0.2\dissers\dissers.csv'
> plot.wordcloud(subset(dissers,substr(Nspec,1,2)=="08")$diser)
Loading required package: Rcpp
Loading required package: RColorBrewer
Ошибка в .Source(readPlain, encoding, nrow(x), FALSE, row.names(x), 0, :
vectorized sources must have positive length
> dissers<-read.csv2("G:\\Program_Files\\R-3.0.2\\dissers\\dissers.csv",header=T)
> plot.wordcloud(subset(dissers,substr(Nspec,1,2)=="08")$diser)
null device
1
> dissers<-read.csv2("G:\\Program_Files\\R-3.0.2\\dissers\\dissers.csv")
> plot.wordcloud(subset(dissers,substr(Nspec,1,2)=="08")$diser)
null device
1
>


Кодировка файла Code.R -- cp1251. И никаких картинок в папке dissers (Об этом и говорит null device 1 Так?). Работает скрипт порядка 5 минут -- так и должно быть?

[Hogfather]

Результат запуска в R ниже

В общем, проблема в следующем. Я делал под MacOS, там локаль стоит UTF-8 и все работает. Windows не даёт сменить локаль, и теряет буквы ч. Изначально у меня терял и Мак, но шаманскими заклинаниями я заставил его корректно прочитать файл. Чтобы работало под Windows загрузку делаем вот так
dissers<-read.csv2("dissers.csv",header=T)

Тогда корректно закачиваются все данные в R. А дальше глючат две вещи:
xkcd.corpus <- tm_map(xkcd.corpus, removePunctuation)
тут в регулярном выражении должны были удалиться знаки препинания, но заодно удаляются буквы "ч" и "я". Что побеждается перекодировкой в UTF-8.

И, что самое страшное, глючит вот эта команда
tdm <- TermDocumentMatrix(xkcd.corpus)


Добавлено через 11 минут
Работает скрипт порядка 5 минут -- так и должно быть?
По 08 специальности да. Там много документов.

Добавлено через 3 часа 36 минут
Итак, только в нашем цирке: TermDocumentMatrix который работает под Windows. В данном примере анализируются диссертации по двум специальностям 05.02.22 и 05.02.23, строится общее облако из всех терминов (объединение), а также анализируется, какой термин где что больше используется, а также пересечение терминов.


library(tm)
library(wordcloud)
library(RColorBrewer)
library(SnowballC)

dissers<-read.csv2("dissers.csv",header=T) # Данные читаем
strsplit_space_tokenizer <- function(x) unlist(strsplit(x, "[[:space:]]+"))


x1<-enc2utf8(tolower(paste(as.vector(subset(dissers,Ns pec=="05.02.22")$diser),collapse = " ")))
x2<-enc2utf8(tolower(paste(as.vector(subset(dissers,Ns pec=="05.02.23")$diser),collapse = " ")))


docs <- data.frame(docs = c(x1,x2),
row.names = c("05.02.22", "05.02.23"))

ds <- DataframeSource(docs,encoding ="UTF8")

ds.corpus <- Corpus(ds,readerControl = list(reader = readPlain, language = "ru"))
ds.corpus <- tm_map(ds.corpus, removePunctuation)
ds.corpus<- tm_map(ds.corpus, removeWords,stopwords("russian"))
ds.corpus <- tm_map(ds.corpus, stemDocument,language = "russian")
ds.corpus <- tm_map(ds.corpus, stripWhitespace)
#ds.corpus <- tm_map(ds.corpus, tolower)

tdm <- TermDocumentMatrix(ds.corpus,control = list(tokenize=strsplit_space_tokenizer))
m <- as.matrix(tdm)
v <- sort(rowSums(m),decreasing=TRUE)
d <- data.frame(word = names(v),freq=v)
pal <- brewer.pal(9, "BuGn")
pal <- pal[-(1:2)]
pal2 <- brewer.pal(8,"Dark2")
## Рисунок 1
png("wordcloud1.png", width=600,height=600)
wordcloud(d$word,d$freq, scale=c(6,.2),min.freq=5,
max.words=Inf, random.order=FALSE, rot.per=.15, colors=pal2)
dev.off()

## Рисунок 2
png("wordcloud2.png", width=600,height=600)
comparison.cloud(m, colors = pal2, title.size=2, max.words=500)
dev.off()

## Рисунок 3
png("wordcloud3.png", width=500,height=500)
commonality.cloud(m, colors = pal2, max.words=500)
dev.off()


http://aspirantura.spb.ru/forum/picture.php?albumid=158&pictureid=1438

http://aspirantura.spb.ru/forum/picture.php?albumid=158&pictureid=1439

http://aspirantura.spb.ru/forum/picture.php?albumid=158&pictureid=1440

Добавлено через 18 минут
# Дендрограмма
# оставляем только самы частые слова (9 дециль)
wf = rowSums(m)
m1 = m[wf>quantile(wf,probs=0.9), ]
# удаляем пустые колонки
m1 = m1[,colSums(m1)!=0]
# преобразуем в двоичный вид
m1[m1 > 1] = 1

# матрица двоичных дистанций
m1dist = dist(m1, method="binary")

# кластер с использованием объединения по методу Варда
clus1 = hclust(m1dist, method="ward")

# дендрограмка
plot(clus1, cex=0.7)

http://aspirantura.spb.ru/forum/picture.php?albumid=158&pictureid=1441

Добавлено через 32 минуты
Своды документов можно объединять командой c()

# Расширяем нашу коллекцию
x3<-enc2utf8(tolower(paste(as.vector(subset(dissers,Ns pec=="08.00.05")$diser),collapse = " ")))

docs <- data.frame(docs = x3,
row.names = "08.00.05")

ds1 <- DataframeSource(docs,encoding ="UTF8")

ds1.corpus <- Corpus(ds1,readerControl = list(reader = readPlain, language = "ru"))
ds1.corpus <- tm_map(ds1.corpus, removePunctuation)
ds1.corpus<- tm_map(ds1.corpus, removeWords,stopwords("russian"))
ds1.corpus <- tm_map(ds1.corpus, stemDocument,language = "russian")
ds1.corpus <- tm_map(ds1.corpus, stripWhitespace)

tdm <- TermDocumentMatrix(c(ds.corpus,ds1.corpus),control = list(tokenize=strsplit_space_tokenizer))

## Рисунок 4
m <- as.matrix(tdm)
png("wordcloud4.png", width=1000,height=1000)
comparison.cloud(m, colors = pal2, title.size=2, max.words=500)
dev.off()

http://aspirantura.spb.ru/forum/picture.php?albumid=158&pictureid=1442

Добавлено через 6 часов 56 минут
Поскольку просто так потратить кучу времени на изучение пакета было бы обидно, пишу статью в нашу вузовскую "Мурзилку" по мотивам анализа данных. Коню понятно, что "играть" на форуме -- это одно, а более-менее серьезное исследование -- это другое, то возник вопрос, что делать с дубликатами защит. А они есть и искажают картину качественно и количественно. Так вот, в R все давно придумано до нас и не надо подключать никаких дополнительных пакетов. Есть чудесная команда duplicated, которая возвращает логический вектор.

Вызов duplicated(dissers[,c(1,3,5)], fromLast = TRUE) проверяет наличие строк, в которых дублированы специальность, вид защиты и фамилия. Логика следующая: дубликаты, в большинстве случаев, возникают когда переносят защиту или когда есть ошибка в названии АР. Обратите внимание: убрано 1185 дублирующихся записей. Теперь анализ проводить гораздо корректнее. Напоминаю, что команда dim выводит размер матрицы или таблицы данных (число строк, число столбцов). По сути, нужна только команда dissers<-dissers[!duplicated(dissers[,c(1,3,5)], fromLast = TRUE),]


> dim(dissers)
[1] 43522 6
> dissers<-dissers[!duplicated(dissers[,c(1,3,5)], fromLast = TRUE),]
> dim(dissers)
[1] 42337 6


Добавлено через 2 часа 36 минут
Трансформация данных и построение графика


library(reshape2)
md<-data.frame(spec=substr(dissers$Nspec,1,2),time=pas te(substr(dissers$date,7,10),substr(dissers$date,4 ,5),sep="."),cnt=1)
res<-dcast(md,time~spec,sum,value.var="cnt")
row.names(res)<-res$time
res<-as.matrix(t(res[,-1]))
oldpar<-par(mai=c(1.36,1.09,1.09,1.56))
barplot(as.matrix(res)[,-1],col=rainbow(18),legend.text=rownames(res),las=2,a rgs.legend=list(x=35),main="Динамика защит")
par(oldpar)


http://aspirantura.spb.ru/forum/picture.php?albumid=158&pictureid=1444

Добавлено через 11 часов 17 минут
Продолжаем издеваться над теми же данными. А вот пример диаграммы Венна. Очень меня интересуют слова-заглушки: инновации, оптимизация, механизм и т.п. Сказано -- сделано. Строим табличку совпаданий, а затем строим диаграмму.

library(gplots)
dis.df<-dissers # Тут может быть любая выборка
Test.df<-data.frame("Качество"=grepl("качеств",tolower(dis.df$diser)),
"Инновации"=grepl("инновац",tolower(dis.df$diser)),
"Оптимизация"=grepl("оптимизац",tolower(dis.df$diser)),
"Эффективность"=grepl("эффективн",tolower(dis.df$diser)),
"Механизм"=grepl("механизм",tolower(dis.df$diser)))
rm(dis.df)
venn(Test.df)

http://aspirantura.spb.ru/forum/picture.php?albumid=158&pictureid=1445

Радуемся "оптимизации эффективности качества" и т.п.

Добавлено через 13 часов 28 минут
Следующий пример визуализации данных несколько из другой области. Чисто экономические заморочки -- бридж диаграмма, также называемая waterfall chart (http://en.wikipedia.org/wiki/Waterfall_chart). Так вот, в R есть и она. Диаграмма может быть полезна визуализации экономического эффекта с разбивкой по факторам. Можно строить в классическом виде, а можно как в журнале The Economist. Собственно, не растекаясь мысью по древу, код и картинки.


> library(waterfall)
> data(rasiel) # Пример данных, имеющихся в пакете
> rasiel
label value subtotal
1 Net Sales 150 EBIT
2 Expenses -170 EBIT
3 Interest 18 Net Income
4 Gains 10 Net Income
5 Taxes -2 Net Income
> waterfallchart(value~label, data=rasiel, groups=subtotal, main="P&L")
> asTheEconomist(waterfallchart(value~label, data=rasiel, groups=subtotal, main="P&L"))

http://aspirantura.spb.ru/forum/picture.php?albumid=158&pictureid=1448

http://aspirantura.spb.ru/forum/picture.php?albumid=158&pictureid=1449

Также можно воспользоваться пакетом ggplot2 и нарисовать графики с использованием его инструментов. Подробности можно посмотреть вот тут (http://learnr.wordpress.com/2010/05/10/ggplot2-waterfall-charts/).

Добавлено через 5 часов 47 минут
Еще немного визуализации нечисловых данных. Анализируем пропорции авторефератов на сайте ВАК по экономике.

library(reshape2)
library(vcd)
tmp<-subset(dissers,substr(Nspec,1,2)=="08")[,2:4]
tmp$date<-substr(tmp$date,7,10)
for(i in 1:3) tmp[,i]=factor(tmp[,i])
mosaic(~TypeOfDisser+Nspec+date, data=tmp,expected=~TypeOfDisser:Nspec + TypeOfDisser:date + Nspec:date,legend=FALSE, gp=shading_binary,pop=FALSE,
labeling_args=list(rot_labels=c(right=0),gp_labels =gpar(fontsize=8)))



http://aspirantura.spb.ru/forum/picture.php?albumid=158&pictureid=1450

Возможности пакета vsd этим не исчепывается. Подробнее в документации к пакету (http://cran.r-project.org/web/packages/vcdExtra/vignettes/vcd-tutorial.pdf).

[Bronepoezd]

Причина: Окончательная версия, работающая на всех машинах

Благодарю! Теперь действительно работает под windows. Полезный скрипт.

[Дмитрий В.]

Благодарю! Теперь действительно работает под windows. Полезный скрипт.
Присоединяюсь, большое спасибо!

[Hogfather]

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

[Дмитрий В.]

Hogfather, интересно должно получиться, я думаю.

[Hogfather]

Еще немного о работе с текстами. Берем послания к федеральному собранию за последние 4 года и рисуем облако. Данный пример показывает, как удобно грузятся данные из текста в UTF-8 в заданном каталоге.Теоретически, R также читает doc и pdf, но через внешние конвертеры в txt.


library(tm)
library(wordcloud)
library(RColorBrewer)
library(SnowballC)

strsplit_space_tokenizer <- function(x) unlist(strsplit(x, "[[:space:]]+"))
ds<-DirSource(directory=".\\texts\\",encoding="UTF-8",pattern="*.txt")
ds.corpus <- Corpus(ds,readerControl = list(reader = readPlain, language = "ru"))
ds.corpus <- tm_map(ds.corpus, removePunctuation)
ds.corpus<- tm_map(ds.corpus, removeWords,stopwords("russian"))
ds.corpus <- tm_map(ds.corpus, stemDocument,language = "russian")
ds.corpus <- tm_map(ds.corpus, stripWhitespace)
ds.corpus <- tm_map(ds.corpus, tolower)

tdm <- TermDocumentMatrix(ds.corpus,control = list(tokenize=strsplit_space_tokenizer))
m <- as.matrix(tdm)
v <- sort(rowSums(m),decreasing=TRUE)
d <- data.frame(word = names(v),freq=v)
pal1 <- rainbow(20)
pal2 <- brewer.pal(8,"Dark2")


# Смотрим слова, которые встречаются более 200 раз
png("poslan.png", width=1000,height=1000)
comparison.cloud(m, colors = pal2, title.size=1.5,scale=c(6,.2), max.words=400)
dev.off()


http://aspirantura.spb.ru/forum/picture.php?albumid=158&pictureid=1451

Исходные данные во вложении, так что экспериментируйте на здоровье.

Самые распространенные слова в посланиях.

> findFreqTerms(tdm, 50)
[1] "важн" "возможн" "вопрос" "год"
[5] "государств" "государствен" "дет" "должн"
[9] "друг" "ещё" "задач" "котор"
[13] "люд" "нам" "наш" "необходим"
[17] "нов" "нужн" "обществ" "прав"
[21] "правительств" "программ" "работ" "работа"
[25] "развит" "регион" "решен" "росс"
[29] "российск" "сам" "сво" "сдела"
[33] "систем" "современ" "созда" "социальн"
[37] "стран" "сфер" "счита" "так"
[41] "такж" "уважа" "цел" "числ"
[45] "экономик" "экономическ" "это"


Интересный инструмент -- поиск ассоциаций с заданным уровнем корреляции.

> findAssocs(tdm,"стратегическ",0.98)
источник проблем флот инструмент одн результат самоуправлен
1.00 1.00 1.00 0.99 0.99 0.99 0.99
страхов послан привест разрешен район
0.99 0.98 0.98 0.98 0.98
> findAssocs(tdm,"дет",0.98)
воспитан окружа семейн суд увеличен дом размер сем 2008
1.00 1.00 1.00 1.00 1.00 0.99 0.99 0.99 0.98
детств зал насил партнёрств применя умн
0.98 0.98 0.98 0.98 0.98 0.98


NB! Обсуждение картинки вот тут (http://aspirantura.spb.ru/forum/showthread.php?t=10607&page=39).

[Дмитрий В.]

Hogfather,
зал
- можно ли посмотреть, какого слова (каких слов) корень? Т.е. это просто слово "зал", или, может, ошибочно определилось...

[Hogfather]

Дмитрий В., надо при обработке корпуса не выделять корни. Других способов нет.

[Alena1]

Здравствуйте, помогите пожалуйста написать Garch модель в R

[Hogfather]

Здравствуйте, помогите пожалуйста написать Garch модель в R
См. http://www.aspirantura.spb.ru/forum/showpost.php?p=396715&postcount=35

Пока не вижу, в чем надо помочь. Вопроса как такового нет. А если "помогите пожалуйста"="напишите за меня", то предупреждаю, что я этим заниматься не буду.

[Hogfather]

Нам пишут

В качестве стат программы я привык работать с SPSS и в принципе она меня устраивала с точки зрения тех параметров, которые она дает возможность анализировать
Однако, видимо, R дает больше возможностей (а так ли это? в чем принципиальная разница подходов для гуманитария?)
С чего начать и в какой программной среде работать?

1. Больше ли возможности R?
R даёт больше возможностей хотя бы тем, что для неё существует бесчисленное множество пакетов, покрывающих практически любую извращенную фантазию исследователя. Также определенное удобство дает то, что это скриптовый язык и у исследователя имеется возможность делать типовые действия и потом их повторять, также удобно обмениваться данными и кодом с другими исследователями. Например, с помощью Github или Rpubs (http://rpubs.com/)

2. В чем разница подходов?
Принципиальная разница подходов в том, что это не ориентированный на кляцкание мышкой инструмент. Надо учить "волшебные слова", т.е. команды и вести диалог с компьютером, объясняя ему на языке R, что вам от него надо, поэтому потребуется некоторое время на то, чтобы перестроиться.

3. С чего лучше начать и в какой среде работать?
Начать лучше с начала, самое разумное, с бесплатных курсов (https://www.coursera.org/specialization/jhudatascience/1?utm_medium=listingPage).В первом курсе как раз рассказывают про инструменты исследователя. Также есть блог, который посвящен именно R и там есть список учебников http://r-analytics.blogspot.ru/p/blog-page_20.html
Ну, и само собой, необходимо установить R и RStudio (http://www.rstudio.com/products/rstudio/). На мой взгляд это самое удачное сочетание. Тут и работать.

Если перед вами стоят задачи классического DataMining (кластерный и регрессионный анализ, машинное обучение), то стоит обратить внимание на пакет rattle (http://rattle.togaware.com/RattleBrochure.pdf). Его я рекомендую студентам, которым нужно решать такие задачи тыкая мышкой в кнопки. Что хорошо, в результате создается скрипт-протокол, который можно будет использовать потом уже не тыкая мышкой по экрану.

Если кратко, то вот так.

[Hogfather]

После обновления RStudio (http://www.rstudio.com/) опять началась ерунда с путями, содержащими русские буквы под Windows.

Чтобы не забыть пишу, как я с этим борюсь.
1. Добавил переменную окружения в меню:
Панель управления-Система-Дополнительно-Переменные среды
R_LIBS_USER со значением C:\RLIB (как вариант)

Перенес туда все пользовательские библиотеки из папки C:\Users\Хогфазер\Documents\R\win-library\3.1

2. Тем путям к папкам с текущим R проектом, которые используют русские буквы, поставил в соответствие виртуальный диск, для чего написал простенький командный файл runme.cmd, который разметил на рабочем столе:


subst x: "C:\Users\Хогфазер\Dropbox\Статьи\201 4\КГАМ\Совсем уж ерунда"
pause
subst x: /D


Напоминаю, что кодировка командного файла должна быть DOS (CP866).

[Hogfather]

О R под андроид.

Есть три варианта работать с R под Android.
Самый простой: поставить приложение R Console Free (https://play.google.com/store/apps/details?id=com.appsopensource.R)
Из недостатков: проблемы с графикой. Т.е. графики не рисует. Но работает.

Второй вариант: развернуть дистрибутив Debian и под ним поставить R
Описание вот тут (http://www.r-ohjelmointi.org/?p=1434)
Теоретически, можно заставить и графики рисовать (в комментариях есть рецепт), но у меня не настолько хорошее зрение, чтобы на телефоне что-то рассмотреть.

Третий вариант. http://www.r-fiddle.org/
Все делать в облаках. C русским языком проблемы, зато рисует графики.

[Hogfather]

"Давненько не держал в руках я шашек".

Ну, что за Новый год, да без ёлочки (http://www.r-bloggers.com/geom_christmas_tree-a-new-geom-for-ggplot2-v2-0/).

Опять же, что интересно людям. Людям интересен козёл и тигр курс доллара, да цена на нефть.

Итак, в сегодняшнем коде показываем.
1. Как получить данные с Яндекс (недокументировано)
2. Как слить данные в одну таблицу
3. Как воспользоваться преимуществами новой версии ggplot и использовать возможность определить пользовательские символы для отображения графиков.
4. Как удовлетворить строгих редакторов журналов, устанавливающих размер шрифта и размер картинки.

Итак, строим график за 2015 год зависимости курса ЦБ от цены на нефть.


library(ggplot2)
require(XML)

# Загружаем цену нефти
download.file("https://news.yandex.ru/quotes/graph_1006.xml","graph_1006.xml")
data <- xmlParse("graph_1006.xml")
xml_data <- xmlToList(data)
graphData1<-data.frame(Dat=format(as.Date(as.numeric(unlist(st rsplit(xml_data$x,";")))/3600/24, origin = "1970-01-01"),"%Y-%m-%d"),
Brent=as.numeric(unlist(strsplit(xml_data$y$text,";"))))

# Загружаем курс USD в рублях
download.file("https://news.yandex.ru/quotes/graph_1.xml","graph_1.xml")
data <- xmlParse("graph_1.xml")
xml_data <- xmlToList(data)
graphData2<-data.frame(Dat=format(as.Date(as.numeric(unlist(st rsplit(xml_data$x,";")))/3600/24, origin = "1970-01-01"),"%Y-%m-%d"),
USD=as.numeric(unlist(strsplit(xml_data$y$text,";"))))

# Объединяем
graphData<-merge(graphData1,graphData2,by="Dat")
graphData$year=as.factor(format(as.Date(graphData$ Dat),"%Y"))
graphData$month=as.factor(format(as.Date(graphData $Dat),"%m"))
graphData$wd=as.factor(format(as.Date(graphData$Da t),"%u"))

# "В лесу родилась ёлочка..."
# См. http://www.r-bloggers.com/geom_christmas_tree-a-new-geom-for-ggplot2-v2-0/
# Нижеследующий код нужен для построения елочек

GeomChristmasTree <- ggproto("GeomChristmasTree", Geom,
required_aes = c("x", "y"),
default_aes = aes(shape = 19, colour = "black",
fill = "green4", size = 3,
linetype = 1, alpha = 1,
fontsize = 1),
draw_key = draw_key_polygon,

draw_panel = function(data, panel_scales, coord) {
coords <- coord$transform(data, panel_scales)

# each tree has 4*branches + 3 points
if (length(coords$size) == 1) {
tsize <- rep(pmax(1, round(coords$size)), length(coords$x))
theight <- rep(pmax(0, round(coords$size)), length(coords$x))
} else {
tsize <- pmax(1, round(coords$size))
theight <- pmax(0, coords$size)
}

# scale factors
r01x <- diff(range(coords$x))/100
r01y <- diff(range(coords$y))/100

# coords
longx <- unlist(lapply(seq_along(coords$x), function(i) {
if (tsize[i] == 1) {
dx <- -c(0.3, 0.3, 1.2, 0, -1.2, -0.3, -0.3)
} else {
dx <- -c(0.3, 0.3, rep(c(1.2,0.3), tsize[i]-1), 1.2, 0, -1.2, rep(c(-0.3,-1.2), tsize[i]-1), -0.3, -0.3)
}
r01x*dx + coords$x[i]
}))
longy <- unlist(lapply(seq_along(coords$y), function(i) {
if (tsize[i] == 1) {
dy <- c(-0.5, 0, 0, theight[i], 0, 0, -0.5)
} else {
dy <- c(-0.5, 0, 0, rep(1:(tsize[i]-1), each=2), theight[i], rep((tsize[i]-1):1, each=2), 0, 0, -0.5)
}
r01y*dy + coords$y[i]
}))
longid <- unlist(sapply(seq_along(coords$y), function(i) {
rep(i, each=4*tsize[i]+3)
}))

grid::polygonGrob(
longx,
longy,
id = longid,
gp = grid::gpar(col = coords[,"colour"],
fill = coords[,"fill"],
fontsize = 10)
)
}
)

geom_christmas_tree <- function(mapping = NULL, data = NULL, stat = "identity",
position = "identity", na.rm = FALSE, show.legend = NA,
inherit.aes = TRUE, ...) {
layer(
geom = GeomChristmasTree, mapping = mapping, data = data, stat = stat,
position = position, show.legend = show.legend, inherit.aes = inherit.aes,
params = list(na.rm = na.rm, ...)
)
}
# "В лесу она росла..."

# Строим график зависимости курса от цены на нефть


ggplot(subset(graphData,year=="2015"), aes(x=Brent, y=USD)) +
stat_smooth()+
geom_christmas_tree(aes(fill=month,size=wd)) +
labs(x="Цена на Нефть Brent (ICE.Brent), USD/баррель",y="Курс USD ЦБ РФ, руб.")+
scale_fill_discrete(name = "Месяц")+
scale_size_discrete(guide=FALSE)+
theme_bw()+theme(text=element_text(size=14))

# Сохраняем картинку. Размер в дюймах.
ggsave("Pic1.png",width=6,height=6,dpi=300)


Результат

http://www.aspirantura.spb.ru/forum/picture.php?albumid=158&pictureid=1786

Цвет ёлочки определяется месяцем, а размер -- днем недели. Обижать художника прошу вот здесь (http://www.aspirantura.spb.ru/forum/showthread.php?p=560638#post560638).


ggplot(subset(graphData,year=="2015"), aes(x=Brent, y=USD)) +
stat_density_2d(aes(color=month))+
geom_christmas_tree(aes(fill=month,size=wd)) +
labs(x="Цена на Нефть Brent (ICE.Brent), USD/баррель",y="Курс USD ЦБ РФ, руб.")+
scale_fill_discrete(name = "Месяц")+
scale_size_discrete(guide=FALSE)+
scale_color_discrete(guide=FALSE)+
theme_bw()+theme(text=element_text(size=14))
ggsave("Pic1.png",width=6,height=6,dpi=300)
http://www.aspirantura.spb.ru/forum/picture.php?albumid=158&pictureid=1788

[Hogfather]

1. Регулярные выражения в R (http://www.johndcook.com/blog/r_language_regex/) (POSIX 1003.2 и Perl)

2. Применение эмоджи на ggplot (http://www.r-bloggers.com/emojis-in-ggplot-graphics/). Что-то в таком стиле:
http://revolution-computing.typepad.com/.a/6a010534b1db25970b01bb08974934970d-800wi

3. Интересная статья про получение дисперсии и матожидания из медианы и размаха (http://bmcmedresmethodol.biomedcentral.com/articles/10.1186/1471-2288-5-13). Идея в следующем.

Матожидание:
https://s0.wp.com/latex.php?latex=%5Cbar%7Bm%7D+%3D+%5Cfrac%7Ba%2B2+ m%2Bb%7D%7B4%7D+%2B%5Cfrac%7Ba-2+m%2Bb%7D%7B4+n%7D&bg=ffffff&%23038;fg=333333&%23038;s=0

Дисперсия:
https://s0.wp.com/latex.php?latex=%5Cfrac%7B1%7D%7Bn-1%7D+%5CBig%28a%5E2%2Bm%5E2%2Bb%5E2%2B%5Cfrac%7Bn-3%7D%7B2%7D+%5Cfrac%7B%28a%2Bm%29%5E2%2B%28b%2Bm%2 9%5E2%7D%7B4%7D-n+%5Cbar%7Bm%7D+%5CBig%29&bg=ffffff&%23038;fg=333333&%23038;s=0

f<-function(a,m,b,n)
{
mn<-(a+2*m+b)/4+(a-2*m+b)/(4*n)
s=sqrt((a*a+m*m+b*b+(n-3)*((a+m)^2+(m+b)^2)/8-n*mn*mn)/(n-1))
c(mn,s)
}

4. Симуляция Цепей Маркова с дискретным временем в R (http://www.r-bloggers.com/a-discrete-time-markov-chain-dtmc-sir-model-in-r/) (DTMC).

5. Для R-Studio разработаны шпаргалки (https://www.rstudio.com/resources/cheatsheets/)
https://www.rstudio.com/wp-content/uploads/2016/01/rstudio-IDE-cheatsheet.png

6. Задачи связанные с методом анализа иерерхии (МАИ) теперь легко решать в R (http://www.r-bloggers.com/analytic-hierarchy-process-ahp-with-the-ahp-package/) с использованием пакета ahp (https://github.com/gluc/ahp). Также в этом пакете можно использовать функцию предпочтения (http://www.r-bloggers.com/analytic-hierarchy-process-ahp-using-preferencefunction-in-ahp/).

7. Избегаем наложение тектовой информации на графиках с помощью пакета ggrepel (https://cran.r-project.org/web/packages/ggrepel/vignettes/ggrepel.html)

https://pbs.twimg.com/media/CX-6IhkWcAUIWYW.png

8. Статья, как обнаружить гетероскедатичность и убрать её из моделей (http://www.r-bloggers.com/how-to-detect-heteroscedasticity-and-rectify-it/).

9. Обновился пакет для построения тернарных графиков ggtern (http://www.r-bloggers.com/ggtern-2-0-now-available/). Тернарные графики полезны, когда исследуюшь взаимосвязь трех переменных. Выглядят они, например, вот так:

http://www.ggtern.com/wp-content/uploads/2016/01/ggtern_basis.png

[Hogfather]

10. Revolution R переименован в Microsoft R и доступен бесплатно для разработчиков и студентов (http://habrahabr.ru/post/275113/)

[Hogfather]

Столкнулся с проблемой: установлен R, выхода в Интернет в вузе на учебных компьютерах нет, но нужно устанавливать дополнительные пакеты для R. Поможет пакет miniCRAN, который создает уменьшенную копию репозитория со всеми зависящими пакетами. Короткий код, чисто для примера: создание репозитория для самых известных пакетов и последующая установка пакета из полученного репозитория (под windows).

library("miniCRAN")
tags <- c("ggplot2", "data.table", "plyr", "knitr", "shiny", "xts", "lattice")
pkgList <- pkgDep(tags, suggests = TRUE, enhances=FALSE, availPkgs = cranJuly2014)
dir.create(pth <- file.path("miniCRAN"))
makeRepo(pkgList, path=pth, type="source")
updateRepoIndex(path, type = "source")
list.files(pth, recursive = TRUE)
makeRepo(pkgList, path=pth, type="win.binary")
updateRepoIndex(path, type = "win.binary")
pkgAvail(repos=pth)
install.packages("ggplot2", repos ="file:///test/miniCRAN")