corr.bino <- function(modelo=fit.model) { # # Descrição e detalhes: # Esta função calcula o coeficiente de correlação intra-classe (ou intra-grupo, rho) que usa os resíduos de Pearson # de um modelo com distribuição binomial em que as unidades experimentais são os g grupos. # # A função também calcula o valor da estatística de escore desenvolvida por Paula e Artes (2000) para testar # H0: rho=0 contra H1: rho>0 num modelo logístico para dados agrupados. Pode-se ver mais da teoria envolvida em Paula # (2003, págs.158-162) e em McCullagh e Nelder (1989, págs.124-128). # # Os dados devem estar disponíveis pelo comando attach( ). # # Argumentos obrigatórios: # modelo: deve-se informar o objeto onde está o ajuste do modelo com distribuição binomial. # # A função retorna os seguintes valores: a correlação intra-classe (rho), estatística escore (escore) e o # nível descritivo (pvalue). # # Autor: Frederico Zanqueta Poleto , arquivo disponível em http://www.poleto.com # # Referências: # MCCULLAGH, P. e NELDER, J. A. (1989). Generalized Linear Models. 2ª ed. Chapman and Hall, London. # PAULA, G. A. (2003). Modelos de Regressão com apoio computacional. IME-USP, São Paulo. [Não publicado, # disponível em http://www.ime.usp.br/~giapaula/Book.pdf] # PAULA, G. A. e ARTES, R. (2000). One-sided test to assess correlation in logistic linear models using estimating # equations. Biometrical Journal 42, 701-714. # # Exemplo: # corr.bino(ajuste) # if(class(modelo)[1] != "glm") { stop(paste("\nA classe do objeto deveria ser glm e nao ",class(modelo),"!!!\n")) } if(modelo$family[[1]] != "Binomial" & modelo$family[[1]] != "binomial") { stop(paste("\nA familia do objeto deveria ser Binomial e nao ",modelo$family[[1]],"!!!\n")) } y<-modelo$y n<-length(y) p <- length(modelo$coef) if(is.null(modelo$prior.weights)) { ntot<-rep(1,n) } else { ntot<-modelo$prior.weights } y<-round(y*ntot,0) fr<-fitted(modelo) nt1 <- 0.5*sum(ntot*(ntot-1)) sum1 <- (0.5*y*(y-1) - fr*(ntot-1)*y + 0.5*fr*fr*ntot*(ntot-1))/(fr*(1-fr)) sum2 <- sum(sum1*sum1) sum1 <- sum(sum1) rho <- sum1/(nt1 - p) escore <- sum1/sqrt(sum2) pvalue <- 1-pnorm(escore) list(rho=rho,escore=escore,pvalue=pvalue) }