zibreg: Regression model for binomial data with unkown group of...

View source: R/zib.R

zibregR Documentation

Regression model for binomial data with unkown group of immortals

Description

Regression model for binomial data with unkown group of immortals (zero-inflated binomial regression)

Usage

zibreg(
  formula,
  formula.p = ~1,
  data,
  family = stats::binomial(),
  offset = NULL,
  start,
  var = "hessian",
  ...
)

Arguments

formula

Formula specifying

formula.p

Formula for model of disease prevalence

data

data frame

family

Distribution family (see the help page family)

offset

Optional offset

start

Optional starting values

var

Type of variance (robust, expected, hessian, outer)

...

Additional arguments to lower level functions

Author(s)

Klaus K. Holst

Examples


## Simulation
n <- 2e3
x <- runif(n,0,20)
age <- runif(n,10,30)
z0 <- rnorm(n,mean=-1+0.05*age)
z <- cut(z0,breaks=c(-Inf,-1,0,1,Inf))
p0 <- lava:::expit(model.matrix(~z+age) %*% c(-.4, -.4, 0.2, 2, -0.05))
y <- (runif(n)<lava:::tigol(-1+0.25*x-0*age))*1
u <- runif(n)<p0
y[u==0] <- 0
d <- data.frame(y=y,x=x,u=u*1,z=z,age=age)
head(d)

## Estimation
e0 <- zibreg(y~x*z,~1+z+age,data=d)
e <- zibreg(y~x,~1+z+age,data=d)
compare(e,e0)
e
PD(e0,intercept=c(1,3),slope=c(2,6))

B <- rbind(c(1,0,0,0,20),
           c(1,1,0,0,20),
           c(1,0,1,0,20),
           c(1,0,0,1,20))
prev <- summary(e,pr.contrast=B)$prevalence

x <- seq(0,100,length.out=100)
newdata <- expand.grid(x=x,age=20,z=levels(d$z))
fit <- predict(e,newdata=newdata)
plot(0,0,type="n",xlim=c(0,101),ylim=c(0,1),xlab="x",ylab="Probability(Event)")
count <- 0
for (i in levels(newdata$z)) {
  count <- count+1
  lines(x,fit[which(newdata$z==i)],col="darkblue",lty=count)
}
abline(h=prev[3:4,1],lty=3:4,col="gray")
abline(h=prev[3:4,2],lty=3:4,col="lightgray")
abline(h=prev[3:4,3],lty=3:4,col="lightgray")
legend("topleft",levels(d$z),col="darkblue",lty=seq_len(length(levels(d$z))))

kkholst/lava documentation built on Feb. 22, 2024, 4:07 p.m.