Scott's rule of thumb

Description

Calculates Scott's rule of thumb bandwidth vector.

Usage

bandwidth.scott(x, kernel = "biweight", product = TRUE)

Arguments

Details

The default bandwidth vector is computed by Scott's rule of thumb for the Gaussian kernel and adapted to the chosen kernel function.

Value

d x 1 bandwidth vector used for calculation

Author(s)

Marlene Mueller

References

Scott, D.W. (1992). Multivariate Density Estimation: Theory, Practice, and Visualization. New York, Chichester: Wiley.

See Also

kernel.function, kde

Examples

  ## two-dimensional data
  n <- 1000
  u <- runif(n)
  thresh <- 0.4
  x1 <- rnorm(n)*(u<thresh) +rnorm(n,mean=3)*(u>=thresh)
  x2 <- rnorm(n)*(u<thresh) +rnorm(n,mean=9)*(u>=thresh)
  bandwidth.scott( cbind(x1,x2) )

Kernel convolution

Description

Calculates the convolution of data with a kernel function.

Usage

convol(x, h = 1, grid = NULL, y = 1, w = 1, p = 2, q = 2,
       product = TRUE, sort = TRUE)

Arguments

Details

The kernel convolution which is calculated is \sum_i K_h(x_i - grid_{j})\,y_i\,w_{ij} for i=1,...,n and j=1,...,m. The kernel function is determined by the kernel parameters p and q, see kernel.function. The default kernel is the biweight (quartic) kernel function. Note that the DLL requires the data matrix to be sorted by its first column.

Value

m x c matrix

Author(s)

Marlene Mueller

See Also

kernel.function, kde, kreg

Examples

  n <- 100
  x <- rnorm(n)
  convol(x,h=0.8,grid=-3:3)/n  ## estimates density of x at points -3:3

Create a grid for kernel estimation

Description

Helps to define a grid for kernel denity or regression estimates (univariate or multivariate).

Usage

create.grid(grid.list, sort=TRUE)

Arguments

Details

This function allows easily to define grids for the "gplm" package. If the data are d-dimensional and the grid vector lengths are n1, ... nd, then the output is a (n1*...*nd) x d matrix with each row corresponding to one d-dimensional data point at which the function estimate is to be calculated.

Value

m x d grid matrix

Author(s)

Marlene Mueller

See Also

expand.grid, kde, kreg

Examples

  v1 <- 1:5
  v2 <- 3:1
  grid <- create.grid(list(v1,v2))

  x <- matrix(rnorm(60),30,2)
  v1 <- seq(min(x[,1]),max(x[,1]),length=10)
  v2 <- seq(min(x[,2]),max(x[,2]),length=5)
  grid <- create.grid(list(v1,v2))

Link function for GLM

Description

Defines the link function for a GLM.

Usage

glm.inverse.link(mu, family="gaussian", link="identity", k=1)

Arguments

Value

n x 1, vector eta (predictors)

Author(s)

Marlene Mueller

See Also

glm.ll, glm.lld, glm.link

Examples

  glm.inverse.link(c(0.25,0.5), family="bernoulli", link="logit")

(Inverse) Link function for GLM

Description

Defines the inverse link function for a GLM.

Usage

glm.link(eta, family="gaussian", link="identity", k=1)

Arguments

Value

n x 1, vector mu (responses)

Author(s)

Marlene Mueller

See Also

glm.ll, glm.lld, glm.inverse.link

Examples

  glm.link(c(-1,2), family="bernoulli", link="logit")

Log-likelihood for GLM

Description

Calculates the log-likelihood function of a GLM. Currently only the gaussian and the bernoulli family are implemented.

Usage

glm.ll(mu, y, phi=1, family="gaussian", k=1)

Arguments

Details

Implemented are the "gaussian" family (with links "identity" and "log"), the "bernoulli" family (with links "logit" and "probit"), the "gamma" family (with link "inverse"), the "poisson" family (with link "log"), the "inverse.gaussian" family (with link "inverse.squared") and the "negative.binomial" (with its canonical "log" type link).

The default value k=1 leads to the geometric distribution (as a special case of the negative binomial).

Value

log-likelihood value

Author(s)

Marlene Mueller

See Also

glm.lld, glm.link

Examples

  glm.ll(rep(0.4,2), c(0,1), family="bernoulli")

Log-likelihood derivatives for GLM

Description

Computes first and second derivatives of the individual log-likelihood with respect to the linear predictor. Currently only the gaussian (with identity link) and the bernoulli family (with logit and probit links) are implemented.

Usage

glm.lld(eta, y, family="gaussian", link="identity", k=1)

Arguments

Details

See details for glm.ll.

Value

List with components:

Author(s)

Marlene Mueller

See Also

glm.ll, glm.link

Examples

  glm.lld(c(-1,2), c(0,1), family="bernoulli", link="logit")

Backfitting for an additive model using kernel regression

Description

Implements kernel-based backfitting in an additive model, optional with a partial linear term.

Usage

kbackfit(t, y, h, x = NULL, grid = NULL, weights.conv = 1,
           offset = 0, method = "generic",
           max.iter = 50, eps.conv = 1e-04, m.start = NULL,
           kernel = "biweight")

Arguments

Value

List with components:

Author(s)

Marlene Mueller

See Also

kernel.function, kreg

Kernel density estimation

Description

Calculates a kernel density estimate (univariate or multivariate).

Usage

kde(x, bandwidth = NULL, grid = TRUE, kernel = "biweight",
    product = TRUE, sort = TRUE)

Arguments

Details

The kernel density estimator is calculated as \frac{1}{n} \sum_i K_h(x_i - grid_{j}) for i=1,...,n and j=1,...,m. The default bandwidth vector is computed by Scott's rule of thumb (adapted to the chosen kernel function).

Value

List with components:

Author(s)

Marlene Mueller

See Also

kernel.function, convol, kreg

Examples

  n <- 1000
  x <- rnorm(n)
  plot(kde(x), type="l")

  ## mixed normal data
  n <- 1000
  u <- runif(n)
  thresh <- 0.4
  x <- rnorm(n)*(u<thresh) +rnorm(n,mean=3)*(u>=thresh)
  h <- 1
  fh <- kde(x,bandwidth=h)
  plot(kde(x,bandwidth=h),type="l",lwd=2); rug(x)
  lines(kde(x,bandwidth=h*1.2),col="red")
  lines(kde(x,bandwidth=h*1.4),col="orange")
  lines(kde(x,bandwidth=h/1.2),col="blue")
  lines(kde(x,bandwidth=h/1.4),col="cyan")

  ## two-dimensional data
  n <- 1000
  u <- runif(n)
  thresh <- 0.4
  x1 <- rnorm(n)*(u<thresh) +rnorm(n,mean=3)*(u>=thresh)
  x2 <- rnorm(n)*(u<thresh) +rnorm(n,mean=9)*(u>=thresh)

  grid1 <- seq(min(x1),max(x1),length=20)  ## grid for x1
  grid2 <- seq(min(x2),max(x2),length=25)  ## grid for x2

  fh <- kde( cbind(x1,x2), grid=create.grid(list(grid1,grid2)) )
  o <- order(fh$x[,2],fh$x[,1])
  density <- (matrix(fh$y[o],length(grid1),length(grid2)))
  
  par(mfrow=c(2,2))
  plot(kde(x1),type="l",main="x1"); rug(x1)
  plot(kde(x2),type="l",main="x2"); rug(x2)
  persp(grid1,grid2,density,main="KDE",
        theta=30,phi=30,expand=0.5,col="lightblue",shade=0.5)
  contour(grid1,grid2,density, main="KDE Contours")
  points(x1,x2,col="red",pch=18,cex=0.5)
  par(mfrow=c(1,1))

Kernel constants

Description

Calculates several constants of a (product) kernel function.

Usage

kernel.constants(kernel = "biweight", d = 1, product = TRUE)

Arguments

Details

The constants which are calculated are the second moment, the square norm and the canonical bandwidth of the kernel (only the two latter terms depend on the dimension d).

Value

List with components:

Author(s)

Marlene Mueller

See Also

kernel.function

Examples

  kernel.constants()                  ## default (biweight), d=1
  kernel.constants("epanechnikov",1)  ## epanechnikov, d=1
  kernel.constants("epanechnikov",2)  ## product epanechnikov, d=2

Kernel function

Description

Calculates several kernel functions (uniform, triangle, epanechnikov, biweight, triweight, gaussian).

Usage

kernel.function(u, kernel = "biweight", product = TRUE)

Arguments

Details

The kernel parameter is a text string specifying the univariate kernel function which is either the gaussian pdf or proportional to (1-|u|^p)^q. Possible text strings are "triangle" (p=q=1), "uniform" (p=1, q=0), "epanechnikov" (p=2, q=1), "biweight" or "quartic" (p=q=2), "triweight" (p=2, q=3), "gaussian" or "normal" (gaussian pdf).

The multivariate kernels are obtained by a product of unvariate kernels K(u_1)...K(u_d) or by a spherical (radially symmetric) kernel proportional to K(||u||). (The resulting kernel is a density, i.e. integrates to 1.)

Value

n x 1 vector of kernel weights

Author(s)

Marlene Mueller

Examples

  kernel.function(0)                         ## default (biweight)
  kernel.function(0, kernel="epanechnikov")  ## epanechnikov
  kernel.function(0, kernel="gaussian")      ## equals dnorm(0)

Generalized partial linear model

Description

Fits a generalized partial linear model (kernel-based) using the (generalized) Speckman estimator or backfitting (in the generalized case combined with local scoring) for two additive component functions.

Usage

kgplm(x, t, y, h, family, link,
          b.start=NULL, m.start=NULL, grid = NULL, 
          offset = 0, method = "speckman", sort = TRUE, weights = 1,
          weights.trim = 1, weights.conv = 1, max.iter = 25, eps.conv = 1e-8,
          kernel = "biweight", kernel.product = TRUE, verbose = FALSE)

Arguments

Value

List with components:

Author(s)

Marlene Mueller

References

Mueller, M. (2001). Estimation and testing in generalized partial linear models – A comparative study. Statistics and Computing, 11:299–309.

Hastie, T. and Tibshirani, R. (1990). Generalized Additive Models. London: Chapman and Hall.

See Also

kernel.function, kreg

Examples

  ## data
  n <- 1000; b <- c(1,-1); rho <- 0.7
  m <- function(t){ 1.5*sin(pi*t) }
  x1 <- runif(n,min=-1,max=1); u  <- runif(n,min=-1,max=1)
  t  <- runif(n,min=-1,max=1); x2 <- round(m(rho*t + (1-rho)*u))
  x  <- cbind(x1,x2)
  y  <- x %*% b + m(t) + rnorm(n)

  ## partial linear model (PLM)
  gh <- kgplm(x,t,y,h=0.25,family="gaussian",link="identity")
  o <- order(t)
  plot(t[o],m(t[o]),type="l",col="green")
  lines(t[o],gh$m[o]); rug(t)

  ## partial linear probit model (GPLM)
  y <- (y>0)
  gh <- kgplm(x,t,y,h=0.25,family="bernoulli",link="probit")

  o <- order(t)
  plot(t[o],m(t[o]),type="l",col="green")
  lines(t[o],gh$m[o]); rug(t)

  ## data with two-dimensional m-function 
  n <- 1000; b <- c(1,-1); rho <- 0.7
  m <- function(t1,t2){ 1.5*sin(pi*t1)+t2 }
  x1 <- runif(n,min=-1,max=1); u  <- runif(n,min=-1,max=1)
  t1 <- runif(n,min=-1,max=1); t2 <- runif(n,min=-1,max=1)
  x2 <- round( m( rho*t1 + (1-rho)*u , t2 ) )
  x  <- cbind(x1,x2); t  <- cbind(t1,t2)
  y  <- x %*% b + m(t1,t2) + rnorm(n)

  ## partial linear model (PLM)
  grid1 <- seq(min(t[,1]),max(t[,1]),length=20)
  grid2 <- seq(min(t[,2]),max(t[,2]),length=25)
  grid  <- create.grid(list(grid1,grid2))

  gh <- kgplm(x,t,y,h=0.5,grid=grid,family="gaussian",link="identity")

  o <- order(grid[,2],grid[,1])
  est.m  <- (matrix(gh$m.grid[o],length(grid1),length(grid2)))
  orig.m <- outer(grid1,grid2,m)
  par(mfrow=c(1,2))
  persp(grid1,grid2,orig.m,main="Original Function",
        theta=30,phi=30,expand=0.5,col="lightblue",shade=0.5)
  persp(grid1,grid2,est.m,main="Estimated Function",
        theta=30,phi=30,expand=0.5,col="lightblue",shade=0.5)
  par(mfrow=c(1,1))

Kernel regression

Description

Calculates a kernel regression estimate (univariate or multivariate).

Usage

kreg(x, y, bandwidth = NULL, grid = TRUE, kernel = "biweight",
     product = TRUE, sort = TRUE)

Arguments

Details

The estimator is calculated by Nadaraya-Watson kernel regression. Future extension to local linear (d>1) or polynomial (d=1) estimates is planned. The default bandwidth is computed by Scott's rule of thumb for kde (adapted to the chosen kernel function).

Value

List with components:

Author(s)

Marlene Mueller

See Also

kernel.function, convol, kde

Examples

  n <- 1000
  x <- rnorm(n)
  m <- sin(x)
  y <- m + rnorm(n)
  plot(x,y,col="gray")
  o <- order(x); lines(x[o],m[o],col="green")
  lines(kreg(x,y),lwd=2)

  ## two-dimensional
  n <- 100
  x <- 6*cbind(runif(n), runif(n))-3
  m <- function(x1,x2){ 4*sin(x1) + x2 }
  y <- m(x[,1],x[,2]) + rnorm(n)
  mh <- kreg(x,y)##,bandwidth=1)

  grid1 <- unique(mh$x[,1])
  grid2 <- unique(mh$x[,2])
  est.m  <- t(matrix(mh$y,length(grid1),length(grid2)))
  orig.m <- outer(grid1,grid2,m)
  par(mfrow=c(1,2))
  persp(grid1,grid2,orig.m,main="Original Function",
        theta=30,phi=30,expand=0.5,col="lightblue",shade=0.5)
  persp(grid1,grid2,est.m,main="Estimated Function",
	theta=30,phi=30,expand=0.5,col="lightblue",shade=0.5)
  par(mfrow=c(1,1))
  
  ## now with normal x, note the boundary problem,
  ## which can be somewhat reduced by a gaussian kernel
  n <- 1000
  x <- cbind(rnorm(n), rnorm(n))
  m <- function(x1,x2){ 4*sin(x1) + x2 }
  y <- m(x[,1],x[,2]) + rnorm(n)
  mh <- kreg(x,y)##,p="gaussian")

  grid1 <- unique(mh$x[,1])
  grid2 <- unique(mh$x[,2])
  est.m  <- t(matrix(mh$y,length(grid1),length(grid2)))
  orig.m <- outer(grid1,grid2,m)
  par(mfrow=c(1,2))
  persp(grid1,grid2,orig.m,main="Original Function",
        theta=30,phi=30,expand=0.5,col="lightblue",shade=0.5)
  persp(grid1,grid2,est.m,main="Estimated Function",
	theta=30,phi=30,expand=0.5,col="lightblue",shade=0.5)
  par(mfrow=c(1,1))

Generalized partial linear model

Description

Fits a generalized partial linear model (based on smoothing spline) using the (generalized) Speckman estimator or backfitting (in the generalized case combined with local scoring) for two additive component functions. In contrast to kgplm, this function can be used only for a 1-dimensional nonparametric function.

Usage

sgplm1(x, t, y, spar, df=4, family, link,
       b.start=NULL, m.start=NULL, grid = NULL, offset = 0, 
       method = "speckman", weights = 1, weights.trim = 1, 
       weights.conv = 1, max.iter = 25, eps.conv = 1e-8,
       verbose = FALSE, ...)

Arguments

Value

List with components:

Note

This function is mainly implemented for comparison. It is not really optimized for performance, however since it is spline-based, it should be sufficiently fast. Nevertheless, there might be several possibilities to improve for speed, in particular I guess that the sorting that smooth.spline performs in every iteration is slowing down the procedure quite a bit.

Author(s)

Marlene Mueller

References

Mueller, M. (2001) Estimation and testing in generalized partial linear models – A comparative study. Statistics and Computing, 11:299–309.

Hastie, T. and Tibshirani, R. (1990) Generalized Additive Models. London: Chapman and Hall.

See Also

kgplm

Examples

  ## generate data
  n <- 1000; b <- c(1,-1); rho <- 0.7
  mm <- function(t){ 1.5*sin(pi*t) }
  x1 <- runif(n,min=-1,max=1); u  <- runif(n,min=-1,max=1)
  t  <- runif(n,min=-1,max=1); x2 <- round(mm(rho*t + (1-rho)*u))
  x  <- cbind(x1,x2)
  y  <- x %*% b + mm(t) + rnorm(n)

  ## fit partial linear model (PLM)
  k.plm <- kgplm(x,t,y,h=0.35,family="gaussian",link="identity")
  s.plm <- sgplm1(x,t,y,spar=0.95,family="gaussian",link="identity")

  o <- order(t)
  ylim <- range(c(mm(t[o]),k.plm$m,s.plm$m),na.rm=TRUE)
  plot(t[o],mm(t[o]),type="l",ylim=ylim)
  lines(t[o],k.plm$m[o], col="green")
  lines(t[o],s.plm$m[o], col="blue")
  rug(t); title("Kernel PLM vs. Spline PLM")

  ## fit partial linear probit model (GPLM)
  y <- (y>0)
  k.gplm <- kgplm(x,t,y,h=0.35,family="bernoulli",link="probit")
  s.gplm <- sgplm1(x,t,y,spar=0.95,family="bernoulli",link="probit")

  o <- order(t)
  ylim <- range(c(mm(t[o]),k.gplm$m,s.gplm$m),na.rm=TRUE)
  plot(t[o],mm(t[o]),type="l",ylim=ylim)
  lines(t[o],k.gplm$m[o], col="green")
  lines(t[o],s.gplm$m[o], col="blue")
  rug(t); title("Kernel GPLM vs. Spline GPLM (Probit)")