pneumo = transform(pneumo, let=log(exposure.time))
(fit = vglm(cbind(normal, mild, severe) ~ let,
cumulative(parallel=TRUE, reverse=TRUE), pneumo))
fit@y # These are sample proportions
weights(fit, type="prior", matrix=FALSE) # Number of observations
# Look at the working residuals
n = nrow(model.matrix(fit, type="lm"))
M = ncol(predict(fit))
temp = weights(fit, type="working", deriv=TRUE)
wz = m2adefault(temp$weights, M=M) # In array format
wzinv = array(apply(wz, 3, solve), c(M,M,n))
wresid = matrix(NA, n, M) # working residuals
for(i in 1:n)
wresid[i,] = wzinv[,,i,drop=TRUE] %*% temp$deriv[i,]
max(abs(c(resid(fit, type="w")) - c(wresid))) # Should be 0
z = predict(fit) + wresid # Adjusted dependent vector
z
Run the code above in your browser using DataLab