Nothing
# file NominalLogisticBiplot/R/Generators.R
# copyright (C) 2012-2013 J.C. Hernandez and J.L. Vicente-Villardon
#
# 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 information about the tesselation, centers, real an dummy points
#, hide categories, etc... from the variable of study.
#----------------------Parameters--------------
#beta: matrix of parameters for the multinomial logistic regression calculated for one specific dependent variable.
Generators <- function(beta) {
ngrup = nrow(beta) + 1
a = matrix(0, ngrup, ngrup)
b = matrix(0, ngrup, ngrup)
# Calculate the straight lines separating each pair of categories
for (i in 1:(ngrup - 1)) {
a[1, (i + 1)] = -1 * beta[i, 1]/beta[i, 3]
b[1, (i + 1)] = -1 * beta[i, 2]/beta[i, 3]
}
if (ngrup > 2) {
for (i in 1:(ngrup - 2)) for (j in (i + 1):(ngrup - 1)) {
a[(i + 1), (j + 1)] = (beta[j, 1] - beta[i, 1])/(beta[i, 3] - beta[j, 3])
b[(i + 1), (j + 1)] = (beta[j, 2] - beta[i, 2])/(beta[i, 3] - beta[j, 3])
}
}
# Calculate the coordinates of the points in which three lines intersect
# Those are the candidates to form the voronoi diagram
triplepoints = (ngrup - 2) * (ngrup - 1) * ngrup/6
xp = matrix(0, triplepoints, 1)
yp = matrix(0, triplepoints, 1)
indices = matrix(0, triplepoints, 3)
if (ngrup > 3) {
l = 0
for (i in 1:(ngrup - 2)) for (j in (i + 1):(ngrup - 1)) for (k in (j + 1):(ngrup)) {
l = l + 1
xp[l] = (a[i, k] - a[i, j])/(b[i, j] - b[i, k])
yp[l] = a[i, j] + b[i, j] * xp[l]
indices[l, ] = c(i, j, k)
}
}else if(ngrup == 3){
xp[1] = (a[1, 3] - a[1, 2])/(b[1, 2] - b[1, 3])
yp[1] = a[1, 2] + b[1, 2] * xp[1]
indices[1, ] = c(1,2,3)
}
# Decide if the candidates are real or virtual using the probabilities given by the model
# Probabilities for each category at the candidate points.
# Real points will be in the diagram, virtual points will be used to determine the dummy points
xn = cbind(matrix(1, triplepoints, 1), xp, yp)
probab = matrix(0, triplepoints, ngrup)
for (i in 1:triplepoints) {
suma = 1
for (j in 1:(ngrup - 1)) suma = suma + exp(sum(beta[j, ] * xn[i, ]))
for (j in 1:(ngrup - 1)) probab[i, (j + 1)] = exp(sum(beta[j, ] * xn[i, ]))/suma
probab[i, 1] = 1/suma
}
# When a point is too far away the probability of one category is 1 and the rest are 0, but sometimes
# the calculation results in a NaN
NoNumbers =which(is.nan(probab))
for (i in NoNumbers) probab[i]=1
# Calculating the status (real or virtual) for each point. A point is real when one its three
# common categories has the highest probability, and virtual otherwise
# the procedure is good when there are more than 3 categories
status = matrix(1, triplepoints, 1) #1 for real, 0 for virtual
if (ngrup > 3) {
l = 0
nvirtual = 0
for (i in 1:(ngrup - 2))
for (j in (i + 1):(ngrup - 1))
for (k in (j + 1):(ngrup)){
l = l + 1
if (max(probab[l, ]) > max(c(probab[l, i], probab[l, j], probab[l, k]))) {
status[l] = 0
nvirtual = nvirtual + 1
}
}
# Separating the real and virtual points, and the corresponding indices
nreal = triplepoints - nvirtual
coorreal = matrix(0, nreal, 2)
coorvirt = matrix(0, nvirtual, 2)
IndReal = matrix(1, nreal, 3)
IndVirt = matrix(1, nvirtual, 3)
nr = 0
nv = 0
for (i in 1:triplepoints) {
if (status[i] == 1) {
nr = nr + 1
coorreal[nr, ] = matrix(c(xp[i], yp[i]), 2, 1)
IndReal[nr, ] = indices[i, ]
}
if (status[i] == 0) {
nv = nv + 1
coorvirt[nv, ] = c(xp[i], yp[i])
IndVirt[nv, ] = indices[i, ]
}
}
}else if(ngrup == 3){
nvirtual = 0
nreal = 1
coorreal = matrix(0, nreal, 2)
IndReal = matrix(1, nreal, 3)
nr = 1
nv = 0
coorreal[nr, ] = matrix(c(xp[1], yp[1]), 2, 1)
IndReal[nr, ] = indices[1, ]
}
# Calculating the tesellation #############
# Each real point should be connected to other three (either real or virtual)
# at least another one has to be real.
# the virtual points are connected just to one real point.
Joint = matrix(0, nr, nr)
# This is the incidence matrix, square matrix with the real points in rows and columns
# the elements are 1 when te points are connected and 0 otherwise
if((ngrup > 3) & (nr > 1)){
for (i in 1:(nr - 1))
for (j in (i + 1):nr)
if (length(intersect(IndReal[i, ], IndReal[j, ])) == 2) {
Joint[i, j] = 1
Joint[j, i] = 1
}
}
# variables for the coordinates of the dummy points
dummy.x = NULL
dummy.y = NULL
JointDummy = NULL
IndDummy = NULL
# n1, n2 & n3 contain the nodes to which each real point is connected
# when the node is dummy, the number is negative
n1 = rep(0, nr)
n2 = rep(0, nr)
n3 = rep(0, nr)
ndummy = 0
for (i in 1:nr) {
if (sum(Joint[i, ]) == 3) {
# When a real point is connected to other three real points, search which are the nodes and store them
neighbors = which(Joint[i, ] == 1)
n1[i] = neighbors[1]
n2[i] = neighbors[2]
n3[i] = neighbors[3]
}
if (sum(Joint[i, ]) == 2) {
# When a real point is connected to other 2 real points, search which are the nodes and store them
# then we need a dummy point. The dummy point is obtained using the virtual points.
# The virtual point is located on the remaining line (at each real point intersect exactly 3 lines and we have already used 2)
# on the remaining line there should be several virtual points (depending on the number of categories), one of them
# is enough to calculate the dummy point. Calculating the parameters of the straight line passing through the
# real and the virtual points, the dumme point is located on the opposite direction to the virtual
# point relative to the real one. The dummy point is placed symmetrically, i. e., the distance
# from the real to the virtual and the distance to the dummy are the same
ndummy = ndummy + 1
ColJointDummy = matrix(0, nr, 1)
ColJointDummy[i] = 1
JointDummy = cbind(JointDummy, ColJointDummy)
neighbors = which(Joint[i, ] == 1)
n1[i] = -1 * ndummy
n2[i] = neighbors[1]
n3[i] = neighbors[2]
A = intersect(IndReal[i, ], IndReal[neighbors[1], ])
B = intersect(IndReal[i, ], IndReal[neighbors[2], ])
C = setdiff(union(A, B), intersect(A, B))
IndDummy = rbind(IndDummy, sort(C))
for (j in 1:nv) if (length(intersect(C, IndVirt[j, ])) == 2) {
pend = (coorvirt[j, 2] - coorreal[i, 2])/(coorvirt[j, 1] - coorreal[i, 1])
const = coorreal[i, 2] - pend * coorreal[i, 1]
cx = coorreal[i, 1] - 3 * (coorvirt[j, 1] - coorreal[i, 1])
cy = const + pend * cx
dummy.x = c(dummy.x, cx)
dummy.y = c(dummy.y, cy)
break
}
}
if (sum(Joint[i, ]) == 1) {
# When a real point is connected to just another real point, we need two virtual points on the non used lines
# The procedure to calculate the dummy points is as before
n3[i] = which(Joint[i, ] == 1)
neighbors = which(Joint[i, ] == 1)
A = intersect(IndReal[i, ], IndReal[neighbors[1], ])
C = setdiff(IndReal[i, ], A)
D = c(C, A[1])
E = c(C, A[2])
IndDummy = rbind(IndDummy, sort(D))
IndDummy = rbind(IndDummy, sort(E))
for (j in 1:nv) if (length(intersect(D, IndVirt[j, ])) == 2) {
ndummy = ndummy + 1
ColJointDummy = matrix(0, nr, 1)
ColJointDummy[i] = 1
JointDummy = cbind(JointDummy, ColJointDummy)
n1[i] = -1 * ndummy
pend = (coorvirt[j, 2] - coorreal[i, 2])/(coorvirt[j, 1] - coorreal[i, 1])
const = coorreal[i, 2] - pend * coorreal[i, 1]
cx = coorreal[i, 1] - 3 * (coorvirt[j, 1] - coorreal[i, 1])
cy = const + pend * cx
dummy.x = c(dummy.x, cx)
dummy.y = c(dummy.y, cy)
break
}
for (j in 1:nv) if (length(intersect(E, IndVirt[j, ])) == 2) {
ndummy = ndummy + 1
ColJointDummy = matrix(0, nr, 1)
ColJointDummy[i] = 1
JointDummy = cbind(JointDummy, ColJointDummy)
n2[i] = -1 * ndummy
pend = (coorvirt[j, 2] - coorreal[i, 2])/(coorvirt[j, 1] - coorreal[i, 1])
const = coorreal[i, 2] - pend * coorreal[i, 1]
cx = coorreal[i, 1] - 3 * (coorvirt[j, 1] - coorreal[i, 1])
cy = const + pend * cx
dummy.x = c(dummy.x, cx)
dummy.y = c(dummy.y, cy)
break
}
}
if (sum(Joint[i, ]) == 0) {
IndDummy = matrix(0,3,2)
IndDummy[1,] = c(IndReal[1,1],IndReal[1,2])
IndDummy[2,] = c(IndReal[1,1],IndReal[1,3])
IndDummy[3,] = c(IndReal[1,2],IndReal[1,3])
for(i in 1:3){
if(i == 1){n1[1] = -1}
if(i == 2){n2[1] = -2}
if(i == 3){n3[1] = -3}
xnTR = cbind(matrix(1, 1, 1), xn[1,2] + 2, a[IndDummy[i,1],IndDummy[i,2]] + b[IndDummy[i,1],IndDummy[i,2]]*(xn[1,2] + 2))
xnTL = cbind(matrix(1, 1, 1), xn[1,2] - 2, a[IndDummy[i,1],IndDummy[i,2]] + b[IndDummy[i,1],IndDummy[i,2]]*(xn[1,2] - 2))
pisubiR = matrix(0, 1, ngrup)
denom = 1
for (j in 1:(ngrup - 1)) denom = denom + exp(sum(beta[j, ] * xnTR[1, ]))
for (j in 1:(ngrup - 1)) pisubiR[1, (j + 1)] = exp(sum(beta[j, ] * xnTR[1, ]))/denom
pisubiR[1, 1] = 1/denom
pisubiL = matrix(0, 1, ngrup)
denom = 1
for (j in 1:(ngrup - 1)){
denom = denom + exp(sum(beta[j, ] * xnTL[1, ]))
}
for (j in 1:(ngrup - 1)) pisubiL[1, (j + 1)] = exp(sum(beta[j, ] * xnTL[1, ]))/denom
pisubiL[1, 1] = 1/denom
if ((pisubiL[1,IndDummy[i,1]] > pisubiR[1,IndDummy[i,1]]) &
(pisubiL[1,IndDummy[i,2]] > pisubiR[1,IndDummy[i,2]])) {
cx = xn[1,2] - 2
cy = a[IndDummy[i,1],IndDummy[i,2]] + b[IndDummy[i,1],IndDummy[i,2]]*(xn[1,2] - 2)
}else{
cx = xn[1,2] + 2
cy = a[IndDummy[i,1],IndDummy[i,2]] + b[IndDummy[i,1],IndDummy[i,2]]*(xn[1,2] + 2)
}
dummy.x = c(dummy.x, cx)
dummy.y = c(dummy.y, cy)
}
ndummy = 3
}
}
if((ngrup > 3)&(nr > 1)){
Borders = matrix(0, sum(sum(Joint))/2, 2)
l = 0
for (i in 1:(nr - 1)) for (j in (i + 1):nr) if (Joint[i, j] == 1) {
l = l + 1
Borders[l, ] = sort(intersect(IndReal[i, ], IndReal[j, ]))
}
Borders = rbind(Borders, IndDummy)
nborders = sum(sum(Joint))/2 + ndummy
}else if((ngrup == 3) | (nr == 1)){
Borders = matrix(0,3,2)
Borders[1,] = c(IndReal[1,1],IndReal[1,2])
Borders[2,] = c(IndReal[1,1],IndReal[1,3])
Borders[3,] = c(IndReal[1,2],IndReal[1,3])
nborders = 3
}
hc = HideCategories(IndReal,ngrup,nreal)
hc_tras = t(hc)
nhide = sum(hc_tras)
ngrupvisible = ngrup - nhide
equivRegiones = matrix(0,2,ngrup)
contVisible = 1
for(i in 1:ngrup){
equivRegiones[1,i] = i
if(hc[i,1] == 0){
equivRegiones[2,i] = contVisible
contVisible = contVisible + 1
}
}
BordersWH = Borders
for(i in 1:nborders){
l = Borders[i, 1]
m = Borders[i, 2]
ind_l = which(equivRegiones[1,]==l)
ind_m = which(equivRegiones[1,]==m)
BordersWH[i, 1] = equivRegiones[2,ind_l]
BordersWH[i, 2] = equivRegiones[2,ind_m]
}
Coord = InvertTesselation(beta,Borders,BordersWH,nborders,ngrupvisible)
Centres = matrix(0, ngrupvisible, 2)
for (i in 1:ngrupvisible) {
Centres[i, 1] = Coord[2 * i - 1]
Centres[i, 2] = Coord[2 * i]
}
if(nr==1)
{
bP12 = (IntersectSegments(Coord[1,1],Coord[2,1],Coord[3,1],Coord[4,1],coorreal[1,1],coorreal[1,2],dummy.x[1],dummy.y[1])
| IntersectSegments(Coord[1,1],Coord[2,1],Coord[3,1],Coord[4,1],coorreal[1,1],coorreal[1,2],dummy.x[2],dummy.y[2])
| IntersectSegments(Coord[1,1],Coord[2,1],Coord[3,1],Coord[4,1],coorreal[1,1],coorreal[1,2],dummy.x[3],dummy.y[3]))
bP13 = (IntersectSegments(Coord[1,1],Coord[2,1],Coord[5,1],Coord[6,1],coorreal[1,1],coorreal[1,2],dummy.x[1],dummy.y[1])
| IntersectSegments(Coord[1,1],Coord[2,1],Coord[5,1],Coord[6,1],coorreal[1,1],coorreal[1,2],dummy.x[2],dummy.y[2])
| IntersectSegments(Coord[1,1],Coord[2,1],Coord[5,1],Coord[6,1],coorreal[1,1],coorreal[1,2],dummy.x[3],dummy.y[3]))
bP23 = (IntersectSegments(Coord[3,1],Coord[4,1],Coord[5,1],Coord[6,1],coorreal[1,1],coorreal[1,2],dummy.x[1],dummy.y[1])
| IntersectSegments(Coord[3,1],Coord[4,1],Coord[5,1],Coord[6,1],coorreal[1,1],coorreal[1,2],dummy.x[2],dummy.y[2])
| IntersectSegments(Coord[3,1],Coord[4,1],Coord[5,1],Coord[6,1],coorreal[1,1],coorreal[1,2],dummy.x[3],dummy.y[3]))
if((bP12 == FALSE) | (bP13 == FALSE) | (bP23 == FALSE)){
CentresF = matrix(0, ngrupvisible, 2)
for (i in 1:ngrupvisible) {
CentresF[i,1]= 2*coorreal[1,1] - Centres[i, 1]
CentresF[i,2]= 2*coorreal[1,2] - Centres[i, 2]
}
bP12n = (IntersectSegments(CentresF[1,1],CentresF[1,2],CentresF[2,1],CentresF[2,2],coorreal[1,1],coorreal[1,2],dummy.x[1],dummy.y[1])
| IntersectSegments(CentresF[1,1],CentresF[1,2],CentresF[2,1],CentresF[2,2],coorreal[1,1],coorreal[1,2],dummy.x[2],dummy.y[2])
| IntersectSegments(CentresF[1,1],CentresF[1,2],CentresF[2,1],CentresF[2,2],coorreal[1,1],coorreal[1,2],dummy.x[3],dummy.y[3]))
bP13n = (IntersectSegments(CentresF[1,1],CentresF[1,2],CentresF[3,1],CentresF[3,2],coorreal[1,1],coorreal[1,2],dummy.x[1],dummy.y[1])
| IntersectSegments(CentresF[1,1],CentresF[1,2],CentresF[3,1],CentresF[3,2],coorreal[1,1],coorreal[1,2],dummy.x[2],dummy.y[2])
| IntersectSegments(CentresF[1,1],CentresF[1,2],CentresF[3,1],CentresF[3,2],coorreal[1,1],coorreal[1,2],dummy.x[3],dummy.y[3]))
bP23n = (IntersectSegments(CentresF[2,1],CentresF[2,2],CentresF[3,1],CentresF[3,2],coorreal[1,1],coorreal[1,2],dummy.x[1],dummy.y[1])
| IntersectSegments(CentresF[2,1],CentresF[2,2],CentresF[3,1],CentresF[3,2],coorreal[1,1],coorreal[1,2],dummy.x[2],dummy.y[2])
| IntersectSegments(CentresF[2,1],CentresF[2,2],CentresF[3,1],CentresF[3,2],coorreal[1,1],coorreal[1,2],dummy.x[3],dummy.y[3]))
if((bP12n == TRUE) & (bP13n == TRUE) & (bP23n == TRUE)){
Centres[,1]=CentresF[,1]
Centres[,2]=CentresF[,2]
}
}
if(sqrt((coorreal[1,1]-Centres[i,1])^2+(coorreal[1,2]-Centres[i,2])^2) < 0.3){ #Si la distancia entre el punto real y los de la inversion es pequeĊa
CentresC = matrix(0, ngrupvisible, 2)
for(i in 1:ngrupvisible){
pte = (Centres[i,2]-coorreal[1,2])/(Centres[i,1]-coorreal[1,1])
oo = coorreal[1,2] - coorreal[1,1]*pte
a = pte*pte + 1
b = -2*coorreal[1,1]+2*pte*(oo-coorreal[1,2])
c=coorreal[1,1]^2 + oo^2 + coorreal[1,2]^2 -2*coorreal[1,2]*oo - 1
xs = Eq2gSolve(a,b,c)
inSegment = ((Centres[i,1] > coorreal[1,1])&(Centres[i,1] < xs[1,1]))
if(inSegment == TRUE) {
CentresC[i,1] = xs[1,1]
CentresC[i,2] = xs[1,1]*pte + oo
}else{
CentresC[i,1] = xs[1,2]
CentresC[i,2] = xs[1,2]*pte + oo
}
}
Centres = CentresC
}
}
result = list()
result$x = coorreal[, 1]
result$y = coorreal[, 2]
result$node = rep(TRUE, nr)
result$n1 = n1
result$n2 = n2
result$n3 = n3
result$dummy.x = dummy.x
result$dummy.y = dummy.y
result$ndummy = ndummy
result$IndReal = IndReal
result$Centers = Centres
result$hideCat = hc_tras
result$equivRegiones = equivRegiones
class(result) <- "voronoiprob"
return(result)
}
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.