[R] Problem to understand calculation of loglikelihood in the ramps package

[Bert Gunter]

You may get an answer here, but this appears to be something that you
should address to the package maintainer or package author.

Cheers,
Bert

On Mon, Oct 21, 2013 at 7:04 AM, zhe zhao <zzgt2006 at gmail.com> wrote:
> Dear R,
>
> We are trying to understand the calculation of loglikelihood in the ramps
> package. Our calculations do not agree with the package's. Can anyone
> explain why not?
>
> Here's an example using a small data set.
>
> # create small dataset
> library(ramps)
>
> data(NURE)  ## NURE is included in ramps package
> attach(NURE)
> data <- NURE[lat>41.9 & lon > -72.5 & !is.na(lat),]
> dim(data)  # 20 data points
>
> # ramps analysis
> set.seed(22)
>
> NURE.ctrl1 <- ramps.control(
>   iter = 10,
>   beta = param(0, "flat"),
>   sigma2.e = param(1, "invgamma", shape = 2.0, scale = 0.1, tuning = 0.75),
>   phi = param(10, "uniform", min = 0, max = 35, tuning = 0.50),
>   sigma2.z = param(1, "invgamma", shape = 2.0, scale = 0.1)
> )
>
> NURE.fit2 <- georamps(log(ppm) ~ 1,
>   correlation = corRExp(form = ~ easting + northing),
>
>   data = data,
>   control = NURE.ctrl1
> )
> rampsll <- NURE.fit2$loglik[1:4] # ramps loglikelihood for iterations 1:4
>
>
> # our calculation
> # We should be able to use the parameter values from a single MCMC iteration
> # to calculate the loglikelihood using dmvnorm
>
> library(mvtnorm)
> cor.exp <- function(x, range = 1, p = 1) # copied from corStruct.R in the
> ramps package
>
> {
>    if (range <= 0 || p <= 0)
>       stop("Exponential correlation parameter must be > 0")
>
>    if (p == 1) exp(x / (-1 * range))
>    else exp(-1 * (x / range)^p)
> }
>
> # Compute the covariance matrix
>
> sig.fn <- function(itno,data){
>   dist <- as.matrix ( dist ( data[,c("easting","northing")] ) )
>   npts <- nrow(data)
>   sig1 <- diag(npts)
>   for ( i in 1:(npts-1) )
>     for ( j in (i+1):npts )
>       sig1[j,i] <- sig1[i,j] <- cor.exp ( dist[i,j], range =
> NURE.fit2$params[itno,"phi"] )
>
>   Sigma <- sig1 * NURE.fit2$params[itno, "sigma2.z"] +
>     diag(npts) * NURE.fit2$params[itno, "sigma2.e"]
>   return(Sigma)
> }
>
> # Calculate the loglikelihood for a single MCMC iteration
> loglik.fn <- function(itno,data){
>
>   loglik <- dmvnorm ( x = log(data[,"ppm"]),
>     mean = rep ( NURE.fit2$params[itno,"(Intercept)"], nrow(data) ),
>     sigma = sig.fn(itno,data),
>     log = TRUE
>   )
>   return(loglik)
> }
>
> # Collect the loglikelihood from iterations 1:4
> myll <- c ( loglik.fn(1,data), loglik.fn(2,data), loglik.fn(3,data),
> loglik.fn(4,data) )
>
>
> But myll does not agree with rampsll.
> Can anyone tell us why not?
>
> Thanks very much for your help.
> Zhe
>
>         [[alternative HTML version deleted]]
>
> ______________________________________________
> R-help at 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.



-- 

Bert Gunter
Genentech Nonclinical Biostatistics

(650) 467-7374