earg = list(cutpoint=2)
polf("p", earg=earg, short=FALSE)
polf("p", earg=earg, tag=TRUE)
p = seq(0.01, 0.99, by=0.01)
y = polf(p, earg=earg)
y. = polf(p, earg=earg, deriv=1)
max(abs(polf(y, earg=earg, inv=TRUE) - p)) # Should be 0
par(mfrow=c(2,1), las=1)
plot(p, y, type="l", col="blue", main="polf()")
abline(h=0, v=0.5, col="red", lty="dashed")
plot(p, y., type="l", col="blue",
main="(Reciprocal of) first POLF derivative")
# Rutherford and Geiger data
ruge = data.frame(yy = rep(0:14,
times=c(57,203,383,525,532,408,273,139,45,27,10,4,0,1,1)))
with(ruge, length(yy)) # 2608 1/8-minute intervals
cutpoint = 5
ruge = transform(ruge, yy01 = ifelse(yy <= cutpoint, 0, 1))
earg = list(cutpoint=cutpoint)
fit = vglm(yy01 ~ 1, binomialff(link="polf", earg=earg), ruge)
coef(fit, matrix=TRUE)
exp(coef(fit))
# Another example
pdat = data.frame(x2 = sort(runif(nn <- 1000)))
pdat = transform(pdat, x3 = runif(nn))
pdat = transform(pdat, mymu = exp( 3 + 1 * x2 - 2 * x3))
pdat = transform(pdat, y1 = rpois(nn, lambda=mymu))
cutpoints = c(-Inf, 10, 20, Inf)
pdat = transform(pdat, cuty = Cut(y1, breaks=cutpoints))
with(pdat, plot(x2, x3, col=cuty, pch=as.character(cuty)))
with(pdat, table(cuty) / sum(table(cuty)))
fit = vglm(cuty ~ x2 + x3, fam = cumulative(link="polf",
reverse=TRUE, parallel=TRUE, intercept.apply=TRUE,
mv=TRUE, earg=list(cutpoint=cutpoints[2:3])),
pdat, trace=TRUE)
head(fit@y)
head(fitted(fit))
head(predict(fit))
coef(fit)
coef(fit, matrix=TRUE)
constraints(fit)
fit@misc$earg
Run the code above in your browser using DataLab