Я пытаюсь смоделировать набор данных, который соответствует имеющимся у меня эмпирическим данным, но я не уверен, как оценить ошибки в исходных данных. Эмпирические данные включают гетероскедастичность, но я не заинтересован в ее преобразовании, а скорее использую линейную модель с ошибочным термином для воспроизведения моделирования эмпирических данных.
Например, допустим, у меня есть эмпирический набор данных и модель:
n=rep(1:100,2)
a=0
b = 1
sigma2 = n^1.3
eps = rnorm(n,mean=0,sd=sqrt(sigma2))
y=a+b*n + eps
mod <- lm(y ~ n)
используя plot(n,y)
мы получаем следующее.
Однако, если я попытаюсь смоделировать данные, simulate(mod)
гетероскедастичность будет удалена и не будет учтена моделью.
Я могу использовать обобщенную модель наименьших квадратов
VMat <- varFixed(~n)
mod2 = gls(y ~ n, weights = VMat)
это обеспечивает лучшее соответствие модели на основе AIC, но я не знаю, как имитировать данные, используя выходные данные.
Мой вопрос заключается в том, как мне создать модель, которая позволит мне моделировать данные в соответствии с исходными эмпирическими данными (n и y выше). В частности, мне нужен способ оценить sigma2, ошибку, используя любой из них с использованием модели?
источник
Ответы:
Чтобы смоделировать данные с изменяющейся дисперсией ошибки, необходимо указать процесс создания данных для дисперсии ошибки. Как было отмечено в комментариях, вы сделали это, когда генерировали исходные данные. Если у вас есть реальные данные и вы хотите попробовать это, вам просто нужно определить функцию, которая определяет, как остаточная дисперсия зависит от ваших ковариат. Стандартный способ сделать это состоит в том, чтобы соответствовать вашей модели, проверить ее обоснованность (кроме гетероскедастичности) и сохранить остатки. Эти остатки становятся переменной Y новой модели. Ниже я сделал это для вашего процесса генерации данных. (Я не вижу, где вы устанавливаете случайное семя, так что это не будут буквально те же данные, но должны быть похожими, и вы можете точно воспроизвести мои, используя мое семя.)
Обратите внимание , что
R
«s ? Plot.lm даст вам участок (ср, здесь ) квадратный корень из абсолютных значений разностей, услужливо наложенный с lowess приступом, который является именно то , что вам нужно. (Если у вас есть несколько ковариат, возможно, вы захотите оценить их по каждому ковариате отдельно.) Существует малейший намек на кривую, но, похоже, что прямая линия хорошо подходит для подгонки данных. Итак, давайте точно подгоним эту модель:Нам не нужно беспокоиться о том, что остаточная дисперсия, похоже, увеличивается на графике расположения шкалы и для этой модели - что, по сути, должно произойти. Снова есть малейший намек на кривую, поэтому мы можем попытаться подогнать квадратное слагаемое и посмотреть, помогает ли это (но не помогает):
Если мы удовлетворены этим, теперь мы можем использовать этот процесс в качестве дополнения для имитации данных.
Обратите внимание, что этот процесс не гарантирует более точного поиска истинного процесса генерирования данных, чем любой другой статистический метод. Вы использовали нелинейную функцию для генерации SD ошибок, и мы аппроксимировали ее линейной функцией. Если вы на самом деле знаете истинный процесс генерации данных априори (как в этом случае, потому что вы имитировали исходные данные), вы можете также использовать его. Вы можете решить, достаточно ли приближенное приближение для ваших целей. Однако мы, как правило, не знаем истинного процесса генерации данных и, основываясь на бритве Оккама, используем простейшую функцию, которая адекватно соответствует данным, которые мы предоставили, с учетом объема доступной информации. Вы также можете попробовать сплайны или более интересные подходы, если хотите. Двусторонние распределения выглядят достаточно похожими на меня,
источник
Вам нужно смоделировать гетероскедастичность. Один из подходов - использование R-пакета (CRAN)
dglm
, обобщенной дисперсионной линейной модели. Это расширение glm, которое, в дополнение к обычномуglm
, подходит для второго glm для рассеивания остатков от первого glm. У меня нет опыта работы с такими моделями, но они кажутся многообещающими ... Вот код:Смоделированный участок показан ниже:
Сюжет выглядит так, как будто в симуляции использовалась оценочная дисперсия, но я не уверен, поскольку у функции simulate () нет методов для dglm's ...
(Еще одна возможность изучить использование
R
пакетаgamlss
, который использует другой подход к моделированию дисперсии как функции ковариабельных переменных.)источник