Страница 1 из 1
Рассчет множественной регрессии для серии файлов с данными
Добавлено: 12 фев 2016, 18:09
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) # пытаемся записать информацию в файл
}
Первая проблема: цикл не работает, хотя если в тело цикла добавить строки:
результат выводится.
Вторая проблема: не получается записать
fit в файл. Посмотрел структуру
fit, там полно переменных. Переменную
$coefficients удалось записать, при подстановке остальных выдает либо ошибку, либо печатает два пустых байта в файл.
А я планировал получит запись, аналогичную выводу функции
summary(fit).
Подскажите, пожалуйста, что я делаю неправильно.
Re: Рассчет множественной регрессии для серии файлов с данны
Добавлено: 12 фев 2016, 19:31
Иван Стрельников
Если поможет, текстовый вывод функции можно перенаправить в файл.
Пример:
Код: Выделить всё
sink(file=z, append = TRUE) # Устанавливаем соединение с файлом
summary(fit) # Перенаправляем вывод в файл
sink(NULL)
Re: Рассчет множественной регрессии для серии файлов с данны
Добавлено: 12 фев 2016, 22:06
gamm
Нужно выводить коэфициенты не из fit, а из summary, записав результат summary в переменную. А потом посмотреть, где они лежат командой str и выводить коэфициенты командой write.table(var.with.summary$var.with.coef, file=z,append=TRUE)
[ Сообщение с мобильного устройства ]
Re: Рассчет множественной регрессии для серии файлов с данны
Добавлено: 12 фев 2016, 23:47
Иван Стрельников
И, если позволите вставить пять копеек, дописывать что-то поверх таблиц в .txt выгляди как-то не кошерно. Я бы сделал отдну табличку для всех датафайлов с примерно следующими колонками:
- Имя файла с данными;
Вызов функции (если он разный);
Пересечение;
Коэффициент для переменной 1;
Коэффициент для переменной 2;
и т.д. для всех переменных …
Значимость коэффициентов (то же для всех переменных);
Скорректированный R2;
+ все, что пожелаете, например RMSE;
А потом сохранить это в отдельный табличный файл. Так, у Вас все показатели будут в одном месте и с общей структурой. Но это на Ваше усмотрение.
Re: Рассчет множественной регрессии для серии файлов с данны
Добавлено: 13 фев 2016, 05:09
gamm
Иван Стрельников писал(а):дописывать что-то поверх таблиц в .txt выгляди как-то не кошерно
почему поверх? дописывание идет последовательно, получается вполне читаемый текстовый документ
Коэффициент для переменной 1;
Коэффициент для переменной 2;
и т.д. для всех переменных …
Значимость коэффициентов (то же для всех переменных);
это все есть в табличке, которую строит summary()
Re: Рассчет множественной регрессии для серии файлов с данны
Добавлено: 13 фев 2016, 18:55
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
Re: Рассчет множественной регрессии для серии файлов с данны
Добавлено: 13 фев 2016, 20:51
gamm
womak писал(а):Только не пойму почему код не работает в теле цикла?
патамучта
нужно делать
. Я больше скажу, и для графиков в цикле нужно делать print(), а для графиков из пакета lattice и прочих подобных еще больше иногда изгаляться. И не забывайте писать
Возможно она рассчитывает этот параметр налету, используя данные $ 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()
Re: Рассчет множественной регрессии для серии файлов с данны
Добавлено: 15 фев 2016, 08:01
womak
Спасибо,
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
Так и табличку, включающие различные формулы регрессии, недолго насобирать....
Re: Рассчет множественной регрессии для серии файлов с данны
Добавлено: 18 фев 2016, 03:31
womak
Коллеги, хотелось бы узнать, существует ли возможность проверить, как ведет себя подбираемая по формуле множественной регрессии зависимость за пределами диапазона исходных данных, участвующих в определении этой зависимости? Простите за мой математический...
Поясню.
Зависимость М определена по трем переменным Х1[-250;250] , Х2[0;6], Х3[0;7], а исходные данные, которые будут участвовать в расчетах лежат в диапазонах Х1[-350;350] , Х2[0;8], Х3[-1;9]. Ну и, соответственно, при реальном счете комбинации переменных могут быть самые разные, в результате чего искомый параметр М выходит за допустимый диапазон.
Я понимаю, что это результат того, что я увлекся полиномами и преобразованиями зависимой переменной.
Re: Рассчет множественной регрессии для серии файлов с данны
Добавлено: 18 фев 2016, 07:31
gamm
womak писал(а):как ведет себя подбираемая по формуле ... зависимость ... увлекся полиномами и преобразованиями
это зависит от того, путаете вы fit и modeling, или нет. Судя по "подбирает" и "полиномами" - первое. Больше ничего, не зная задачи и не видя ни данных, ни того, что вы делаете, сказать нельзя. Из общих советов - использовать gam()/gamm() из пакета mgcv вместо полиномов, перекрестную проверку, эзотерику (типа Random Forest, если данных много - он по крайней мере за диапазон не выйдет) и пытаться таки делать modeling исходя из природы явления, а не из данных, поскольку плохая модель все лучше хорошей регрессии ... и да, волшебную кнопку еще не изобрели, нужно будет повозиться, особенно с моделью
Re: Рассчет множественной регрессии для серии файлов с данны
Добавлено: 18 фев 2016, 15:53
Иван Стрельников
Основная проблема при экстраполяции -- возрастающая "неуверенность" модели. Чем дальше Вы уходите от центра интервала, на котором обучали модель, тем выше ожидаемая ошибка. Причем растет она не линейно, а по степенному закону. Попробуйте сделать predict вашей модели с параметром interval="confidence". Наложите полученные интервалы (+-) на график Ваших данных и получите приблизительное представление о том, как у Вас растет ожидаемая ошибка. Можете прикинуть, на глазок, как она будет расти при экстраполяции за границы.
Если говорить строго, то экстраполировать, вообще не следует. Если на практике, то нужно смотреть на конкретный вариант. Например, если Ваша модель описывает связь между объемом металлического шарика и его массой, то, логично, что экстраполировать можно как угодно (выше 0, конечно). Если это, например, связь роста и веса человека, то экстраполяция может выдать весьма странные результаты. Я как-то слышал доклад на студенческой конференции, в котором девочка уверяла, что через пару десятков лет автомобильный транспорт вообще перестанет производить выбросы. Она посмотрела, что с вводом нормативов от Евро-0 до Евро-5, выбросы уменьшались почти линейно. И естественно, в будущем эта линия пересечет 0. Вот она ничтоже сумняшеся и сделала такой вывод. Это, конечно, клиника, но суть та же.
Re: Рассчет множественной регрессии для серии файлов с данны
Добавлено: 18 фев 2016, 18:24
gamm
Иван Стрельников писал(а):Основная проблема при экстраполяции -- возрастающая "неуверенность" модели. Чем дальше Вы уходите от центра интервала, на котором обучали модель, тем выше ожидаемая ошибка. Причем растет она не линейно, а по степенному закону.
это не так, если модель соответствует физической природе явления. Тогда ошибка вообще не растет ... а вот если вместо адекватной модели (или на худой конец "правильных" сплайнов) использовать полиномы, то да - получим "свечку". Поэтому при невозможности построить адекватную модель используют что-то консервативное, типа shrinking splines, функций с конечным носителем (типа MBA), и т.д.
поэтому главная проблема - сочинить адекватную модель