Множественная регрессия. Краткий справочник.
Попробуем резюмировать.
Для подгонки модели МНК, используем код Код:
> fit <- lm(lnZ ~ x + y, data=myData) Код:
> coefficients(fit) # коэффициенты модели Код:
> fit1 <- lm(lnZ ~ y, data=myData) Автоматическая подгонка модели посредством исключения незначимых членов. Код:
> library(MASS) |
Краткий ликбез по моделям, или как читать то, что мы получаем в результате подгонки.
Хотя меня особо не спрашивают, но имея опыт работы со студентами постараюсь кратко и без лишнего занудства изложить, что же Пытливый Ум может увидеть в сборище букв на языке эвентуального противника, когда строит модель.
Рассмотрим пример, который был приведен выше. Код:
> summary(fit) 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) Второе, на что мы смотрим, это на остатки модели (Residuals). Код:
Residuals: Третье, на что мы смотрим -- коэффициенты. Код:
Coefficients: Дальше, 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) А теперь смотрите на график 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) Код:
> library(MASS) 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)]) Делаем преобразование и строим новую модель. Код:
> lm1<-lm(I(y^lambda)~x,data=test) Код:
> plot(lm1,which=1) Воот! Теперь совсем другое дело. Следует отметить особые случаи в преобразовании Бокса-Кокса. Если λ=1, то преобразование не требуется. Если λ=0, то следует использовать логарифмическое преобразование log(y). Пример. Код:
> test<-data.frame(x<-1:20) Скрипт для графика
Код:
> par(mfrow=c(1,2)) http://aspirantura.spb.ru/forum/pict...pictureid=1020 Добавлено через 18 часов 12 минут Кореллограммы В рассмотренной задаче (#10) рисовалась корреляционная матрица в виде соответствующих диаграмм рассеяния. Также можно получить корреляционную матрицу в виде коэффициентов корреляции в текстовом виде. Код:
> # Инициализация данных Код:
> symnum(cor(myData[,c("x","y","z")])) Код:
> symnum(cor(myData[,c("x","y","lnZ")])) Однако, есть и другой способ сделать это, да и к тому же красиво, нарисовав коррелограмму. Такой рисунок несказанно обогатит статью и поразит всех читателей. Для этого нужна установленная библиотека corrgram. Примеры вызовов функции приведены ниже. Существует 4 возможных отображения корреляции: panel.pts -- диаграмма разброса panel.pie -- круговая диаграмма, меняет цвет и заполненность в зависимость от коэффициента корреляции (отрицательная -- красная, положительная -- синяя). panel.ellipse -- рисует красоту, как показана на рисунке 3. В общем, границы и аппроксимацию. panel.shade -- закрашивает прямоугольник в зависимости от коэффициента корреляции. Их присваивают, соответственно lower.panel или upper.panel, в зависимости от того, где нужен рисунок: сверху или снизу. Также можно выводить название переменных и минимум/максимум значений, как это показано ниже. Код:
> Код:
> corrgram(myData[,c("x","y","lnZ")], lower.panel=panel.pts, Код:
> corrgram(log(myData[,c("x","y","z")]), order=TRUE, lower.panel=panel.shade, Понятно, что на трех переменных это не сильно показательно, но если их будет десяток, то картинка будет гораздо симпатичнее. |
Регресионные деревья
Чем хорош R, что в нём можно много всякого забавного сделать. например, захотелось поиграть в регрессионные деревья -- никто слова плохого не скажет.
Рассмотрим ещё раз нашу задачу. Код:
Результат, как видите следующий. Если y<912.5, то ln(z) равен 1.319, а если y>=912.5, то 3.499. Смысла для данной задачи никакого, разве что лишний раз убедились, что Икс у нас не при делах. Как это должно работать
В данном случае я немного модифицировал штатный пример из документации В качестве исходных данных берется таблица данных по кифозу у детей из 81 строк and 4 столбцов, описывающая детей у которых были корректирующие операции на позвоночнике Таблица содержит следующие колонки: Kyphosis - качественные данные (factor) с вариантами "absent" и "present", показывающая был ли кифоз после операции. Age- возраст в месяцах Number - число затронутых позвонков Start - номер верхнего позвонка, на котором происходила операция Код:
> # Первые строки таблицы |
Hogfather,
мне бывает нужен МНК на очень подробных сетках (1000\times 1000, например). Потянет это оболочка? (без ограничения производительности железа) |
Цитата:
Обратите внимание, я специально выводил время. 13 секунд на модель. Дольше таблицу инициализировал... Процесcор старый, Intel Core Duo T7250 2 ГГц. Озу 3 Гб Скрытый текст
Код:
> # Инициализируем матрицу 1000х1000 для x1..x1000, с матожидаением i |
Hogfather, спасибо!
Синтаксис не точно понимаю, но я понял, что R справился. Как он отдает результаты? |
Цитата:
|
Цитата:
Но может понадобиться что-нибудь сугубо универсальное, типа текстового файла, для демонстрации результатов через всякие Excel'и с Maple'ами. (ведь R, как я понял, -- GNU, а те далеко нет) |
Цитата:
Дублирую табличку из темы про GRETL. [table 1 5 1]Откуда берем данные|Формат обмена|Куда перегоняем|Примечание Excel|xls, xlsx или csv|Gretl| Читается прямо в формате Excel Excel|csv|GNU R |Можно и напрямую, но с доп. библиотекой, а также через clipboard (см. read.table("clipboard",header=T) ) 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" [/table] |
Hogfather, merci за консультацию.
|
Текущее время: 04:04. Часовой пояс GMT +3. |
Powered by vBulletin® Version 3.8.8
Copyright ©2000 - 2024, vBulletin Solutions, Inc. Перевод: zCarot
© 2001—2024, «Аспирантура. Портал аспирантов»