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

[zhe zhao]

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