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