# R/circfit.R In circmax: Circular Regression with Maximum Likelihood Estimation and Regression Trees

```circfit <- function(y, x = NULL, start = NULL, subset = NULL, na.action = NULL,
weights = NULL, offset = NULL, ...,
vcov = TRUE, estfun = TRUE, object = FALSE, fit_control = circfit_control()) {

# TODO: Why is circfit quite often called with object TRUE and what are the start values used for?!
# if(object) cat("now object = TRUE\n")

# Check unsupported arguments
if(!(is.null(x) || NCOL(x) == 0L)) warning("no regression coefficients
are currently taken into account..")
if(!is.null(subset)) warning("'subset' currently not supported and therefore not used")
if(!is.null(na.action)) warning("'na.action' currently not supported and therefore not used")
if(!is.null(offset)) warning("'offset' currently not supported and therefore not used")

## Convenience variables
n <- NROW(y)
allequy <- (length(unique(y)) == 1)
cl <- match.call()

## Weights
if(is.null(weights) || (length(weights)==0L)) weights <- as.vector(rep.int(1, n))
if(unique(weights) == 0L) stop("weights are not allowed to be all zero")
if(length(weights) != n) stop("number of observations and length of weights are not equal")

## Control parameters
solve_kappa <- fit_control\$solve_kappa
useC <- fit_control\$useC
ncores <- fit_control\$ncores
fit_control\$solve_kappa <- fit_control\$useC <- fit_control\$ncores <- NULL

## Get Von Mises family functions (for interchangeability the same distlist is used)
circfam <- dist_vonmises(useC = useC, ncores = ncores)

## MLE according to Bettina Gruen
eta <- circfam\$startfun(y, weights = weights, solve_kappa = solve_kappa)

## Compute negative loglik
nll <- -circfam\$ddist(y, eta, log = TRUE, weights = weights, sum = TRUE)

# FIXME: Check if estfun and vcov are correct (at least they are identical to Lisa's).
if(estfun) {
if(allequy) {
ef <- matrix(0, ncol = length(eta), nrow = n)
} else {
ef <- as.matrix(weights * circfam\$sdist(y, eta, sum = FALSE))
}
} else {
ef <- NULL
}

if(vcov) {
hess <- circfam\$hdist(y, eta, weights = weights)

## Convert to matrix with nice column names
hess <- as.matrix(hess)
colnames(hess) <- rownames(hess) <- names(eta)

vc <- try(solve(-hess), silent = TRUE)
if(inherits(vc, "try-error")) {
vc <- try(qr.solve(-hess), silent = TRUE)
if(inherits(vc, "try-error")) {
vc <- try(chol2inv(chol(-hess)))
if(inherits(vc, "try-error")) {
warning("hessian matrix is not invertible ('numerically' singular)")
}
}
}

## Convert to matrix with nice column names
vc <- as.matrix(vc)
colnames(vc) <- rownames(vc) <- colnames(hess)

} else {
hess <- NULL
vc <- NULL
}

# FIXME: Which model should I return, is there a mistake in the mobster?!
if(vcov) object <- TRUE
model <- list(
npar = length(par),
y = y,
ny = n,
weights = weights,
family = circfam,
start = start,
par = par,
eta = eta,
hess = hess,
vcov = vc,
loglik = -nll,
call = cl,
estfun = ef
)
class(model) <- "circfit"

list(coefficients = par,
objfun = nll,
estfun = ef,
object = if(object) model else NULL)
}

## Different methods for circfit
# FIXME: All models go for the fit\$object: good, bad?!

nobs.circfit <- function(object, ...) {
return(object\$ny)
}

coef.circfit <- function(object, type = "parameter" , ...) {
if(type == "parameter") return(object\$par)
}

predict.circfit <- function(object, type = "parameter", ...){
if(type == "parameter") {
return(object\$par)
} else {
# FIXME: Implement type = 'response'
stop("currently just type = 'parameter' implemented")
}
}

vcov.circfit <- function(object, type = "link", ...) {
if(type == "parameter"){
## delta method
colnames(delta.m) <- rownames(delta.m) <- names(object\$par)
return(delta.m %*% object\$vcov %*% delta.m)
}
}

estfun.circfit <- function(object, ...) {
return(object\$estfun)
}

logLik.circfit <- function(object, ...) {
structure(object\$loglik, df = object\$npar, class = "logLik")
}

type <- match.arg(type)
if(type == "parameter") return(vcov(object, type = "parameter") * object\$ny)
if(type == "link") return(object\$vcov * object\$ny)
}

print.circfit <- function(x, digits = max(3, getOption("digits") - 3), ...) {
cat("Fitted ", x\$family\$family.name, "\n\n")
if(length(x\$par)) {
cat("Distribution parameter(s):\n")
print.default(format(x\$par, digits = digits), print.gap = 2, quote = FALSE)
cat("\n")
} else {
cat("No parameters \n\n")
}
cat(paste("Log-likelihood: ", format(x\$loglik, digits = digits), "\n", sep = ""))
if(length(x\$npar)) {
cat(paste("Df: ", format(x\$npar, digits = digits), "\n", sep = ""))
}
cat("\n")

invisible(x)
}
```

## Try the circmax package in your browser

Any scripts or data that you put into this service are public.

circmax documentation built on May 2, 2019, 5:16 p.m.