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 = c(57,203,383,525,532,408,273,139,45,27,10,4,0,1,1)
y = 0:14
yy = rep(y, times=ruge)
length(yy) # 2608 1/8-minute intervals
cutpoint = 5
yy01 = ifelse(yy <= cutpoint, 0, 1)
earg = list(cutpoint=cutpoint)
fit = vglm(yy01 ~ 1, binomialff(link="polf", earg=earg))
coef(fit, matrix=TRUE)
exp(coef(fit))
# Another example
nn = 1000
x2 = sort(runif(nn))
x3 = runif(nn)
mymu = exp( 3 + 1 * x2 - 2 * x3)
y1 = rpois(nn, lambda=mymu)
cutpoints = c(-Inf, 10, 20, Inf)
cuty = Cut(y1, breaks=cutpoints)
plot(x2, x3, col=cuty, pch=as.character(cuty))
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])),
trace=TRUE)
fit@y[1:5,]
fitted(fit)[1:5,]
predict(fit)[1:5,]
coef(fit)
coef(fit, matrix=TRUE)
constraints(fit)
fit@misc$earg
Run the code above in your browser using DataLab