[R] significant difference between Gompertz hazard parameters?

[T-Bone]

Hi, David.

Thank you for the suggestions.  If it makes a difference, I've included more
information on the data and how the Gompertz model was fit to each sample
below.  I'm sure my methods are inefficient, and I should utilizing R's
process of vectorization.  Thanks, again.

--Trey


David Winsemius wrote
> 
> If you supply some data, possibly generated through simulations, you can
> test for differences in fit using parametric fits.
> 

The data come from two assemblages of skeletons (each from a different
cemetery), a type of dataset common for anthropologists/archaeologists.  The
model was fit to counts of skeletons assigned to one of four age-range
categories (based on observed skeletal morphology).  For example:

pop1 <-
structure(c(15,20,35,50,20,35,50,Inf,50,166,88,24),.Dim=as.integer(c(4,3)),.Dimnames=list(NULL,c("age.start","age.end","dead")))    
## age-structure of the cemetery

pop1_Gomp <- function(x,deaths=pop1)
	{
		a3=x[1]
		b3=x[2]
		shift<-15
		nrow<-NROW(deaths)
		
		S.t<-function(t)
			{
				return(exp(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.01, 0.01),pop1_Gomp,control=list(fnscale=-1))

## output...
$par
[1] 0.03286343    0.04271132

$value
[1] -386.2922

$counts
function gradient 
     129       NA 

$convergence
[1] 0

$message
NULL
##################

In the original message, "pop2," representing a separate cemetery, had the
following age structure:

pop2 <-
structure(c(15,20,35,50,20,35,50,Inf,46,310,188,101),.Dim=as.integer(c(4,3)),.Dimnames=list(NULL,c("age.start","age.end","dead")))

...and the Gompertz parameters estimated using the function:

pop2_Gomp <- function(x,deaths=pop2)	
	{
		a3=x[1]
		b3=x[2]
		shift<-15
		nrow<-NROW(deaths)
		
		S.t<-function(t)
			{
				return(exp(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.01, 0.01),pop2_Gomp,control=list(fnscale=-1))
#######################
$par
[1] 0.02207778     0.04580059

$value
[1] -781.7549

$counts
function gradient 
      71       NA 

$convergence
[1] 0

$message
NULL


-----
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-tp4635018p4635736.html
Sent from the R help mailing list archive at Nabble.com.