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.
