Brian Ripley wrote: > > BTW, your example cannot be pasted in as 'sdat' self-references. It could > be fixed, but I gave up at that point. > Oh dear, I'm very sorry. I forgot to run rm(list=ls(all=TRUE)) before testing. The corrected code is:
#Data: S <- c(0, 0, 0, 0, 28, 0, 1, 7, 0, 0, 39, 2, 0, 0, 40, 0, 0, 0, 0, 0, 0, 15, 0, 0, 3, 0, 0, 0, 2, 0, 3, 0, 30, 0, 20, 0, 1, 0, 0, 1, 21, 0, 0, 4, 14, 0, 0, 0, 1, 1, 1, 0, 0, 0, 1, 1, 0, 3, 0, 2, 5, 0, 0, 0, 0, 0, 0, 0, 25, 0, 5, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 10, 0, 0, 0, 0, 0, 9, 1, 1, 1, 1, 0, 0, 3, 0, 27, 7, 0, 0, 0, 0, 0, 1, 1, 4, 2, 1, 2, 4, 1, 4, 6, 12, 4, 6, 3, 4, 0, 4, 0, 6, 1, 3, 1, 4, 4, 1, 1, 2, 1, 6, 1, 0, 0, 1, 0, 1, 0, 0, 6, 0, 0, 0, 0, 2, 2, 3, 1, 6, 2, 2, 1, 1, 4, 4, 3, 3, 7, 1, 3, 5, 6, 4, 0, 1, 4, 2, 0, 1, 1, 1, 0, 1, 1, 1, 0, 0, 1, 0, 4, 0, 0, 1, 0, 2, 0, 0, 1, 2, 0, 4, 0, 2, 0, 3, 0, 2, 2, 0, 4, 0, 2, 0, 1, 2, 3, 0, 0, 0, 0, 3, 1, 0, 0, 0, 3, 2, 1, 2, 1, 1, 2, 0, 0, 3, 3, 1, 0, 1, 1, 2, 0, 1, 0, 1, 0, 1, 3, 2, 0, 0, 2, 1, 0, 1, 0, 0, 0, 0, 0, 0, 4, 2, 4, 0, 2, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 1, 3, 2, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 4, 1, 0, 0, 1) M <- 620+c(0,cumsum(S[-328])) sdat <- data.frame(S,M) #S is the number of N individuals that irreversibly change state in a time #interval t. The data provided are a subset of the full data set. #M is the cumulative sum of individuals that have changed state up to t-1. #Assume that the rate of state change is constant (S ~ kM), but the #distribution of S is clustered. #N can be estimated by fitting: qpglm <- glm(S ~ M, family = quasipoisson(link = identity), sdat) summary(qpglm) N.est <- -coef(qpglm)[1]/coef(qpglm)[2] N.est #i.e. x-intercept is minus intercept/slope #To estimate confidence limits on N.est, fit models without intercept to #N.est - M + x, where x is an integer. The model will have the lowest deviance #when x = 0. x <- 0 Mdif <- N.est - M + x qpglm2 <- glm(S ~ -1 + Mdif, family = quasipoisson(link = identity), sdat) summary(qpglm2) #Use analysis of deviance to estimate confidence limits on N.est: disp <- summary(qpglm)$dispersion dfres <- df.residual(qpglm) dev.res <- deviance(qpglm) #From MASS4, p. 210, assume that changes in deviance scaled by #dispersion as |x| increases have an F distribution dev.crit <- dev.res+qf(0.95,1,dfres)*disp dev.crit #values of x for which the deviance = dev.crit give approximate 95% confidence limits #on N.est. #The error occurs when x <= -91.7: x <- -91.7 sdat$Mdif <- N.est - sdat$M + x strt <- coef(glm(S ~ -1 + I(Mdif+1), family = quasipoisson(link = identity), sdat)) qpglm2 <- glm(S ~ -1 + Mdif, family = quasipoisson(link = identity), start=strt, sdat) ______________________________________________ [email protected] mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
