Ronaldo Reis Jr.
2002-Jul-02 16:31 UTC
[R] error in plot residuals in a glm with iterations.
Hi, I make this model:> moths.m6 <-glm(nind~metros+especie+habitat+metros:habitat+habitat:especie,family=poisson)> anova(moths.m6,test="F")Analysis of Deviance Table Model: poisson, link: log Response: nind Terms added sequentially (first to last) Df Deviance Resid. Df Resid. Dev F Pr(>F) NULL 81 488.32 metros 1 0.01 80 488.32 0.0083 0.9276077 especie 1 13.40 79 474.92 13.3967 0.0002521 *** habitat 7 185.93 72 288.99 26.5615 < 2.2e-16 *** metros:habitat 6 34.89 66 254.10 5.8143 4.535e-06 *** especie:habitat 7 99.24 59 154.86 14.1775 < 2.2e-16 *** --- Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1 I try to make a residuals plot with plot(moths.m6), the follow Warning occur:> plot(moths.m6)Warning message: NaNs produced in: sqrt(1 - hii) But the plots are produced. Then I try to plot a envelope function and its dont work.> envpois(moths.m6)Error in solve.default(t(X) %*% W %*% X) : singular matrix `a' in solve Without iterations this work well. What is the problem??? Thanks for all. The envelope function is: envpois <- function(fit.model) { #par(mfrow=c(1,1)) X <- model.matrix(fit.model) n <- nrow(X) p <- ncol(X) w <- fit.model$weights W <- diag(w) H <- solve(t(X)%*%W%*%X) H <- sqrt(W)%*%X%*%H%*%t(X)%*%sqrt(W) h <- diag(H) td <- resid(fit.model,type="deviance")/sqrt(1-h) e <- matrix(0,n,100) # for(i in 1:100){ nresp <- rpois(n, fitted(fit.model)) fit <- glm(nresp ~ X , family=poisson) w <- fit$weights W <- diag(w) H <- solve(t(X)%*%W%*%X) H <- sqrt(W)%*%X%*%H%*%t(X)%*%sqrt(W) h <- diag(H) e[,i] <- sort(resid(fit,type="deviance")/sqrt(1-h))} # e1 <- numeric(n) e2 <- numeric(n) # for(i in 1:n){ eo <- sort(e[i,]) e1[i] <- eo[5] e2[i] <- eo[95]} # med <- apply(e,1,mean) faixa <- range(td,e1,e2) #par(pty="s") qqnorm(td,xlab="Theoretical Quantiles", ylab="Std. deviance resid.", ylim=faixa, pch=1) par(new=T) # qqnorm(e1,axes=F,xlab="",ylab="",type="l",ylim=faixa,lty=1) par(new=T) qqnorm(e2,axes=F,xlab="",ylab="", type="l",ylim=faixa,lty=1) par(new=T) qqnorm(med,axes=F,xlab="", ylab="", type="l",ylim=faixa,lty=2) } -- They also surf who only stand on waves. -- | //|\\ [*****************************][*******************] || ( ? ? ) [Ronaldo Reis J?nior ][PentiumIII-600 ] | V [ESALQ/USP-Entomologia, CP-09 ][HD: 30 + 10 Gb ] || / l \ [13418-900 Piracicaba - SP ][RAM: 128 Mb ] | /(lin)\ [Fone: 19-429-4199 r.229 ][Video: SiS620-8Mb ] ||/(linux)\ [chrysopa at insecta.ufv.br ][Modem: Pctel-onboar] |/ (linux) \[ICQ#: 5692561 ][SO: CL 7.0 (2.2.19)] || ( x ) [*****************************][*******************] ||| _/ \_Powered by Conectiva Linux 7.0 D+:) | Lxuser#: 205366 -.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.- r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html Send "info", "help", or "[un]subscribe" (in the "body", not the subject !) To: r-help-request at stat.math.ethz.ch _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._