Hi Simon,
thank you for your help. You suggest to use 'alpha_zero' and 'beta_zero' in
the 'L' function. But that's not exactly what I'm trying to do. Maybe it
helps if I show where the 'Likelihood_cov' model developed from. The basic
likelihood look like this:
### Likelihood function ###
Likelihood <- function(params, x, tx, T) {
r <- params[1]
alpha <- params[2]
s <- params[3]
beta <- params[4]
f <- function(x, tx, T)
{
g <- function(y)
(y + alpha)^(-( r + x))*(y + beta)^(-(s + 1))
integrate(g, tx, T)$value
}
integral <- mdply(data, f)
L <-
exp(lgamma(r+x)-lgamma(r)+r*(log(alpha)-log(alpha+T))-x*log(alpha+T)+s*(log(beta)-log(beta+T)))+exp(lgamma(r+x)-lgamma(r)+r*log(alpha)+log(s)+s*log(beta)+log(integral$V1))
f <- -sum(log(L))
return (f)
}
Parameters to be estimated are r, alpha, s and beta. Works fine so far. Now
I intent to incorporate the covariate 'IS' into this model. I do this by
facilitating a proportional hazard approach for the parameters 'alpha' and
'beta'. So alpha should be replaced by alpha_zero*exp(-gamma_1*IS) and
likewise beta shall become beta_zero*exp(-gamma_2*IS). So for that extended
model the parameters to be estimated are r, alpha_zero, s, beta_zero,
gamma_1 and gamma_2.
Therefore alpha and beta shall also be replaced in the L function by
alpha_zero*exp(-gamma_1*IS)
and beta_zero*exp(-gamma_2*IS) as well as with the integrate function. I
thought that would work by assigning:
data$alpha <- alpha_zero*exp(-gamma_1*IS)
data$beta <- beta_zero*exp(-gamma_2*IS)
as I did in the code I posted previously. But as you pointed out alpha and
beta do not exist when calling them in the L function. Unfortunately i do
not understand why they do not exist and how to fix this. Maybe you could
help me out here?
Thanks in advance and best regards,
Carlos
2013/9/3 Simon Zehnder <szehn...@uni-bonn.de>
> Hi Carlos,
>
> your problem is a wrong definition of your Likelihood function. You call
> symbols in the code (alpha, beta) which have no value assigned to. When L
> the long calculation in the last lines is assigned to L alpha and beta do
> not exist. The code below corrects it. But you have a problem with a
> divergent integral when calling integrate. A problem you can surely fix <as
> you know what your function is doing.
>
> Likelihood_cov <- function(params, x, tx, T, IS) {
> r <- params[1]
> alpha_zero <- params[2]
> s <- params[3]
> beta_zero <- params[4]
> gamma_1 <- params[5]
> gamma_2 <- params[6]
> data$alpha <- alpha_zero*exp(-gamma_1*IS)
> data$beta <- beta_zero*exp(-gamma_2*IS)
> f <- function(x, tx, T, alpha, beta)
> {
> g <- function(y)
> (y + alpha)^(-( r + x))*(y + beta)^(-(s + 1))
> integrate(g, tx, T)$value
> }
> integral <- mdply(data, f)
> L <-
>
> exp(lgamma(r+x)-lgamma(r)+r*(log(alpha_zero)-log(alpha_zero+T))-x*log(alpha_zero+T)+s*(log(beta_zero)-log(beta_zero+T)))+exp(lgamma(r+x)-lgamma(r)+r*log(alpha_zero)+log(s)+s*log(beta_zero)+log(integral$V1))
> f <- -sum(log(L))
> return (f)
> }
>
>
> Best
>
> Simon
>
>
> On Sep 3, 2013, at 1:28 PM, Carlos Nasher <carlos.nas...@googlemail.com>
> wrote:
>
> > Dear R helpers,
> >
> > I have problems to properly define a Likelihood function. Thanks to your
> > help my basic model is running quite well, but I have problems to get the
> > enhanced version (now incorporating covariates) running.
> >
> > Within my likelihood function I define a variable 'alpha'. When I want to
> > optimize the function I get the error message:
> >
> > 'Error in fn(ginv(par), ...) : object 'alpha' not found'
> >
> > I think it's actually not a problem with the optimization function
> (nmkb),
> > but with the Likelihood function itself. I do not understand why 'alpha'
> is
> > a missing object. 'alpha' should be part of the dataframe 'data' (as
> 'beta'
> > should be too), like 'x', 'tx', ''T. But it obviously isn't.
> >
> > Here's a minimum example which reproduces my problem:
> >
> > ##################################################################
> >
> > library(plyr)
> > library(dfoptim)
> >
> > ### Sample data ###
> > x <- c(3, 0, 2, 5, 1, 0, 0, 1, 0, 2)
> > tx <- c(24.57, 0.00, 26.86, 34.57, 2.14, 0.00, 0.00, 8.57, 0.00, 14.29)
> > T <- c(33.29, 30.71, 31.29, 34.57, 36.00, 35.43, 31.14, 33.86, 35.71,
> 35.86)
> > IS <- c(54.97, 13.97, 122.33, 110.84, 30.72, 14.96, 30.72, 20.74, 29.16,
> > 83.00)
> > data <- data.frame(x=x, tx=tx, T=T)
> > rm(x, tx, T)
> >
> > ### Likelihood function ###
> > Likelihood_cov <- function(params, x, tx, T, IS) {
> > r <- params[1]
> > alpha_zero <- params[2]
> > s <- params[3]
> > beta_zero <- params[4]
> > gamma_1 <- params[5]
> > gamma_2 <- params[6]
> > data$alpha <- alpha_zero*exp(-gamma_1*IS)
> > data$beta <- beta_zero*exp(-gamma_2*IS)
> > f <- function(x, tx, T, alpha, beta)
> > {
> > g <- function(y)
> > (y + alpha)^(-( r + x))*(y + beta)^(-(s + 1))
> > integrate(g, tx, T)$value
> > }
> > integral <- mdply(data, f)
> > L <-
> >
> exp(lgamma(r+x)-lgamma(r)+r*(log(alpha)-log(alpha+T))-x*log(alpha+T)+s*(log(beta)-log(beta+T)))+exp(lgamma(r+x)-lgamma(r)+r*log(alpha)+log(s)+s*log(beta)+log(integral$V1))
> > f <- -sum(log(L))
> > return (f)
> > }
> >
> > ### ML optimization ###
> > params <- c(0.2, 5, 0.2, 5, -0.02, -0.02)
> > fit <- nmkb(par=params, fn=Likelihood_cov, lower=c(0.0001, 0.0001,
> 0.0001,
> > 0.0001, -Inf, -Inf), upper=c(Inf, Inf, Inf, Inf, Inf, Inf), x=data$x,
> > tx=data$tx, T=data$T, IS=IS)
> >
> > ##################################################################
> >
> >
> > Maybe you could give me a hint were the flaw in my code is. Many thanks
> in
> > advance.
> > Carlos
> >
> >
> > -----------------------------------------------------------------
> > Carlos Nasher
> > Buchenstr. 12
> > 22299 Hamburg
> >
> > tel: +49 (0)40 67952962
> > mobil: +49 (0)175 9386725
> > mail: carlos.nas...@gmail.com
> >
> > [[alternative HTML version deleted]]
> >
> > ______________________________________________
> > 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.
>
>
--
-----------------------------------------------------------------
Carlos Nasher
Buchenstr. 12
22299 Hamburg
tel: +49 (0)40 67952962
mobil: +49 (0)175 9386725
mail: carlos.nas...@gmail.com
[[alternative HTML version deleted]]
______________________________________________
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.
