Когнитивная психология и эмоции
Субъективные заметки постдока-психолога
Записи в рубрике «регрессия»

Описание результатов glm из R в MS Word

20 июня 2011, заметки Метки: , , , ,

Суть проблемы: glm при использовании не-гауссовских распределений (пуассоновское, биномиальное) не выдает информации о коэффициенте псевдо-R2. У этого есть определенные основания (которые я пока до конца не понял), но, тем не менее, обычно этот коэффициент требуется в публикациях. Для его вычисления можно использовать функцию pR2. Ниже мини-функция, которая вставляет в Word описательную таблицу регрессии и статистики самой модели.

library(R2wd) #пакет для работы с Word
library(pscl) #пакет для pR2
 
wdGet() #подключение к Word
 
# сама функция
describe.glm<-function (fit){
  cb<-cbind(LogLik=logLik(fit),deviance=fit$null.deviance-fit$deviance, df=fit$df.null-fit$df.residual, p=1-pchisq(-fit$deviance+fit$null.deviance, -fit$df.residual+fit$df.null ),"Pseudo R-squared" = pR2(fit)[4])
  summ<-cbind(coef(summary(fit)),NA)
  summ<-round(summ,digits=3)
  summ[,5]<-symnum(summ[,4], cutpoints=c(0,0.001,0.01,0.05,0.1,1),symbols=c('***', '**', '*', '.','' ), legend=F)
  summ[,4]<-ifelse(summ[,4]==0,"< 0.001",summ[,4])
  wdTable(summ,align=c("l",rep("r",ncol(summ)-1),"l"))
  wdNormal(paste("LL = ",round(cb[1],digits=2),", /chi2(",cb[3],") = ",round(cb[2],digits=2),", p = ",round(cb[4],digits=3),", /r2 = ",round(cb[5],digits=2),"",sep=""))
}
 
#пример с логистической регрессией из Venables and Ripley (2002, pp. 190-2.)
ldose <- rep(0:5, 2)
numdead <- c(1, 4, 9, 13, 18, 20, 0, 2, 6, 10, 12, 16)
sex <- factor(rep(c("M", "F"), c(6, 6)))
SF <- cbind(numdead, numalive=20-numdead)
budworm.lg <- glm(SF ~ sex*ldose, family=binomial)
summary(budworm.lg)
 
#Call:
#glm(formula = SF ~ sex * ldose, family = binomial)
#
#Deviance Residuals: 
#     Min        1Q    Median        3Q       Max  
#-1.39849  -0.32094  -0.07592   0.38220   1.10375  
#
#Coefficients:
#            Estimate Std. Error z value Pr(>|z|)    
#(Intercept)  -2.9935     0.5527  -5.416 6.09e-08 ***
#sexM          0.1750     0.7783   0.225    0.822    
#ldose         0.9060     0.1671   5.422 5.89e-08 ***
#sexM:ldose    0.3529     0.2700   1.307    0.191    
#---
#Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
#
#(Dispersion parameter for binomial family taken to be 1)
#
#    Null deviance: 124.8756  on 11  degrees of freedom
#Residual deviance:   4.9937  on  8  degrees of freedom
#AIC: 43.104
#
#Number of Fisher Scoring iterations: 4
#
 
describe.glm(budworm.lg)

Результат в Word (на самом деле, он там даже симпатичнее):

  Estimate Std. Error z value Pr(>|z|)  
(Intercept) -2.994 0.553 -5.416 < 0.001 ***
sexM 0.175 0.778 0.225 0.822  
ldose 0.906 0.167 5.422 < 0.001 ***
sexM:ldose 0.353 0.27 1.307 0.191  

LL = -17.55, /chi2(3) = 119.88, p = < 0.001, /r2 = 0.77 -------------- /chi2 потом меняется на нормальную запись χ2 с верхней двойкой, /r2 — на R2.
Еще б понять, почему после таблицы пустая строка вылезает.

Логистическая регрессия в R

23 октября 2010, наука Метки: , , ,

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

продолжение »