Составьте профиль отклонения профиля для соответствия GLM в R - PullRequest
Новая
       15

Составьте профиль отклонения профиля для соответствия GLM в R

5 голосов
/ 20 марта 2012

Я хотел бы иметь возможность построить отклонение профиля для оценки параметра, установленной с использованием функции glm() в R. Отклонение профиля - это функция отклонения для различных значений рассматриваемой оценки параметра после оценки всех других параметров. , Мне нужно построить график отклонения для нескольких значений вокруг подогнанного параметра, чтобы проверить предположение о квадратичной функции отклонения.

Моя модель предсказывает повторное осуждение преступников. Формула имеет вид: reconv ~ [other variables] + sex, где reconv - бинарный фактор да / нет, а sex - бинарный мужской / женский фактор. Я хотел бы построить отклонение профиля параметра, оцененного для пола = женщины (пол = мужчина - это контрольный уровень).

Функция glm() оценивает параметр как -0,22 с ошибкой стандартного отклонения 0,12.

[Я задаю этот вопрос, потому что не было никакого ответа, который я мог найти, но я решил его и хотел опубликовать решение, которое будет полезно для других. Конечно, дополнительная помощь приветствуется. : -)]

Ответы [ 3 ]

6 голосов
/ 20 марта 2012

См. profileModel пакет от Иоанниса Космидиса.У него была статья в R Journal (появятся R News), иллюстрирующая пакет:

Иоаннис Космидис.Пакет profilemodel R: цели профилирования для моделей с линейными предикторами.R News, 8 (2): 12-18, октябрь 2008 г.

PDF-файл здесь (вся новостная рассылка).

5 голосов
/ 20 марта 2012

См. ?profile.glm (и example("profile.glm")) в пакете MASS - я думаю, что он будет делать все, что вы хотите (по умолчанию он не загружается, но упоминается в ?profile, который мог быть первое место, которое вы посмотрели ...) (Обратите внимание, что профили обычно отображаются в масштабе со знаком квадратного корня, так что действительно квадратичный профиль будет выглядеть как прямая линия.)

0 голосов
/ 23 марта 2012

Способ, который я нашел, включает использование функции offset () (как подробно описано в Pawitan, Y. (2001) «По всей вероятности», стр. 172). Ответы, данные @BenBolker и @GavinSimpson, лучше, чем эти, поскольку они ссылаются на пакеты, которые сделают все, что делают, и многое другое. Я пишу об этом, потому что это совсем по-другому, а также, что зарисовка «вручную» иногда полезна для обучения. Это многому меня научило. 

sexi <- as.numeric(data.frame$sex)-1      #recode a factor as 0/1 numeric

beta <- numeric(60)              #Set up vector to Store the betas
deviance <- numeric(60)          #Set up vector to Store the deviances

for (i in 1:60){

  beta[i] <- 0.5 - (0.01*i)  
  #A vector of values either side of the fitted MLE (in this case -0.22)

  mod <- update(model,
                   .~. - sex             #Get rid of the fitted variable
                   + offset(   I(sexi*beta[i])   )   #Replace with offset term.
                )
  deviance[i] <- mod$deviance                        #Store i'th deviance
}

best <- which.min(deviance)                   
#Find the index of best deviance. Should be the fitted value from the model

deviance0 <- deviance - deviance[best]         
#Scale deviance to zero by subtracting best deviance

betahat <- beta[best]    #Store best beta. Should be the fitted value.
stderror <- 0.12187      #Store the std error of sex, found in summary(model)

quadratic <- ((beta-betahat)^2)*(1/(stderror^2))  
#Quadratic reference function to check quadratic assumption against

x11()                                    
plot(beta,deviance0,type="l",xlab="Beta(sex)",ylim=c(0,4))    
lines(beta,quadratic,lty=2,col=3)           #Add quadratic reference line
abline(3.84,0,lty=3)                #Add line at Deviance = 3.84
...