R/llrwtfrp.R

llrwtfrp <-
function(x,y)
{
# Makes the weighted forward response and residual plots for loglinear regression.
# Workstation: need to activate a graphics
# device with command "X11()" or "motif()."

#
#   If q is changed, change the formula in the glm statement.
	q <- 5
# change formula to x[,1]+ ... + x[,q] with q
	out <- glm(y ~ x[, 1] + x[, 2] + x[, 3] +
			x[, 4] + x[,5], family = poisson)
	ESP <- x %*% out$coef[-1] + out$coef[1]
        Z <- y
        Z[y<1] <- Z[y<1] + 0.5
        out2<-lsfit(x,y=log(Z),wt=Z)
        #WRES <- sqrt(Z)*(log(Z) - x%*%out2$coef[-1] - out2$coef[1])
        WRES <- out2$res
        WFIT <- sqrt(Z)*log(Z) - WRES
        MWRES <- sqrt(Z)*(log(Z) - x%*%out$coef[-1] - out$coef[1])
        MWFIT <- sqrt(Z)*log(Z) - MWRES
        par(mfrow=c(2,2))
        plot(WFIT,sqrt(Z)*log(Z))
        abline(0,1)
        title("a) Weighted Forward Response Plot")
        plot(WFIT,WRES)
        title("b) Weighted Residual Plot")
        plot(MWFIT,sqrt(Z)*log(Z))
        abline(0,1)
        title("c) WFRP Based on MLE")
        plot(MWFIT,MWRES)
        title("d) WRP Based on MLE")
                                }
musto101/wilcox_R documentation built on May 23, 2019, 10:52 a.m.