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

[womak]

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

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

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)  # пытаемся записать информацию в файл

}

Первая проблема: цикл не работает, хотя если в тело цикла добавить строки:

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

m
f[m]

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

[Иван Стрельников]

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

Пример:

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

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

[gamm]

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

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

[Иван Стрельников]

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

• Имя файла с данными;
  Вызов функции (если он разный);
  Пересечение;
  Коэффициент для переменной 1;
  Коэффициент для переменной 2;
  и т.д. для всех переменных …
  Значимость коэффициентов (то же для всех переменных);
  Скорректированный R2;
  + все, что пожелаете, например RMSE;

А потом сохранить это в отдельный табличный файл. Так, у Вас все показатели будут в одном месте и с общей структурой. Но это на Ваше усмотрение.

[gamm]

дописывать что-то поверх таблиц в .txt выгляди как-то не кошерно

почему поверх? дописывание идет последовательно, получается вполне читаемый текстовый документ

Коэффициент для переменной 1;
Коэффициент для переменной 2;
и т.д. для всех переменных …
Значимость коэффициентов (то же для всех переменных);

это все есть в табличке, которую строит summary()

[womak]

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

Спойлер

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

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"

Записал коэффициенты командой:

Спойлер

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

"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

Но не могу найти где лежит информация об максимальном и минимальном расхождениях между теоретической моделью и данными, которые выводятся функцией

Спойлер

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,

Спойлер

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

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]

Только не пойму почему код не работает в теле цикла?

патамучта нужно делать

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

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]

Спасибо, 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]

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

[gamm]

как ведет себя подбираемая по формуле ... зависимость ... увлекся полиномами и преобразованиями

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

[Иван Стрельников]

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

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

[gamm]

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

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

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