2013-09-13 5 views
1

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

stackoverflow.com/questions/15275236/least-square-percentage-regression

Однако, это предполагает, что я знаю, насколько каждое наблюдение должно быть взвешенным априори. Я не. Я не знаю, если ошибка процент составляет 1%, 10%, 15% и т.д. То, что я хочу, это модель подходит в качестве

y= b1*x + e 

где термин ошибка моделируется как:

e= b2*x 

b2 будет процентной ошибкой, которая должна быть минимизирована в модели регрессии. Я не смог найти какой-либо пакет или какой-либо код, подходящий для модели такого типа для R. Любая оценка обратной связи о том, как это сделать, будет очень признательна.

+0

Это не эквивалентен использованию лога-преобразование у и продолжения обычная регрессия наименьших квадратов? Чтобы вернуться к нетрансформированному y, вы должны повысить уровень RHS, что приведет к мультипликативным терминам и ошибкам. – zkurtz

+0

@zkurtz log transform предполагает, что отношение принимает вид y = e^x, и поэтому размер эффекта от обратного преобразования будет уменьшаться, если соотношение действительно линейно. Преобразование журнала может быть подходящим для определения значимости, но не для размера эффекта. Я предпочел бы моделировать данные так, как они есть на самом деле, вместо того, чтобы преобразовывать их к не нормальным остаткам за счет оценки правильного размера эффекта. – colin

+0

Чтобы быть более явным, преобразование журнала обеспечило бы нормальное распределение ошибок и касалось гетероседастичности этого шаблона, чтобы его можно было обрабатывать в рамках OLS. Кроме того, в этой ситуации может оказаться оправданным взять размер эффекта из нетрансформированной модели. Однако то, что я предпочел бы сделать, это запустить одну модель, которая обрабатывает все это. Это было бы достигнуто путем выполнения процентной регрессии наименьших квадратов, а не обычной регрессии наименьших квадратов. – colin

ответ

4

Я предполагаю, что вы имеете в виду процент регрессии, как это определено Tofallis (2009).

Используя свой пример:

Sales <- c(6375,11626,14655,21869,26408,32406,35108,40295,70762,80553,95294,101314,116141,122316,141650,175026,230614,293543) 
Expenses <- c(62.5,92.9,178.3,258.4,494.7,1083,1620.6,421.7,509.2,6620.1,3918.6,1595.3,6107.5,4454.1,3163.8,.7,1703.8,9528.2) 

Если мы применим обычный метод наименьших квадратов с продаж в качестве зависимой переменной мы получаем модель продаж = 43942 + 15,00 R & D с р-значения 0,03 и 0,0015 для перехвата и наклона соответственно.

fit1 <- lm(Sales ~ Expenses) 
summary(fit1) 
#    Estimate Std. Error t value Pr(>|t|) 
# (Intercept) 43941.705 18493.079 2.376 0.03033 * 
# Expenses  14.994  3.915 3.830 0.00148 ** 

Если мы делаем это и выполнять обычные наименьших квадратов получаем модель: Ln (Sales) = 10,341 + 0.000198 R & D с р-значениями 0,002 для наклона и по существу ноль для перехвата.

fit2 <- lm(log(Sales) ~ Expenses) 
summary(fit2) 
#    Estimate Std. Error t value Pr(>|t|)  
# (Intercept) 1.034e+01 2.535e-01 40.793 < 2e-16 *** 
# Expenses 1.982e-04 5.366e-05 3.694 0.00197 ** 

Наконец, обратимся к подходу, представленному в данной статье, сводя к минимуму квадратов процентных остатков. Полученная модель оказывается равной после преобразования назад: Sales = 8817 + 17,88 R & D с p-значениями 0,002 и 5 × 10-5 для наклона и перехвата соответственно.

fit3 <- lm(Sales ~ Expenses, weights = 1/Sales^2) 
summary(fit3) 
#    Estimate Std. Error t value Pr(>|t|)  
# (Intercept) 8816.553 2421.644 3.641 0.0022 ** 
# Expenses  17.880  3.236 5.525 4.61e-05 *** 

Таким образом, в конце концов, это взвешенная регрессия.

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

resfun <- function(par) { 
    sum((Sales - par[[1]]*Expenses - par[[2]])^2/Sales^2) 
} 

optim(c(10,1000), resfun) 
# $par 
# [1] 17.87838 8816.44304 

optim(c(10,1000), resfun, method="BFGS") 
# $par 
# [1] 17.97975 8575.71156 

+0

Я имею в виду процент регрессии, как определено Tofallis 2009. Спасибо за добавление ссылки. Ваш ответ и оригинальная статья очень ясны. У меня есть одна проблема, хотя в документе Tofallis, которую вы связываете, предлагается взвешивать на (1/y), в то время как вы рекомендуете взвешивать (1/y^2). y - «Продажи» в вашем примере. Почему вы выбираете вес (1/y^2), а не (1/y), как рекомендовано в источнике? Для меня я, взвешивая обратный квадрат, а не только обратный, дал бы большие наблюдения еще меньшему весу, чего я не хочу делать! – colin

+0

@colin Взвешивание с 1/y^2 (с использованием функции 'lm') означает, что квадратные остатки умножаются на веса (т. Е. Остатки умножаются на 1/y). С быстрым взглядом это согласуется с ссылкой. Если вы не согласны, пожалуйста, укажите соответствующую часть статьи. Вы можете использовать любые веса, которые хотите, если это улучшит модель. – Roland

+0

спасибо. Я решил вес с 1/y, а не 1/y^2. В рукописи описание приведено в разделе «Вывод формул для коэффициентов», которое начинается в нижней части страницы 3, далее на стр. 4. Чтобы следить за проверкой нормальности, я считаю, что в таких ситуациях правильно посмотрите на график qqnorm остатков модели, деленный на квадратный корень из значений y, и то же самое должно быть сделано для остатков по сравнению с установленными значениями? – colin

0

Посмотрите на gls функции в nlme пакете, вместе с одним из varClasses таких как varIdent или varPower.

Возможно модель, как:

gls(y ~ x, data=mydata, weights=varPower(form= ~x)) 
Смежные вопросы