не подходит dbinom для регрессии журнала

Я пытался подогнать регрессию журнала с четырьмя параметрами, используя dbinom. Логарифмическая регрессия с четырьмя параметрами выражается как: F(x) = d+(a/(1+exp((b-(x)/c))), где d=нижняя асимптота (ymin), a+d=верхняя асимптота , b = точка перегиба, c = наклон, Моя переменная ответа представляет собой отношение из набора данных обследования млекопитающих (Patch_richness / Richness_proportion), y принимает значения от 0 до 1 (0,8, 0,4, 0,25 ......) , моя цель - сравнить модель glm.null, glm и регрессию журнала с четырьмя параметрами с помощью AIC, чтобы найти, какой из трех подходит лучше всего. Не было проблем при запуске функции для y = count (Patch_Richness) с использованием dpois, а затем просто заменить значения coeff при построении кривой:

library(bbmle)
cerrado = read.csv("data_stack.csv")
attach(cerrado)
logip = function(p,lambda,x){
  a = p[1]
  b = p[2]
  c = p[3]
  d = p[4]
  Riq1 = d+(a/(1+exp((b-(FOREST500+km))/c)))
  -sum(dpois(x,lambda=Riq1, log=TRUE))
}
parnames(logip) = c("a","b","c","d")

modTR.log = mle2(minuslog = logip, start = c(a = 5,b = 72,c = 3,d = 0.1), data  = list(x = Patch_Richness))
summary(modTR.log)
plot(FOREST500,Patch_Richness, xlab = "Forest cover", ylab = "Patch Richness")#original data
curve (-0.29382+(4.95218/(1+exp((118.34117-x)/60.30478))), add=T)

Но у меня возникли проблемы при попытке подогнать эту функцию для y = proportion (Richness_prop)

logip = function(p, lambda, x){
  a = p[1]
  b = p[2]
  c = p[3]
  d = p[4]
  Riq1 = d+(a/(1+exp((b-(FOREST500 + km))/c)))
  -sum(dbinom(x,18,0.5,log = FALSE))
}
parnames(logip) = c("a","b","c","d")

modTR.log = mle2(minuslog = logip, start = c(a = 1,b = 72,c = 1,d = 0), data = list(x = Richness_prop))

summary(modTR.log)
AIC(modTR.log)

Модель запускается только при log = FALSE (с предупреждениями о нецелочисленных значениях), но итоговый вывод дает в качестве значений coeff точно такое же число, что и начальные начальные значения, независимо от числа в начальных значениях. Так что я думаю, что-то действительно плохо с этим. Правильно ли я устанавливаю параметры dbinom? Почему он работает только с log=FALSE?

данные

Очень ценю некоторую помощь Спасибо!


r logistic-regression poisson mle
person mmr09    schedule 25.03.2021    source источник

Ответы (1)


arrow_upward
0
arrow_downward

Я запустил вашу первую часть, но мне интересно, в порядке ли ваша модель? Я не вижу хорошего соответствия между данными и вашей моделью. Кроме того, у вашей функции logip был ненужный лямбда-аргумент, так как он вычисляется вами. Не могли бы вы использовать этот код и объяснить, почему вы думаете, что он работает? Потому что я вижу модель, у которой два параметра не удалось четко найти?

logip = function(p,x){
  a = p[1]
  b = p[2]
  c = p[3]
  d = p[4]
  Riq1 = d+(a/(1+exp((b-(FOREST500+km))/c)))
  -sum(dpois(x,lambda=Riq1, log=TRUE))

parnames(logip) = c("a","b","c","d")

modTR.log = mle2(minuslog = logip, start=list(a=5, b=72,c=3,d=0.1), data=list(x = 
                 Patch_Richness), vecpar=TRUE)

summary(modTR.log)
coef_fit = coef(modTR.log)
plot(FOREST500,Patch_Richness, xlab = "Forest cover", ylab = "Patch Richness")
curve (coef_fit["d"]+(coef_fit["a"]/(1+exp((coef_fit["b"]-x)/coef_fit["d"]))), add=T)

}

person flow_me_over    schedule 26.03.2021
comment
Спасибо за Ваш ответ. Запустил ваш код, и он работает точно так же, как мой. Тот же вывод. Я понимаю ненужную лямбду и знаю, что она может выглядеть плохо, это всего лишь одна из многих моделей, которые я запускал и сравнивал. Я не совсем понимаю ваш вопрос, но я хочу оценить эти коэффициенты для пропорции: y принимает значения от 0 до 1 (0,8, 0,4, 0,25 ......). Я предположил, что мне пришлось использовать dbinom вместо dpois, но у меня проблемы, и я не уверен, смогу ли я сделать что-то подобное для данных: - person mmr09; 31.03.2021
comment
Вы в основном делаете логистическую регрессию, используя MLE, верно? Я бы сначала использовал некоторые фиктивные данные, сгенерированные из точной модели в семействе моделей, с которым вы хотите соответствовать, чтобы убедиться, что подгонка работает так, как вы хотите. - person flow_me_over; 31.03.2021