[R] significant difference between Gompertz hazard parameters?
Hello, all.
I have co-opted a number of functions that can be used to plot the
hazard/survival functions and associated density distribution for a Gompertz
mortality model, given known parameters. The Gompertz hazard model has been
shown to fit relatively well to the human adult lifespan. For example, if I
wanted to plot the hazard (i.e., mortality) functions:
pop1 <- function (t)
{
x=c(0.03286343, 0.04271132)
a3<-x[1]
b3<-x[2]
shift<-15 # only considering mortality after 15 years
h.t<-a3*exp(b3*(t-shift))
return<-h.t
}
pop2 <- function (t)
{
x=c(0.02207778, 0.04580059)
a3<-x[1]
b3<-x[2]
shift<-15 # only considering mortality after 15 years
h.t<-a3*exp(b3*(t-shift))
return<-h.t
}
ylab.name <- expression(paste(italic(h),"(",italic(a),")"))
plot(seq(15,80,1),pop1(seq(15,80,1)),type='l',ylab=ylab.name,xlab='Age
(years)',ylim=c(0,0.8))
lines(seq(15,80,1),pop2(seq(15,80,1)),lty=2)
How may I test for a significant difference in the hazard parameters that
define the mortality experience for these two populations? Thanks in
advance.
Regards,
Trey
-
Trey Batey---Anthropology Instructor
Division of Social Sciences
Mt. Hood Community College
Gresham, OR 97030
Alt. Email: trey.batey[at]mhcc[dot]edu
--
View this message in context:
http://r.789695.n4.nabble.com/significant-difference-between-Gompertz-hazard-parameters-tp4635018.html
Sent from the R help mailing list archive at Nabble.com.
__
R-help@r-project.org 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.
[R] summing two probability density functions from Gompertz hazard model
Hi, r-help members.
I have a question about summing two density distributions. I have two
samples from which I've estimated hazard parameters for a Gompertz mortality
model. With those parameters, I can calculate the PDF (survival function
times hazard function) of ages-at-death in a birth cohort subject to the
hazard function at each age. I'd like to combine these two density
functions and recalculate the resultant hazard parameters that describe the
combined group. Is there any way to do so? Unfortunately, the two
datasets, from which I estimated the initial parameters, are not in the same
format, so I can't just combine counts of individuals by age or age category
and calculate the parameters that way. I've included my functions for the
two distributions below. Any suggestions are welcome. Thanks.
--Trey
group1<- function (t){
x=c(0.05893007, 0.03339980)
a3<-x[1]
b3<-x[2]
shift<-15
S.t <- exp(a3/b3*(1-exp(b3*(t-shift
h.t <- a3*exp(b3*(t-shift))
return<-S.t*h.t}
plot(seq(15,80,1),group1(seq(15,80,1)),type='l',xlab='Age
years)',ylim=c(0,0.12)) # plot age-at-death distribution for group 1
group2<- function (t){
x=c(0.05920472, 0.01128975)
a3<-x[1]
b3<-x[2]
shift<-15
S.t <- exp(a3/b3*(1-exp(b3*(t-shift
h.t <- a3*exp(b3*(t-shift))
return<-S.t*h.t}
plot(seq(15,80,1),group2(seq(15,80,1)),type='l',xlab='Age
(years)',ylim=c(0,0.12)) # plot age-at-death distribution for group 2
-
Trey Batey---Anthropology Instructor
Division of Social Sciences
Mt. Hood Community College
Gresham, OR 97030
Alt. Email: trey.batey[at]mhcc[dot]edu
--
View this message in context:
http://r.789695.n4.nabble.com/summing-two-probability-density-functions-from-Gompertz-hazard-model-tp4581793p4581793.html
Sent from the R help mailing list archive at Nabble.com.
__
R-help@r-project.org 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.
[R] using a loop with an integration
Hi, all.
I've written a function that returns the survival function for a Gompertz
mortality model. I've specified the two model parameters. Using a simple
integration, I can calculate the life expectancy at any age. Is there a way
I can use a loop with the integration that will quickly return life
expectancy over a range of ages, say 0 to 80, so that I don't have to
manually type in the age in which I'm interested? Please see the code
below. Thanks so much.
--Trey
hk.bothsex_Gompsurv <- function (t)
{
x=c(0.02342671, 0.05837508)
a3<-x[1]
b3<-x[2]
shift<-15
S.t<-exp(a3/b3*(1-exp(b3*(t-shift
return<-S.t
}
integrate(hk.bothsex_Gompsurv,0,Inf)$value/hk.bothsex_Gompsurv(0) # life
expectancy at birth (change lower limit of integral and corresponding "t" in
denominator to calculate life expectancy at any age
-
Trey Batey---Anthropology Instructor
Division of Social Sciences
Mt. Hood Community College
Gresham, OR 97030
Alt. Email: trey.batey[at]mhcc[dot]edu
--
View this message in context:
http://r.789695.n4.nabble.com/using-a-loop-with-an-integration-tp4578752p4578752.html
Sent from the R help mailing list archive at Nabble.com.
__
R-help@r-project.org 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.
[R] Gompertz-Makeham hazard models---test for significant difference
Hi, all.
I'm working with published paleodemographic data (counts of skeletons that
have been assigned into an age-range category, based upon observed
morphological characteristics). For example, the following is the age
distribution from a prehistoric cemetery in Egypt:
naga <-
structure(c(15,20,35,50,20,35,50,Inf,46,43,26,14),.Dim=as.integer(c(4,3)),.Dimnames=list(NULL,c("col1","col2","col3")))
> naga
col1 col2 col3
[1,] 15 20 46
[2,] 20 35 43
[3,] 35 50 26
[4,] 50 Inf 14
So above, the first two columns are the lower and upper limits of the
age-range categories (typically used by paleodemographers), and the third
column is the number of skeletons assigned to each group. I've co-opted
some R script for fitting a Gompertz-Makeham hazard model to this data...
##
GM.naga <- function(x,deaths=naga)
{
a2=x[1]
a3=x[2]
b3=x[3]
shift<-15
nrow<-NROW(deaths)
S.t<-function(t)
{
return(exp(-a2*(t-shift)+a3/b3*(1-exp(b3*(t-shift)
}
d<-S.t(deaths[1:nrow,1])-S.t(deaths[1:nrow,2])
obs<-deaths[,3]
lnlk<-as.numeric(crossprod(obs,log(d)))
return(lnlk)
}
optim(c(0.001, 0.01, 0.1),GM.naga,control=list(fnscale=-1))
And the output:
$par
[1] 0.05486053 0.31290971 -1.87744033
$value
[1] -168.3748
$counts
function gradient
174 NA
$convergence
[1] 0
$message
NULL
Let's say that I have data from another cemetery, for which I would estimate
another set of hazard model parameters. How would I go about comparing the
two to see if the resulting parameters for each (and their resulting
survival, hazard, and density curves) are significantly different? Also,
what if I wanted to include data from even more cemeteries and compare many
sets of estimated hazard parameters? Below, I've included a the
data/results for the another cemetery that I'd like to compare to the first.
Any suggestions are welcome. Thanks so much.
--Trey
#
naqada <-
structure(c(15,20,25,50,19,24,49,Inf,26,45,219,30),.Dim=as.integer(c(4,3)),.Dimnames=list(NULL,c("col1","col2","col3")))
GM.naqada <- function(x,deaths=naqada)
{
a2=x[1]
a3=x[2]
b3=x[3]
shift<-15
nrow<-NROW(deaths)
S.t<-function(t)
{
return(exp(-a2*(t-shift)+a3/b3*(1-exp(b3*(t-shift)
}
d<-S.t(deaths[1:nrow,1])-S.t(deaths[1:nrow,2])
obs<-deaths[,3]
lnlk<-as.numeric(crossprod(obs,log(d)))
return(lnlk)
}
optim(c(0.001, 0.01, 0.1),GM.naqada,control=list(fnscale=-1))
And the output:
$par
[1] 1.544739e-08 1.862670e-02 6.372117e-02
$value
[1] -331.1865
$counts
function gradient
276 NA
$convergence
[1] 0
$message
NULL
Trey Batey---Anthropology Instructor
Division of Social Sciences
Mt. Hood Community College
Gresham, OR 97030
trey.batey[at]mhcc[dot]edu
--
View this message in context:
http://r.789695.n4.nabble.com/Gompertz-Makeham-hazard-models-test-for-significant-difference-tp4560550p4560550.html
Sent from the R help mailing list archive at Nabble.com.
__
R-help@r-project.org 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.
