применить и ошибку функции boxcox в R

У меня проблема с использованием функции apply при вызове функции boxcox.

Набор данных состоит из 35 образцов, 3 процедур в течение 2 дней (некоторые люди повторяются в течение нескольких дней, но всегда в группе лечения) с ~ 1000 измерений / зависимых переменных, я включил сюда только 5 (# 1,2,4,7,9 ).

NUM day trt 1   2   4   7   9  
8145    7   L   0.986   0.423   0.86    0.648   2.031  
8169    7   L   1.013   0.512   1.157   0.633   0.55  
8201    7   L   0.236   6.144   0.604   1.759   1.181  
8212    7   L   0.455   0.707   0.693   0.972   0.615  
8168    7   L   1.261   0.618   1.138   0.943   0.868  
8193    7   L   1.754   0.273   1.224   0.719   0.895  
8145    13  L   1.257   0.36    1.505   0.729   0.735  
8169    13  L   1.549   0.577   1.53    0.771   1.116  
8201    13  L   0.489   0.378   0.842   1.538   1.676  
8212    13  L   0.991   0.34    1.03    1.157   1.076  
8168    13  L   1.355   0.391   1.416   0.953   1.479  
8193    13  L   1.029   0.308   1.027   0.902   1.934  
8214    7   C   1.224   0.298   1.113   1.445   0.218  
8139    7   C   1.277   0.554   1.443   0.895   0.74  
8151    7   C   1.312   2.025   1.197   0.675   0.791  
8160    7   C   1.555   0.405   1.432   0.826   0.501  
8196    7   C   0.938   0.717   0.917   0.801   1.462  
8213    7   C   0.835   1.863   0.942   1.967   0.739  
8139    13  C   0.958   0.275   1.273   1.351   0.842  
8151    13  C   0.864   0.517   0.98    1.368   1.865  
8160    13  C   1.516   0.895   1.318   0.551   0.779  
8239    13  C   1.071   0.194   0.955   1.87    0.68  
8196    13  C   1.299   0.594   1.14    0.877   1.873  
8213    13  C   1.601   0.733   1.375   0.738   1.273  
8231    7   H   1.401   0.483   1.001   1.052   0.548  
8232    7   H   1.292   0.574   1.634   0.641   0.464  
8219    7   H   0.785   0.396   0.886   0.903   1.734  
8177    7   H   0.525   0.252   0.563   0.914   1.174  
8143    7   H   1.398   0.266   0.947   0.94    0.781  
8219    13  H   0.903   6.225   1.109   1.218   1.073  
8143    13  H   1.086   0.435   1.4 0.922   0.925  
8167    13  H   0.83    0.574   1.09    1.338   1.563  
8231    13  H   1.498   0.482   1.375   0.855   0.719  
8232    13  H   1.055   0.811   0.887   0.606   0.634  
8177    13  H   0.857   1.324   0.954   1.635   0.675  

Я пытаюсь запустить функцию boxcox в пакете MASS в R, чтобы оценить, нужно ли преобразовать данные, и если да, преобразовать в соответствии со значением lamda с максимальной вероятностью журнала (или его округлением).

linear.f=function(x){lm(x~day+trt+day*trt, data=data)}
linear.multiple=apply(data[,4:ncol(data)],2,linear.f)

работает нормально

и запуск функции boxcox по отдельности также работает нормально

boxcox(lm(x~day+trt+day*trt,data=data))$x[which.max(boxcox(lm(x~day+trt+day*trt, data=data))$y)]

ввод в функцию для применения:

lamda.f=function(x){boxcox(lm(x~day+trt+day*trt, data=data))$x[which.max(boxcox(lm(x~day+trt+day*trt, data=data))$y)]}

Однако проблема возникает при попытке запустить boxcox для всех столбцов / зависимых переменных с помощью apply:

lamda.multiple=apply(data[,4:ncol(data)],2,lamda.f)

ошибка:

 Error in model.frame.default(formula = x ~ day + trt + day * trt, data = data,  : 
  invalid type (list) for variable 'x'

Я делаю что-то явно неправильное? как boxcox будет работать, если он будет один за другим в каждом столбце данных ..

Я также пробовал это, используя цикл for вместо применения, но безуспешно

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


person Mike    schedule 10.08.2015    source источник


Ответы (1)


Проблема не связана с вызовом apply(). Это связано с тем, что вы используете x в своей формуле, но набор данных, который вы передаете lm(), не имеет столбца x.

Функция lm() немного «снисходительна», возможно, нежелательно. Если символ в данной формуле не может быть найден в данном наборе данных, то lm() позволяет ему связываться с любой переменной, которая может быть найдена в цепочке окружения замыкания формулы. Из https://stat.ethz.ch/R-manual/R-devel/library/stats/html/model.frame.html:

Все переменные в формуле, подмножестве и в ... сначала ищутся в данных, а затем в среде формулы (дополнительные сведения см. В справке по формуле ()) и собираются во фрейм данных.

...

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

Теперь, на самом деле, обе ваши linear.f() и lamda.f() функции имеют параметр функции x, и это позволяет выполнить вызов lm() в обеих функциях. IOW, когда lm() выполняется внутри этих функций, x не может привязаться к какому-либо столбцу в данном наборе данных (data), но затем привязывается к параметру функции в среде закрытия формулы (которая является средой оценки текущей функции оценка). Таким образом, вызов lm() завершается успешно.

Более подробно: когда функция оценивается, ее параметры всегда сохраняются в среде оценки, которая создается для этой конкретной оценки этой конкретной функции. Поскольку вы определили формулу буквально внутри функции, она замыкается вокруг текущей среды оценки, и поэтому при запуске lm() после неудачной попытки привязки к столбцу в наборе данных поиск цели для символа x попадает в эту среду. first и привязывается к записи параметра функции в среде.

Я заметил, что вы указали, что однократный вызов верхнего уровня _17 _ + _ 18_ был успешным. Это любопытно, поскольку нет параметра функции для привязки. Я могу только догадываться, что в то время, когда вы запускали эту строку кода, в глобальной среде у вас была переменная x, которая действительна для модели. В противном случае произойдет сбой, поскольку символ x в формуле не имеет ничего, с чем можно было бы связываться ни в наборе данных, ни в цепочке среды закрытия. Обратите внимание, что lm() должен завершиться ошибкой; Я даже не говорю о boxcox() здесь:

if (exists('x')) rm(x); ## remove global x
lm(x~day+trt+day*trt,data=data);
## Error in eval(expr, envir, enclos) : object 'x' not found

Теперь, когда вы выполнили свой второй вызов apply, мы, вероятно, можем догадаться, что вы перезаписали глобальную переменную x списком, что привело бы к точному сообщению об ошибке, которое вы получили:

x <- list();
lm(x~day+trt+day*trt,data=data);
## Error in model.frame.default(formula = x ~ day + trt + day * trt, data = data,  :
##   invalid type (list) for variable 'x'

Но здесь должно быть что-то еще. Вспомните, как я указывал, что и linear.f(), и lamda.f() имеют параметр функции x, поэтому lm() должен быть привязан к нему, и, поскольку это действительный вектор-предиктор, он должен быть успешным.

Мне кажется, что boxcox() вызывает lm() сам по себе и каким-то образом препятствует выполнению поиска цепочки замыканий в этом случае. Я считаю, что могу доказать эту гипотезу с помощью следующего кода, который использует IIFE (изначально JavaScript термин, но применимый к любому языку с первоклассными выражающими функциями, например R):

if (exists('x')) rm(x); ## remove global x
data2 <- data.frame(y=1:3); ## no x
(function(x) lm(y~x,data2))(1:3); ## x binds to function parameter
##
## Call:
## lm(formula = y ~ x, data = data2)
##
## Coefficients:
## (Intercept)            x
##           0            1
##
(function(x) boxcox(lm(y~x,data2)))(1:3);
## Error in eval(expr, envir, enclos) : object 'x' not found
traceback();
## 14: eval(expr, envir, enclos)
## 13: eval(predvars, data, env)
## 12: model.frame.default(formula = y ~ x, data = data2, drop.unused.levels = TRUE)
## 11: stats::model.frame(formula = y ~ x, data = data2, drop.unused.levels = TRUE)
## 10: eval(expr, envir, enclos)
## 9: eval(mf, parent.frame())
## 8: lm(formula = y ~ x, data = data2, y = TRUE, qr = TRUE)
## 7: eval(expr, envir, enclos)
## 6: eval(call, parent.frame())
## 5: update.default(object, y = TRUE, qr = TRUE, ...)
## 4: update(object, y = TRUE, qr = TRUE, ...)
## 3: boxcox.lm(lm(y ~ x, data2))
## 2: boxcox(lm(y ~ x, data2)) at #1
## 1: (function(x) boxcox(lm(y ~ x, data2)))(1:3)

Как видите, похоже, что boxcox() работает lm(), но, в отличие от прямого вызова lm(), не может выполнить привязку к параметру функции x, даже если это должно быть возможно, поскольку аргумент формулы, переданный в boxcox(), содержит указатель среды закрытия. Полагаю, мы могли бы списать это на слабость функции boxcox().

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

data2 <- data.frame(y=1:3,x=10:12);
boxcox(lm(y~x,data2)); ## succeeds

Однако в вашем случае это не так просто, потому что у вас есть динамическая переменная результата. Хорошим подходом было бы дополнение вашего data.frame переменной результата непосредственно перед передачей ее в вызов lm(), например с cbind(data,x).

К сожалению, и это удивительно, я получаю следующую ошибку:

lamda.f <- function(x) { data.x <- cbind(data,x); boxcox(lm(x~day+trt+day*trt,data=data.x))$x[which.max(boxcox(lm(x~day+trt+day*trt,data=data.x))$y)]; };
lamda.multiple <- apply(data[,4:ncol(data)],2,lamda.f);
## Error in is.data.frame(data) : object 'data.x' not found

Я даже попытался отделить вызов lm() от вызова boxcox(), на случай, если boxcox() выполнил какую-то сумасшедшую нестандартную оценку своих аргументов и, таким образом, попытался оценить data.x в контексте, который предотвратил бы привязку к среде оценки lamda.f():

lamda.f <- function(x) { data.x <- cbind(data,x); m <- lm(x~day+trt+day*trt,data=data.x); b <- boxcox(m); b$x[which.max(b$y)]; };
lamda.multiple <- apply(data[,4:ncol(data)],2,lamda.f);
## Error in is.data.frame(data) : object 'data.x' not found

Мне кажется, что boxcox() очень сильно зависит от глобальной среды. Запуск boxcox() изнутри области функции и в зависимости от каких-либо локальных переменных просто прерывает ее. Я подозреваю, что он проверяет вызов, хранящийся в объекте модели (например, m$call), и пытается напрямую прочитать символы. В любом случае, это действительно странно.

Я думаю, вы можете решить эту проблему, сохранив дополненный data.frame в глобальной среде и убедившись, что он существует там в момент запуска boxcox(). Для этого вы можете использовать оператор сверхназначения:

lamda.f <- function(x) { data.x <<- cbind(data,x); m <- lm(x~day+trt+day*trt,data=data.x); b <- boxcox(m); b$x[which.max(b$y)]; };
lamda.multiple <- apply(data[,4:ncol(data)],2,lamda.f);
lamda.multiple;
##         X1         X2         X4         X7         X9
##  1.2323232 -0.6666667  0.7474747 -0.6666667  0.2222222
person bgoldst    schedule 10.08.2015