R help - Oct 2013 - Problem to understand calculation of loglikelihood in the ramps package

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]]

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