Nothing
# file NominalLogisticBiplot/R/polylogist.R
# copyright (C) 2012-2013 J.L. Vicente-Villardon and J.C. Hernandez
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2 or 3 of the License
# (at your option).
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# A copy of the GNU General Public License is available at
# http://www.r-project.org/Licenses/
#
#Function that calculates an object with the fitted multinomial logistic regression for a nominal variable.
#----------------------Parameters--------------
#y: response nominal variable.
#x: matrix with independent variables.
#penalization: value to correct the separation problem through the ridge regression
#cte : it will be true if the model has an independent term.
#tol : value to decide if the algorith should continue
#maxiter : value to decide if the algorith should continue
#show: boolean parameter if we want to see values in each iteration of the process.
polylogist <- function(y, x, penalization = 0.2, cte = TRUE, tol = 1e-04, maxiter = 200,show=FALSE) {
# POLYLOGIST Fits a polytomous ridge logistic regression
# Produces a vector B of estimates of parameters
# in a polytomoius logistic regression
# The response vector must be constructed with integers starting at 0
# the category labeled 0 is used for comparisons
# The rows of B are for the categories and the column for the variables
# To break the singularity of the information matrix a constant (penalization) is added
# to the diagonal
if(is.matrix(x)){
n <- nrow(x)
}else{
n <- length(x)
}
if (min(y)==1) y= y-1
if (cte){
x <- cbind(matrix(1,n),x)
}
p <- ncol(x)
J = max(y)
Y = matrix(0, n, J)
for (i in 1:n){
if (y[i] > 0){
Y[i, y[i]] = 1
}
}
beta <- matrix(0,J,p)
yy = as.vector(Y)
err = 0.1
iter = 0
m = matrix(0, n * J, n * J)
prob = matrix(0, n, J)
p0 = matrix(0, n, 1)
xx = kronecker(diag(1, J), x)
while ((err > tol) & (iter < maxiter)) {
iter = iter + 1
betaold = beta
b = as.vector(t(beta))
for (i in 1:n) {
suma = 1
for (j in 1:J) suma = suma + exp(sum(beta[j, ] * x[i, ]))
for (j in 1:J) prob[i, j] = exp(sum(beta[j, ] * x[i, ]))/suma
p0[i, 1] = 1/suma
}
P = as.vector(prob)
for (j in 1:J)
for (k in 1:J){
mj = matrix(0, n, n)
for (i in 1:n) {
if (k == j) {
mj[i, i] = prob[i, j] * (1 - prob[i, j])
}
if (!(k == j)) {
mj[i, i] = -1 * prob[i, j] * prob[i, k]
}
}
m[((j - 1) * n + 1):(j * n), ((k - 1) * n + 1):(k * n)] = mj
}
a = (t(xx) %*% m %*% xx) + 2 * diag(J*p) * penalization
u = (t(xx) %*% matrix((yy - P), n * J, 1)) - 2 * b * penalization
b = b + solve(a) %*% u
beta = t(matrix(b, p, J))
err = sum(sum((betaold - beta) * (betaold - beta)))
#if (show) print(c(iter, err))
}
y0=matrix(0,n,1)
for (i in 1:n) if (sum(Y[i,]) == 0) y0[i]=1
cov=solve(a)
model<-list()
model$fitted=cbind(p0,prob)
model$cov=cov
model$Y=cbind(y0,Y)
model$beta=beta
model$stderr=sqrt(t(matrix(diag(cov), p, J)))
model$logLik=sum(log(model$fitted^model$Y))
model$Deviance = -2*sum(log(model$fitted^model$Y))
model$AIC = model$Deviance + 2*nrow(beta)*ncol(beta)
model$BIC = model$Deviance + nrow(beta)*ncol(beta)*log(n)
class(model)='polylogist'
return(model)
}
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.