likelihood.smooth.spline: Calculate the log likelihood of a smoothing spline given the...

In HenrikBengtsson/aroma.light-BioC_release: Light-Weight Methods for Normalization and Visualization of Microarray Data using Only Basic R Data Types

Description

Calculate the (log) likelihood of a spline given the data used to fit the spline, g. The likelihood consists of two main parts: 1) (weighted) residuals sum of squares, and 2) a penalty term. The penalty term consists of a smoothing parameter lambda and a roughness measure of the spline J(g) = \int g''(t) dt. Hence, the overall log likelihood is

\log L(g|x) = (y-g(x))'W(y-g(x)) + λ J(g)

In addition to the overall likelihood, all its seperate components are also returned.

Note: when fitting a smooth spline with (x,y) values where the x's are not unique, smooth.spline will replace such (x,y)'s with a new pair (x,y') where y' is a reweighted average on the original y's. It is important to be aware of this. In such cases, the resulting smooth.spline object does not contain all (x,y)'s and therefore this function will not calculate the weighted residuals sum of square on the original data set, but on the data set with unique x's. See examples below how to calculate the likelihood for the spline with the original data.

Usage

## S3 method for class 'smooth.spline'
likelihood(object, x=NULL, y=NULL, w=NULL, base=exp(1),
  rel.tol=.Machine$double.eps^(1/8), ...)

Arguments

object	The smooth.spline object.
x, y	The x and y values for which the (weighted) likelihood will be calculated. If x is of type xy.coords any value of argument y will be omitted. If x==NULL, the x and y values of the smoothing spline will be used.
w	The weights for which the (weighted) likelihood will be calculated. If NULL, weights equal to one are assumed.
base	The base of the logarithm of the likelihood. If NULL, the non-logged likelihood is returned.
rel.tol	The relative tolerance used in the call to integrate.
...	Not used.

Details

The roughness penalty for the smoothing spline, g, fitted from data in the interval [a,b] is defined as

J(g) = \int_a^b g''(t) dt

which is the same as

J(g) = g'(b) - g'(a)

The latter is calculated internally by using predict.smooth.spline.

Value

Returns the overall (log) likelihood of class SmoothSplineLikelihood, a class with the following attributes:

wrss	the (weighted) residual sum of square
penalty	the penalty which is equal to -lambda*roughness.
lambda	the smoothing parameter
roughness	the value of the roughness functional given the specific smoothing spline and the range of data

Author(s)

Henrik Bengtsson

See Also

smooth.spline and robustSmoothSpline().

Examples

# Define f(x)
f <- expression(0.1*x^4 + 1*x^3 + 2*x^2 + x + 10*sin(2*x))

# Simulate data from this function in the range [a,b]
a <- -2; b <- 5
x <- seq(a, b, length=3000)
y <- eval(f)

# Add some noise to the data
y <- y + rnorm(length(y), 0, 10)

# Plot the function and its second derivative
plot(x,y, type="l", lwd=4)

# Fit a cubic smoothing spline and plot it
g <- smooth.spline(x,y, df=16)
lines(g, col="yellow", lwd=2, lty=2)

# Calculating the (log) likelihood of the fitted spline
l <- likelihood(g)

cat("Log likelihood with unique x values:\n")
print(l)

# Note that this is not the same as the log likelihood of the
# data on the fitted spline iff the x values are non-unique
x[1:5] <- x[1]  # Non-unique x values
g <- smooth.spline(x,y, df=16)
l <- likelihood(g)

cat("\nLog likelihood of the *spline* data set:\n");
print(l)

# In cases with non unique x values one has to proceed as
# below if one want to get the log likelihood for the original
# data.
l <- likelihood(g, x=x, y=y)
cat("\nLog likelihood of the *original* data set:\n");
print(l)

HenrikBengtsson/aroma.light-BioC_release documentation built on May 7, 2019, 1:55 a.m.