Skip to content

significant difference between Gompertz hazard parameters?

2 messages · David Winsemius, T-Bone

David Winsemius
#
On Jul 1, 2012, at 12:04 AM, piltdownpunk wrote:

            

      
      
      
    

        
You cannot test for differences in pre-specified parameters. These are  
by definition, "different". If you supply some data, possibly  
generated through simulations, you can test for differences in fit  
using parametric fits. The survival package offers facilities for  
fitting,  and in the Archives you can find several responses from  
Terry Therneau to questions about fitting data to Gompertz or Gompertz- 
Makeham distributions.

            

      
      
      
    

        

  
    

      
      
    

  


  


    

    
6 days later

    

      
      

        


  

        

      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

            

      
      
      
    

        
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.