harsm | R Documentation |
A user friendly interface to the softmax regression under the Harville model.
harsm(
formula,
data,
group = NULL,
weights = NULL,
fit0 = NULL,
reg_wt = NULL,
reg_zero = 0,
reg_power = NULL,
reg_coef_idx = NULL,
reg_standardize = FALSE,
na.action = na.omit
)
## S3 method for class 'harsm'
vcov(object, ...)
## S3 method for class 'harsm'
print(x, ...)
formula |
an object of class |
data |
an optional data frame, list or environment (or object
coercible by |
group |
the string name of the group variable in the data, or a bare character with the group name. The group indices need not be integers, but that is more efficient. They need not be sorted. |
weights |
an optional vector of weights, or the string or bare name of the
weights in the |
fit0 |
An optional object of class |
reg_wt |
the multiplicative weight(s) of the regularization terms.
must be non-negative.
May be a scalar or vector, but will be recycled to the length of the |
reg_zero |
the ‘zero’ of the regularization terms. This would
usually be zero if you want to shrink to zero, but in some cases you may
with to shrink to 1 for example.
May be a scalar or vector, but will be recycled to the length of the |
reg_power |
the power of the regularization terms, 2 for ridge regression, 1 for lasso. |
reg_coef_idx |
the index of the coefficient which the corresponding
regularization term is applied to. For the Harville model, the indices only refer to
the |
reg_standardize |
if true, the |
na.action |
How to deal with missing values in the outcomes, groups, weights, etc. |
object |
an object of class |
... |
For |
x |
an object used to select a method. |
Performs a softmax regression by groups, via Maximum Likelihood Estimation,
under the Harville model.
We fit \beta
where odds are \eta = x^{\top}\beta
for
independent variables x
.
The probability of taking first place is then \mu=c\exp{\eta}
,
where the c
is chosen so the \mu
sum to one.
Under the Harville model, conditional on the first place finisher
having been observed, the probability model for second
(and successive) places with the probabilities of the remaining
participants renormalized.
The print
method of the harsm
object includes
a display of the R-squared. This measures the improvement
in squared errors of the expected rank from the model
over the null model which posits that all odds are equal.
When the formula includes an offset, a ‘delta R-squared’
is also output. This is the improvement in predicted
ranks over the model based on the offset term.
Note that the expected ranks are only easy to produce
under the Harville model; under the Henery model,
the summary R-squared statistics are not produced.
Note that this computation uses weighted sums of squares,
as the weights contribute to the likelihood term.
However, the square sum computation does not take into
account the standard error of the rank, and so
unlike in linear regression, the softmax regression
does not always give positive R-squareds,
and the statistic is otherwise hard to interpret.
An object of class harsm
, but also of maxLik
with the
fit.
The regularization term is of the form
\sum_i w_i |\nu_{c_i} - z_i|^{p_i},
where w_i
are the reg_wt
weights,
z_i
are the reg_zero
zeroes,
p_i
are the reg_power
powers,
and c_i
are the reg_coef_idx
coefficient indices.
Note that the coefficient indices can be repeated so that the
regularization term can include multiple contributions from one
coefficient element. This allows ‘elasticnet’ regularization.
The \nu
here refer to the regression coefficients \beta
concatenated with the \gamma
coefficients in the Henery model.
The fit functions return an object of type maxLik
even when regularization penalties are applied; the statistical inference
functions are not valid when regularization is used. The user is warned.
This regression may give odd results when the outcomes are tied, imposing an arbitrary order on the tied outcomes. Moreover, no warning may be issued in this case. In future releases, ties may be dealt with differently, perhaps in analogy to how ties are treated in the Cox Proportional Hazards regression, using the methods of Breslow or Efron.
To avoid incorrect inference when only the top performers are recorded, and all others are effectively tied, one should use weighting. Set the weights to zero for participants who are tied non-winners, and one for the rest So for example, if you observe the Gold, Silver, and Bronze medal winners of an Olympic event that had a starting field of 12 participants, set weights to 1 for the medal winners, and 0 for the others. Note that the weights do not attach to the participants, they attach to the place they took.
Since version 0.1.0 of this package, the normalization of weights used in this function have changed under the hood. This is to give correct inference in the case where zero weights are used to signify finishing places were not observed. If in doubt, please confirm inference by simulations, taking as example the simulations in the README.
Steven E. Pav shabbychef@gmail.com
harsmfit
, harsmlik
.
nfeat <- 5
set.seed(1234)
g <- ceiling(seq(0.1,1000,by=0.1))
X <- matrix(rnorm(length(g) * nfeat),ncol=nfeat)
beta <- rnorm(nfeat)
eta <- X %*% beta
y <- rsm(eta,g)
# now the pretty frontend
data <- cbind(data.frame(outcome=y,race=g),as.data.frame(X))
fmla <- outcome ~ V1 + V2 + V3 + V4 + V5
fitm <- harsm(fmla,data,group=race)
eta0 <- rowMeans(X)
data <- cbind(data.frame(outcome=y,race=g,eta0=eta0),as.data.frame(X))
fmla <- outcome ~ offset(eta0) + V1 + V2 + V3 + V4 + V5
fitm <- harsm(fmla,data,group=race)
# with weights
data <- cbind(data.frame(outcome=y,race=g,eta0=eta0),as.data.frame(X))
data$wts <- runif(nrow(data),min=1,max=2)
fmla <- outcome ~ offset(eta0) + V1 + V2 + V3 + V4 + V5
fitm <- harsm(fmla,data,group=race,weights=wts)
# softmax on the Best Picture data
data(best_picture)
df <- best_picture
df$place <- ifelse(df$winner,1,2)
df$weight <- ifelse(df$winner,1,0)
fmla <- place ~ nominated_for_BestDirector + nominated_for_BestActor + Drama
fit0 <- harsm(fmla,data=df,group=year,weights=weight)
# warm start is a thing:
sub_fmla <- place ~ nominated_for_BestDirector + nominated_for_BestActor
fit1 <- harsm(sub_fmla,data=df,group=year,weights=weight,fit0=fit0)
# ridge regression.
fitr2 <- harsm(sub_fmla,data=df,group=year,weights=weight,
reg_wt=rep(1,2), reg_power=2, reg_zero=0, reg_coef_idx=c(1,2))
# l1 regularization regression.
fitr1 <- harsm(sub_fmla,data=df,group=year,weights=weight,
reg_wt=rep(1,2), reg_power=1, reg_zero=0, reg_coef_idx=c(1,2))
# elasticnet regularization regression.
fitr12 <- harsm(sub_fmla,data=df,group=year,weights=weight,
reg_wt=rep(1,4), reg_power=c(1,1,2,2), reg_zero=0, reg_coef_idx=c(1,2,1,2))
# test against logistic regression
if (require(dplyr)) {
nevent <- 10000
set.seed(1234)
adf <- data_frame(eventnum=floor(seq(1,nevent + 0.7,by=0.5))) %>%
mutate(x=rnorm(n()),
program_num=rep(c(1,2),nevent),
intercept=as.numeric(program_num==1),
eta=1.5 * x + 0.3 * intercept,
place=ohenery::rsm(eta,g=eventnum))
# Harville model
modh <- harsm(place ~ intercept + x,data=adf,group=eventnum)
# the collapsed data.frame for glm
ddf <- adf %>%
arrange(eventnum,program_num) %>%
group_by(eventnum) %>%
summarize(resu=as.numeric(first(place)==1),
delx=first(x) - last(x),
deli=first(intercept) - last(intercept)) %>%
ungroup()
# glm logistic fit
modg <- glm(resu ~ delx + 1,data=ddf,family=binomial(link='logit'))
all.equal(as.numeric(coef(modh)),as.numeric(coef(modg)),tolerance=1e-4)
all.equal(as.numeric(vcov(modh)),as.numeric(vcov(modg)),tolerance=1e-4)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.