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

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

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

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

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

[table 1 3 0]3^http://aspirantura.spb.ru/forum/pict...&pictureid=979 |

Внимание![BR]Автор темы совместно с Администрацией портала просит писать в эту тему только относящееся к GNU R. Благодарности, разговоры о погоде будут безжалостно удаляться. Все "чмоки", пожалуйста, во флейме. Вопросы по R прошу задавать в личных сообщениях, чтобы тема была удобна для восприятия. Надеюсь на понимание.

[/table]

Итак, поскольку я показал как пользоваться R для решения задачи подбора распределения, а добро не должно остаться безнаказанным, в личные сообщения посыпались вопросы. Чтобы два раза не вставать, ответы публикую сюда.

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

1. Подгонка распределения. fitdistr из MASS
2. Подгонка распределения. fitdist из fitdistrplus. Лучше!
3. Нелинейный метод наименьших квадратов
4. Гамма-распределение. Вывод графика.
5. Диаграмма "Ящик с усами" для нескольких векторов
6. Упрощение себе жизни с помощью написания функции в R
7. Подгонка МНК

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

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

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

A: Почитать статьи вот тут

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(13,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/pict...&pictureid=973

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

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

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

A: Воспользуйтесь командой windows()

Код:

>windows()

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

Код:

>?windows

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

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

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

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

Вывод в файл

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.

Получение данных

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/pict...&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) Получить из базы данных
Подробнее см. тут

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

Цитата:

Сообщение от SergeyMastitsky (Сообщение 291870)

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

Подтверждаю. Блог полезный, кошерный, рекомендуется к посещению и вдумчивому прочтению.

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

Цитата:

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

Есть ли в R
1)Фурье-анализ

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

fft

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

Цитата:

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

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

Естественно есть, причем разными методами. Я уже рассматривал нелинейный МНК. Трехмерные графики тоже есть. Вот тут есть разбор пакетов.

Классический эффектный пример нормального распределения в двух измерениях [1].

Код:

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~sqrt(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/pict...&pictureid=991

Допиливаем пакет fitdistrplus

Рассмотренная ранее диаграмма Каллена и Фрея (Cullen and Frey graph) всем хороша, но она, во-первых, на языке эвентуального противника, во-вторых, нет возможности менять заголовок. Такое безобразие долго продолжаться не могло и, в процессе работы над статьёй, пришлось подправить функцию и русифицировать её вывод. Для этого идём на CRAN, скачиваем "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/pict...&pictureid=992

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

Обживаем RStudio

Прочитав про RStudio решил её обжить и делюсь решением первых проблем.

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

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

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

Множественная линейная регрессия

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

Цитата:

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

Задача: получить зависимость 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/pict...&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/pict...pictureid=1001
Поясняю, мы от модель вида http://www.texify.com/img/%5Cnormals...%5Cepsilon.gif плавно перешли к http://www.texify.com/img/%5Cnormals...ilon%29%7D.gif. В данном случае эпсилон у нас изображает остатки модели, которые, в идеальном варианте, должны иметь нормальное распределение.

Третий вариант мы можем получить логарифмированием всего и вся. В таком случае имеем http://www.texify.com/img/%5Cnormals...scr%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/pict...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/%5Cnormals...scr%7BE%7D.gif
Едем дальше. Впереди дисперсионный анализ.

Код:

> 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/pict...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/pict...pictureid=1004

Множественная регрессия. Краткий справочник.

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

Код:

> 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) [COLOR="rgb(46, 139, 87)"]# вычисленные значения[/COLOR]

        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) [COLOR="rgb(46, 139, 87)"]# дисперсионный анализ[/COLOR]

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

Прочую красоту смотреть вот здесь. А подробно проверка регрессионной модели описана вот тут.

Краткий ликбез по моделям, или как читать то, что мы получаем в результате подгонки.

Хотя меня особо не спрашивают, но имея опыт работы со студентами постараюсь кратко и без лишнего занудства изложить, что же Пытливый Ум может увидеть в сборище букв на языке эвентуального противника, когда строит модель.

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

Код:

> 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

В данном случае нулевая гипотеза выглядит так: H₀: β₁=β₂=...β_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/pict...pictureid=1016

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

Код:

Residuals:

     Min       1Q   Median       3Q      Max 

-0.61700 -0.30351 -0.03218  0.19850  0.81873

Алгоритм МНК подразумевает, что матожидание остатков должно быть равно 0 и мы надеемся, что остатки имеют нормальное распределение. Напоминаю, что у нормального распределения матожидание и медиана совпадают. В данном случае -- не очень. Поскольку мы видим отрицательную медиану (-0.03218), то распределение остатков незначительно смещено влево. Квартили, по идее, должны быть симметричны относительно 0, здесь же они значительно смещены влево. Максимум и минимум позволяют оценить выбросы (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

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

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

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

Т.е. модель у нас неадекватна, не смотря на милый Эр Квадрат. Такие дела.
Для сравнения посмотрите график на предыдущей странице

Добавлено через 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/pict...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/pict...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/pict...pictureid=1020

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

Кореллограммы

В рассмотренной задаче (#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/pict...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/pict...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/pict...pictureid=1023

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

Регресионные деревья

Чем хорош R, что в нём можно много всякого забавного сделать. например, захотелось поиграть в регрессионные деревья -- никто слова плохого не скажет.

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

Код:



> 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/pict...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/pict...pictureid=1025

Hogfather,

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

Цитата:

Сообщение от Olafson (Сообщение 297627)

мне бывает нужен МНК на очень подробных сетках (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