[R] error in plot residuals in a glm with iterations.

[Ronaldo Reis Jr.]

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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._