Рассчет множественной регрессии для серии файлов с данными

Вопросы по статистическому пакету R. Не обязательно гео.
Ответить
womak
Участник
Сообщения: 73
Зарегистрирован: 13 окт 2006, 06:19
Репутация: 0

Рассчет множественной регрессии для серии файлов с данными

Сообщение womak » 12 фев 2016, 18:09

Уважаемые коллеги!
Передо мной стоит задача рассчитать уравнение множественной регрессии для набора данных. Есть несколько файлов одинаковой структуры с данными для разных районов. Предполагаю организовать цикл, в котором последовательно подгружать файлы, выполнять расчеты и записывать результат в файл с таким же именем, но другим расширением. Составил следующую программку:

Код: Выделить всё


f <- c("_M_.dat","basin.DAT","Land.DAT","Mount.DAT","trogs.DAT","flat.DAT", "island.DAT","sedim.DAT","RF2.DAT") # составляем вектор из имен файлов
for(m in 1:9) {
z <- paste(f[m],"txt",sep=".") # задаем имя выходного файла
x <- read.table(f[m], sep=' ', header= TRUE) # читаем исходные данные
y <- data.frame(M=abs(x$M), X1=x$X1, X2= x$X2, X3=x$X3/-1000, X4=(x$X3/-1000)* x$X1* x$X2) # подготавливаем таблицу для рассчета
row.names(y) <- x$N
fit <- lm(M ~ X1+ X2 + X3 + X4, data=y) # рассчитываем уравнение регрессии
write(f[m],file=z, append = FALSE)
summary(fit)
write(fit$coefficients,file=z,append = TRUE) # пытаемся записать информацию в файл
}

Первая проблема: цикл не работает, хотя если в тело цикла добавить строки: результат выводится.
Вторая проблема: не получается записать fit в файл. Посмотрел структуру fit, там полно переменных. Переменную $coefficients удалось записать, при подстановке остальных выдает либо ошибку, либо печатает два пустых байта в файл.
А я планировал получит запись, аналогичную выводу функции summary(fit).
Подскажите, пожалуйста, что я делаю неправильно.

Иван Стрельников
Интересующийся
Сообщения: 40
Зарегистрирован: 11 авг 2011, 13:23
Репутация: 15

Re: Рассчет множественной регрессии для серии файлов с данны

Сообщение Иван Стрельников » 12 фев 2016, 19:31

Если поможет, текстовый вывод функции можно перенаправить в файл.

Пример:

Код: Выделить всё

sink(file=z, append = TRUE) # Устанавливаем соединение с файлом
summary(fit) # Перенаправляем вывод в файл
sink(NULL)

gamm
Гуру
Сообщения: 2659
Зарегистрирован: 15 окт 2010, 08:33
Репутация: 573
Ваше звание: программист
Откуда: Казань

Re: Рассчет множественной регрессии для серии файлов с данны

Сообщение gamm » 12 фев 2016, 22:06

Нужно выводить коэфициенты не из fit, а из summary, записав результат summary в переменную. А потом посмотреть, где они лежат командой str и выводить коэфициенты командой write.table(var.with.summary$var.with.coef, file=z,append=TRUE)

[ Сообщение с мобильного устройства ]

Иван Стрельников
Интересующийся
Сообщения: 40
Зарегистрирован: 11 авг 2011, 13:23
Репутация: 15

Re: Рассчет множественной регрессии для серии файлов с данны

Сообщение Иван Стрельников » 12 фев 2016, 23:47

И, если позволите вставить пять копеек, дописывать что-то поверх таблиц в .txt выгляди как-то не кошерно. Я бы сделал отдну табличку для всех датафайлов с примерно следующими колонками:
  • Имя файла с данными;
    Вызов функции (если он разный);
    Пересечение;
    Коэффициент для переменной 1;
    Коэффициент для переменной 2;
    и т.д. для всех переменных …
    Значимость коэффициентов (то же для всех переменных);
    Скорректированный R2;
    + все, что пожелаете, например RMSE;
А потом сохранить это в отдельный табличный файл. Так, у Вас все показатели будут в одном месте и с общей структурой. Но это на Ваше усмотрение.

gamm
Гуру
Сообщения: 2659
Зарегистрирован: 15 окт 2010, 08:33
Репутация: 573
Ваше звание: программист
Откуда: Казань

Re: Рассчет множественной регрессии для серии файлов с данны

Сообщение gamm » 13 фев 2016, 05:09

Иван Стрельников писал(а):дописывать что-то поверх таблиц в .txt выгляди как-то не кошерно
почему поверх? дописывание идет последовательно, получается вполне читаемый текстовый документ
Коэффициент для переменной 1;
Коэффициент для переменной 2;
и т.д. для всех переменных …
Значимость коэффициентов (то же для всех переменных);
это все есть в табличке, которую строит summary()

womak
Участник
Сообщения: 73
Зарегистрирован: 13 окт 2006, 06:19
Репутация: 0

Re: Рассчет множественной регрессии для серии файлов с данны

Сообщение womak » 13 фев 2016, 18:55

Коллеги, спасибо за ответы!
Иван Стрельников, как я понял, функция sink перенаправляет вывод с консоли в файл, а sink(NULL) отключает это действие. Только не пойму почему приведенный в Вашем ответе код не работает в теле цикла? Если запускать только его - работает, а в цикле - нет.
Идея про сводную табличку понравилась, но мне пока не хватает знаний ее осуществить. На данном этапе мне проще создавать для каждого файла "*.dat", файл "*.txt" с таким же именем и в нем прописывать все параметры уравнения регрессии.
gamm , изучил структуру переменной tmp <- summary(fit):
str(tmp)Show

Код: Выделить всё

List of 11
 $ call         : language lm(formula = M ~ X1 + X2 + X3 + X4, data = y)
 $ terms        :Classes 'terms', 'formula' length 3 M ~ X1 + X2 + X3 + X4
  .. ..- attr(*, "variables")= language list(M, X1, X2, X3, X4)
  .. ..- attr(*, "factors")= int [1:5, 1:4] 0 1 0 0 0 0 0 1 0 0 ...
  .. .. ..- attr(*, "dimnames")=List of 2
  .. .. .. ..$ : chr [1:5] "M" "X1" "X2" "X3" ...
  .. .. .. ..$ : chr [1:4] "X1" "X2" "X3" "X4"
  .. ..- attr(*, "term.labels")= chr [1:4] "X1" "X2" "X3" "X4"
  .. ..- attr(*, "order")= int [1:4] 1 1 1 1
  .. ..- attr(*, "intercept")= int 1
  .. ..- attr(*, "response")= int 1
  .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv> 
  .. ..- attr(*, "predvars")= language list(M, X1, X2, X3, X4)
  .. ..- attr(*, "dataClasses")= Named chr [1:5] "numeric" "numeric" "numeric" "numeric" ...
  .. .. ..- attr(*, "names")= chr [1:5] "M" "X1" "X2" "X3" ...
 $ residuals    : Named num [1:110] 0.386 -8.718 0.512 2.419 -1.084 ...
  ..- attr(*, "names")= chr [1:110] "4" "106" "209" "296" ...
 $ coefficients : num [1:5, 1:4] 29.57842 -0.00797 -0.98101 -3.09812 0.00418 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : chr [1:5] "(Intercept)" "X1" "X2" "X3" ...
  .. ..$ : chr [1:4] "Estimate" "Std. Error" "t value" "Pr(>|t|)"
 $ aliased      : Named logi [1:5] FALSE FALSE FALSE FALSE FALSE
  ..- attr(*, "names")= chr [1:5] "(Intercept)" "X1" "X2" "X3" ...
 $ sigma        : num 3.41
 $ df           : int [1:3] 5 105 5
 $ r.squared    : num 0.187
 $ adj.r.squared: num 0.156
 $ fstatistic   : Named num [1:3] 6.05 4 105
  ..- attr(*, "names")= chr [1:3] "value" "numdf" "dendf"
 $ cov.unscaled : num [1:5, 1:5] 0.158236 -0.000747 -0.022853 -0.087799 0.000326 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : chr [1:5] "(Intercept)" "X1" "X2" "X3" ...
  .. ..$ : chr [1:5] "(Intercept)" "X1" "X2" "X3" ...
 - attr(*, "class")= chr "summary.lm"
Записал коэффициенты командой:
write.table(tmp$ coefficients,Show

Код: Выделить всё

"Estimate" "Std. Error" "t value" "Pr(>|t|)"
"(Intercept)" 38.1392347130519 3.8014273827735 10.0328720958562 3.00079223536808e-10
"X1" 0.150173362242343 0.111867927570546 1.34241659342121 0.191524889798156
"X2" -1.06198663649355 0.334867081685138 -3.17136766967167 0.00398571218033112
"X3" -6.23140100509138 1.17766745789259 -5.291307799438 1.7571366734917e-05
"X4" -0.00946316959389718 0.00975082867033636 -0.970499012323508 0.341094404485524
Но не могу найти где лежит информация об максимальном и минимальном расхождениях между теоретической моделью и данными, которые выводятся функцией
summary(fit)Show
Call:
lm(formula = M ~ X1 + X2 + X3 + X4, data = y)

Residuals:
Min 1Q Median 3Q Max
-3.8955 -1.5621 -0.2931 1.2948 8.6130


Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 38.139235 3.801427 10.033 3.00e-10 ***
X1 0.150173 0.111868 1.342 0.19152
X2 -1.061987 0.334867 -3.171 0.00399 **
X3 -6.231401 1.177667 -5.291 1.76e-05 ***
X4 -0.009463 0.009751 -0.970 0.34109
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.597 on 25 degrees of freedom
Multiple R-squared: 0.6133, Adjusted R-squared: 0.5514
F-statistic: 9.911 on 4 and 25 DF, p-value: 6.031e-05
Возможно она рассчитывает этот параметр налету, используя данные $ residuals
Чисто теоретически хотел узнать, если потребуется записать в файл переменные $ call или $ terms,
то какой командой это можно сделать?Show

Код: Выделить всё

write.table(tmp$ call, file=z,append=TRUE)
Error in as.data.frame.default(x[[i]], optional = TRUE) : 
  cannot coerce class ""call"" to a data.frame

gamm
Гуру
Сообщения: 2659
Зарегистрирован: 15 окт 2010, 08:33
Репутация: 573
Ваше звание: программист
Откуда: Казань

Re: Рассчет множественной регрессии для серии файлов с данны

Сообщение gamm » 13 фев 2016, 20:51

womak писал(а):Только не пойму почему код не работает в теле цикла?
патамучта :mrgreen: нужно делать

Код: Выделить всё

print(summary(...))
. Я больше скажу, и для графиков в цикле нужно делать print(), а для графиков из пакета lattice и прочих подобных еще больше иногда изгаляться. И не забывайте писать

Код: Выделить всё

sink(z,append=TRUE))
Возможно она рассчитывает этот параметр налету, используя данные $ residuals

Код: Выделить всё

sink(...,append=TRUE)
print(range(my.model$residuals))
sink()
Чисто теоретически хотел узнать, если потребуется записать в файл переменные $ call или $ terms,
terms вам не нужен от слова совсем - если уж очень хочется, выводите начало модельной матрицы

Код: Выделить всё

sink(...,append=TRUE)
head(model.matrix(my.model))
sink()
Для вызова $call просто печать

Код: Выделить всё

sink(...,append=TRUE)
print(my.summary$call)
sink()

womak
Участник
Сообщения: 73
Зарегистрирован: 13 окт 2006, 06:19
Репутация: 0

Re: Рассчет множественной регрессии для серии файлов с данны

Сообщение womak » 15 фев 2016, 08:01

Спасибо, gamm !
gamm писал(а):

Код: Выделить всё

sink(...,append=TRUE)
print(range(my.model$residuals))
sink()
Вроде все просто. Получить запись аналогичную:

Код: Выделить всё

Residuals:
    Min      1Q  Median      3Q     Max 
-7.7468 -1.8409 -0.4698  1.4052  6.9903 
можно:

Код: Выделить всё

z1<-range(fit$residuals)
z2<-quantile(fit$residuals, c(.25, .75))
z3<-median(fit$residuals)
tmp<-c(Res.Min=z1[1],z2[1],Median=z3,z2[2],Res.Max=z1[2])
tmp
  Res.Min        25%     Median        75%    Res.Max 
-7.7468376 -1.8409315 -0.4698456  1.4051511  6.9903028 
Так и табличку, включающие различные формулы регрессии, недолго насобирать....

womak
Участник
Сообщения: 73
Зарегистрирован: 13 окт 2006, 06:19
Репутация: 0

Re: Рассчет множественной регрессии для серии файлов с данны

Сообщение womak » 18 фев 2016, 03:31

Коллеги, хотелось бы узнать, существует ли возможность проверить, как ведет себя подбираемая по формуле множественной регрессии зависимость за пределами диапазона исходных данных, участвующих в определении этой зависимости? Простите за мой математический...
Поясню.
Зависимость М определена по трем переменным Х1[-250;250] , Х2[0;6], Х3[0;7], а исходные данные, которые будут участвовать в расчетах лежат в диапазонах Х1[-350;350] , Х2[0;8], Х3[-1;9]. Ну и, соответственно, при реальном счете комбинации переменных могут быть самые разные, в результате чего искомый параметр М выходит за допустимый диапазон.
Я понимаю, что это результат того, что я увлекся полиномами и преобразованиями зависимой переменной.

gamm
Гуру
Сообщения: 2659
Зарегистрирован: 15 окт 2010, 08:33
Репутация: 573
Ваше звание: программист
Откуда: Казань

Re: Рассчет множественной регрессии для серии файлов с данны

Сообщение gamm » 18 фев 2016, 07:31

womak писал(а):как ведет себя подбираемая по формуле ... зависимость ... увлекся полиномами и преобразованиями
это зависит от того, путаете вы fit и modeling, или нет. Судя по "подбирает" и "полиномами" - первое. Больше ничего, не зная задачи и не видя ни данных, ни того, что вы делаете, сказать нельзя. Из общих советов - использовать gam()/gamm() из пакета mgcv вместо полиномов, перекрестную проверку, эзотерику (типа Random Forest, если данных много - он по крайней мере за диапазон не выйдет) и пытаться таки делать modeling исходя из природы явления, а не из данных, поскольку плохая модель все лучше хорошей регрессии ... и да, волшебную кнопку еще не изобрели, нужно будет повозиться, особенно с моделью :mrgreen:

Иван Стрельников
Интересующийся
Сообщения: 40
Зарегистрирован: 11 авг 2011, 13:23
Репутация: 15

Re: Рассчет множественной регрессии для серии файлов с данны

Сообщение Иван Стрельников » 18 фев 2016, 15:53

Основная проблема при экстраполяции -- возрастающая "неуверенность" модели. Чем дальше Вы уходите от центра интервала, на котором обучали модель, тем выше ожидаемая ошибка. Причем растет она не линейно, а по степенному закону. Попробуйте сделать predict вашей модели с параметром interval="confidence". Наложите полученные интервалы (+-) на график Ваших данных и получите приблизительное представление о том, как у Вас растет ожидаемая ошибка. Можете прикинуть, на глазок, как она будет расти при экстраполяции за границы.

Если говорить строго, то экстраполировать, вообще не следует. Если на практике, то нужно смотреть на конкретный вариант. Например, если Ваша модель описывает связь между объемом металлического шарика и его массой, то, логично, что экстраполировать можно как угодно (выше 0, конечно). Если это, например, связь роста и веса человека, то экстраполяция может выдать весьма странные результаты. Я как-то слышал доклад на студенческой конференции, в котором девочка уверяла, что через пару десятков лет автомобильный транспорт вообще перестанет производить выбросы. Она посмотрела, что с вводом нормативов от Евро-0 до Евро-5, выбросы уменьшались почти линейно. И естественно, в будущем эта линия пересечет 0. Вот она ничтоже сумняшеся и сделала такой вывод. Это, конечно, клиника, но суть та же.

gamm
Гуру
Сообщения: 2659
Зарегистрирован: 15 окт 2010, 08:33
Репутация: 573
Ваше звание: программист
Откуда: Казань

Re: Рассчет множественной регрессии для серии файлов с данны

Сообщение gamm » 18 фев 2016, 18:24

Иван Стрельников писал(а):Основная проблема при экстраполяции -- возрастающая "неуверенность" модели. Чем дальше Вы уходите от центра интервала, на котором обучали модель, тем выше ожидаемая ошибка. Причем растет она не линейно, а по степенному закону.
это не так, если модель соответствует физической природе явления. Тогда ошибка вообще не растет ... а вот если вместо адекватной модели (или на худой конец "правильных" сплайнов) использовать полиномы, то да - получим "свечку". Поэтому при невозможности построить адекватную модель используют что-то консервативное, типа shrinking splines, функций с конечным носителем (типа MBA), и т.д.

поэтому главная проблема - сочинить адекватную модель

Ответить

Вернуться в «R»

Кто сейчас на конференции

Сейчас этот форум просматривают: нет зарегистрированных пользователей и 1 гость