[R] Gompertz-Makeham hazard models---test for significant difference

[piltdownpunk]

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.