From: Camarda, Carlo Giovanni <Camarda_at_demogr.mpg.de>

Date: Tue 28 Feb 2006 - 07:23:04 EST

# fitting the first dataset

fit <- glm(formula = deaths.sim ~ age + offset(my.offset),

# plotting log-hazard

plot(age, log(hazard.sim), cex=1.2)

# plotting residuals

plot(age, res, cex=1.2, lwd=2, ylim=range(res, res1)) points(age, res1, col="red", cex=1.2, lwd=2, pch=2) all(round(res,5)==round(res1,5)) ## FALSE

# looking at the ratio

ratio <- (res/res1)[1]

ratio

all(round(res,10)==round(res1*ratio,10))

# here the scatterplot resid-ratios vs multiplier-constant

R-help@stat.math.ethz.ch mailing list

https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html Received on Tue Feb 28 07:33:47 2006

Date: Tue 28 Feb 2006 - 07:23:04 EST

Dear R-users,

I would like to show you a simple example that gives an overview of one
of my current issue.

Although my working setting implies a different parametric model (which
cannot be framed in the glm), I guess that what I'll get from the
following example it would help for the next steps.
Anyway here it is.

Firstly I simulated from a series of exposures, a series of deaths
(given a model, Gompertz, and a probability distribution, Poisson). Then
a multiply both deaths and exposures by a constant. Finally I fitted
with glm the simulated data sets and I compared the deviance residuals.
They're different by a certain constant, but I could not get a
meaningful relationship between the used constants (a scatter plot of
them is given at the end).

The following example tried to be as completed as possible and I hope
that it will be clearer than my previous words.
Thanks,

Carlo Giovanni Camarda

#### simulation procedure

age <- 50:100 # time range

# parameters

alpha.sim <- 0.00005

beta.sim <- 0.1

# simulated hazard from a Gompertz

hazard.sim <- alpha.sim*(exp(beta.sim*age))

# first dataset

exposures <- c(4748, 4461, 4473, 4326, 4259, 4219, 4049, 3965, 3801,
3818, 3670, 3521, 3537,

3482, 3347, 3180, 3100, 2890, 2755, 2622, 2530, 2502, 2293, 2216, 1986, 1875,

1693, 1550, 1395, 1295, 1104, 952, 792, 755, 726, 608, 523, 419, 338,

246, 205, 148, 112, 75, 52, 32, 23, 15, 7, 6, 3) # just
as example

deaths.sim <- rpois(length(age), exposures*hazard.sim) # simulating
deaths from a poisson distribution (Brillinger, 1986)
my.offset <- log(exposures) # offset for the poisson regression

# new dataset: decupleing the sample size

deaths.sim1 <- deaths.sim * 10 exposures1 <- exposures * 10 my.offset1 <- log(exposures1)

# fitting the first dataset

fit <- glm(formula = deaths.sim ~ age + offset(my.offset),

family = poisson(link = "log"))
res <- residuals(fit, type="deviance")

# fitting the new dataset

fit1 <- glm(formula = deaths.sim1 ~ age + offset(my.offset1),

family = poisson(link = "log"))
res1 <- residuals(fit1, type="deviance")

# estimating hazard

hazard.act <- deaths.sim/exposures # == deaths.sim1/exposures1 hazard.fit <- predict(fit, type="response") / exposures hazard.fit1 <- predict(fit1, type="response") / exposures1

# plotting log-hazard

plot(age, log(hazard.sim), cex=1.2)

lines(age, log(hazard.act), col="red", lwd=2) lines(age, log(hazard.fit), col="blue", lwd=2, lty=2) lines(age, log(hazard.fit1), col="green", lwd=2, lty=3)all(round(hazard.fit,5)==round(hazard.fit1,5)) ## TRUE

# plotting residuals

plot(age, res, cex=1.2, lwd=2, ylim=range(res, res1)) points(age, res1, col="red", cex=1.2, lwd=2, pch=2) all(round(res,5)==round(res1,5)) ## FALSE

# looking at the ratio

ratio <- (res/res1)[1]

ratio

all(round(res,10)==round(res1*ratio,10))

# here the scatterplot resid-ratios vs multiplier-constant

const.a <- c(1:16) # constants which I tried to multiply by const.b <- c(1,0.7071068,0.5773503,0.5,0.4472136,0.4082483,0.3779645,0.3535534, #ratio between the deviance residuals

0.3333333,0.3162278,0.3015113,0.2886751,0.2773501,0.2672612,0.2581989,0.
25)

plot(const.a, const.b, pch=16)

+++++

This mail has been sent through the MPI for Demographic Rese...{{dropped}}

R-help@stat.math.ethz.ch mailing list

https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html Received on Tue Feb 28 07:33:47 2006

*
This archive was generated by hypermail 2.1.8
: Fri 03 Mar 2006 - 03:42:46 EST
*