R/poireg.R

In musto101/wilcox_R: What the package does (short line)

poireg <-
function(x,y,xout=FALSE,outfun=outpro,plotit=FALSE,xlab="X",ylab="Y",
varfun=var,YHAT=FALSE,STAND=TRUE,...){
#
# Perform Poisson regression.
# The predictors are assumed to be stored in the n by p matrix x.
# The y values are typically count data (integers).
#
# xout=T will remove outliers from among the x values and then fit
# the regression line.
#  Default:
# One predictor, a mad-median rule is used.
# With more than one, projection method is used.
#
# outfun=out will use MVE method
#
xy=elimna(cbind(x,y))
x<-as.matrix(x)
x=xy[,1:ncol(x)]
y=xy[,ncol(xy)]
x<-as.matrix(x)
if(xout){
flag<-outfun(x,...)$keep
x<-x[flag,]
y<-y[flag]
x<-as.matrix(x)
}
temp=glm(formula=y~x,family=poisson)
init=summary(temp)
yhat=temp$coef[1]
for(j in 1:ncol(x)){
j1=j+1
yhat=yhat+temp$coef[j1]*x[,j]
}
yhat=exp(yhat)
if(plotit){
x=as.matrix(x)
if(ncol(x)>1)stop("Cannot plot with more than one predictor")
plot(x,y,xlab=xlab,ylab=ylab)
#points(x,yhat,pch=".")
xord=order(x)
lines(x[xord],yhat[xord])
init$coef
}
ex=varfun(yhat)/varfun(y)
str=sqrt(ex)
hatv=NULL
if(YHAT)hatv=yhat
list(results=init,Explanatory.Power=ex,Strength.Assoc=str,yhat=hatv)
}