Nothing
#' Fit a PIM
#'
#' Fit a PIM
#'
#' @aliases pim pim-class glmnetpim glmnetpim-class
#'
#' @param formula Original formula
#' @param data Context where the formula \code{formula} is to be interpreted
#' @param link Name of the link function (defaults to \code{"logit"})
#' @param blocking.variables Character vector holding column names that hold blocking
#' variables.
#' @param poset Matrix of two columns indicating what the original observation number is
#' for the left and right real observation in the pseudo-observation. Alternatively, can
#' be a function like \code{\link{fullposet}} and similar.
#' @param leftsuffix,rightsuffix Suffixes that will be added to the 'left' and 'right'
#' observation's column name in the pseudo-observation. Note: no checking is done that
#' these suffixes are safe, so the wrong suffixes may lead to unexpected behaviour.
#' @param interpretation If \code{"marginal"} parts of the formula are converted to imply
#' marginal pim modeling (see e.g. \code{\link{Mainreplacetext}}). If it is \code{"difference"}
#' (the default), then the design matrix of the PIM is the difference of the
#' design matrices of each part of the pseudo-observations. A popular option is
#' \code{"regular"}, which will interpret unaltered columns as differences. A new option
#' is \code{"symmetric"}, which works the same as \code{"regular"}, but will enforce
#' the symmetry condition by making the sign switch when changing the order (typically,
#' this is achieved by subtracting the inverse for each dummy).
#' @param na.action Defaults to \code{\link{na.fail}}: handles missing data in \code{data}.
#' @param estimator Function like the result of \code{\link{estimator.nleqslv}()} (the
#' default). See there to find the required form of this function or some provided
#' alternatives.
#' @param varianceestimator Function like the result of \code{\link{varianceestimator.sandwich}()}
#' (the default). See there to find the required form of this function or some
#' provided alternatives.
#' @param lhs \code{"PO"}, \code{"<"} or \code{"<="}: Unequality used for the lefthandside of
#' the formula. The default (\code{"PO"}) is the normal probabilistic index.
#' @param keep.data If \code{TRUE} (not the default), the data and design matrices are maintained
#' in the return value
#' @param verbosity The higher this value, the more levels of progress and debug
#' information is displayed (note: in R for Windows, turn off buffered output)
#' @param nicenames Defaults to \code{TRUE}: try to make the column names more readable.
#' @param interactions.difference If \code{TRUE} (note that the default is
#' \code{interpretation!="marginal"}) interaction terms will be interpreted as the
#' differences of the onesided interaction terms (if this is possible at all). This is
#' unsupported if \code{unsupported if} is \code{"marginal"}. Some special interaction
#' terms with calculated columns may lead to unexpected behaviour.
#' @param extra.nicenames Should be a \code{data.frame} containing two character columns:
#' \code{org} and \code{nice}. For each "constructed" column name, provide a nicer one,
#' that will make the results more readable. You may also use parts of constructed column
#' names. Note: make sure to use \code{stringsAsFactor=FALSE} when creating the \code{data.frame}.
#' @param check.symmetric Defaults to \code{TRUE} for \code{interpretation=="regular"},
#' \code{FALSE} otherwise: if the model does not support the
#' symmetry condition, a warning is displayed. Note: this option is on by default for
#' novice users, but it may be time- and resourceconsuming, so turn it off if you know
#' what you are doing!
#' @param threshold When checking the symmetry condition, how much digression is allowed.
#' @param weights Defaults to \code{NULL}: vector of weights for every row of \code{data}, or
#' \code{NULL} if no weighting is required.
#' @param pseudoweights Defaults to \code{\link{pseudoweights.default}}: function that can convert
#' weights by observation to weights per pseudo-observation. Should have the same signature
#' and outcome as \code{\link{pseudoweights.default}}.
#' @return An object of class \code{pim},which is a list with items:
#' \item{coefficients}{Parameter estimates. See \code{\link{estimator.nleqslv}}}
#' \item{morefitinfo }{Parameter estimation info. See \code{\link{estimator.nleqslv}}}
#' \item{vcov }{Covariance estimates. See \code{\link{varianceestimator.sandwich}}}
#' \item{morevarfitinfo }{Covariance estimation info. See \code{\link{varianceestimator.sandwich}}}
#' \item{varestimationerror }{Only present if covariance estimation failed. See \code{\link{pim.fit}}}
#' \item{call }{Call used when running this function.}
#' \item{formula }{Passed in formula.}
#' \item{link }{Name of the link function (defaults to \code{"identity"})}
#' \item{fitted.values }{Predicted pseudo-observation values.}
#' \item{pfd }{Object of class \code{\link{pimfitdata}}. If \code{keep.data} was \code{TRUE}, the
#' \code{X} and \code{Y} items are set to \code{NULL}.}
#' \item{data }{Passed in \code{data} (may be altered by \code{poset}). Only present if
#' \code{keep.data} was \code{TRUE}}
#' @details Some of the items in the return value may be missing if parts of the estimation
#' fail. See \code{\link{pim.fit}} for this.
#'
#' If \code{estimator} was \code{estimator.glmnet}, The return value gets as an additional class
#' \code{"glmnetpim"}.
#'
#' For more details on how to use this function and its parameters, we refer the reader to
#' \code{vignette("pim")}.
#' @seealso \code{\link{estimator.nleqslv}} \code{\link{varianceestimator.sandwich}}
#' @keywords pim fit
#' @examples set.seed(1)
#' iris$out<-factor(sample(2, nrow(iris), replace=TRUE))
#' iris$xord<-as.ordered(iris$Species)
#' pima<-pim(out~Sepal.Length, data=iris, link="logit", interpretation="regular")
#' pimb<-pim(out~I((R(Sepal.Length) - L(Sepal.Length))/sqrt(R(Sepal.Length) * L(Sepal.Length)) ), data=iris,
#' link="logit", interpretation="regular", extra.nicenames=data.frame(
#' org="I((R(Sepal.Length) - L(Sepal.Length))/sqrt(R(Sepal.Length) * L(Sepal.Length)) )",
#' nice="Sepal.Length.WDiff", stringsAsFactors=FALSE))
#' pimc<-pim(out~O(xord), data=iris, link="logit", interpretation="regular")
#' ##And an example with weighted fitting
#' iris50<-iris[sample.int(nrow(iris), 50),]
#' rwt50<-runif(50) #50 random weights
#' pimd<-pim(out~O(xord), data=iris50, link="logit", interpretation="regular",
#' weights = rwt50)
#' @export
pim<-function(formula, data, link=c("logit", "identity", "probit", "inverse", "1/mu^2", "log"),
blocking.variables=character(),
poset=noselfposet, leftsuffix="_L", rightsuffix="_R",
interpretation=c("difference", "regular", "marginal", "symmetric"), na.action=na.fail,
estimator=estimator.nleqslv(), varianceestimator=varianceestimator.sandwich(),
lhs=c("PO", "<", "<="), keep.data=FALSE, verbosity=0,
nicenames=TRUE, interactions.difference=(interpretation!="marginal"),
extra.nicenames=data.frame(org=character(), nice=character(), stringsAsFactors=FALSE),
check.symmetric=(interpretation == "regular"), threshold=1e-6, weights=NULL,
pseudoweights=pseudoweights.default)
{
this.call<-match.call()
link<-match.arg(link)
.fulfillreqs(estimator)
.fulfillreqs(varianceestimator)
if(is.function(poset))
{
poset<-poset(data, formula, weights=weights, verbosity=verbosity-1)
data<-poset$data
weights<-poset$weights
poset<-poset$poset
}
if(verbosity>0) cat("Get design matrix and converted formula\n")
pfd<-pim.fit.prep(formula=formula, data=data, blocking.variables=blocking.variables, poset=poset, leftsuffix=leftsuffix, rightsuffix=rightsuffix,
interpretation=interpretation, na.action=na.action, lhs=lhs, verbosity=verbosity-1, nicenames=nicenames,
interactions.difference=interactions.difference, extra.nicenames=extra.nicenames, check.symmetric=check.symmetric, link=link,
weights=weights, pseudoweights=pseudoweights)
if(verbosity>0) cat("Coefficient and covariance estimates\n")
actualfit<-pim.fit(pfd=pfd, link=link, estimator=estimator, varianceestimator=varianceestimator, verbosity=verbosity-1, nicenames=nicenames)
actualfit$call<-this.call
actualfit$formula<-formula
actualfit$link<-link
actualfit$pfd<-pfd
if(!keep.data)
{
actualfit$pfd$X<-NULL
}
else
{
actualfit$data<-data
}
actualfit$blocking.variables<-blocking.variables
class(actualfit)<-"pim"
#this may not be the cleanest way to do this. Reconsider in the future?
if("glmnet" %in% class(actualfit$morefitinfo)) class(actualfit)<-c("glmnetpim", class(actualfit))
return(actualfit)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.