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.
______________________________________________
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.
