The first thing to do is to construct the "linear predictors" or "logits" for each category for each prediction. So you have your model equation:
$$\eta_{ir}=\sum_{j=1}^{p}X_{ij}\hat{\beta}_{jr}\;\; (i=1,\dots,m\;\; r=1,\dots,R)$$
Where for notational convenience, the above is to be understood to have $\hat{\beta}_{jR}=\eta_{iR}=0$ (as this is the reference category), and $\hat{\beta}_{jr}=0$ if variable $j$ was not included in the linear predictor for class $r$. So you will have an $m\times R$ matrix of logits. You then exponentiate to form predicted odds ratios and renormalise to form predicted probabilities. Note that the predicted odds ratios can be calculated by a simple matrix operation if your data is sufficiently organised:
$$\bf{O}=\exp(\boldsymbol{\eta})=\exp(\bf{X}\boldsymbol{\beta})$$
of the $m\times p$ prediction matrix $\bf{X}\equiv\it{\{X_{ij}\}}$ with the $p\times R$ estimated co-efficient matrix $\boldsymbol{\beta}\equiv\it{\{\hat{\beta}_{jr}\}}$, and $\exp(.)$ is defined component wise (i.e. not the matrix exponential). The matrix $\bf{O}$ is an "odds ratio" matrix, with the last column should be all ones. If we take the $m\times 1$ vector $\bf{T}=O1_{R}$ This gives the normalisation constant for each prediction "row" of odds ratios. Now create the $(m\times m)$ diagonal matrix defined by $W_{kk}\equiv T_{k}^{-1}$, and the predicted probability matrix is given by:
$$\bf{P}=\bf{W}\bf{O}$$
So in the example you give the matrix $\boldsymbol{\beta}$ would look like this:
$$\begin{array}{c|c}
Int:1 & Int:2 & Int:3 & Int:4 & Int:5 & 0 \\
\hline
Dev:1 & Dev:2 & Dev:3 & Dev:4 & Dev:5 & 0 \\
\hline
Vote:1 & Vote:2 & Vote:3 & Vote:4 & Vote:5 & 0 \\
\hline
\end{array}$$
With the (roughly rounded) values plugged in we get:
$$\boldsymbol{\beta}=\begin{array}{c|c}
3.68 & -0.90 & -1.56 & -1.13 & -1.68 & 0 \\
\hline
-0.18 & 0.06 & -0.04 & 0.08 & -0.07 & 0 \\
\hline
0.02 & 0.02 & 0.04 & -0.04 & -0.01 & 0 \\
\hline
\hline
\end{array}$$
And the $\bf{X}$ matrix would look like this:
$$\bf{X}=\begin{array}{c|c}
1 & 0.72 & 2\\ \hline 1 & 1.02 & 5\\ \hline 1 & 1.02 & 4\\ \hline 1 & 1.20 & 7\\
\hline 1 & 0.72 & 10\\ \hline 1 & 1.20 & 27\\ \hline 1 & 0.56 & 5\\
\hline 1 & 1.20 & 9\\ \hline 1 & 0.72 & 17\\ \hline 1 & 0.56 & 16\\
\hline
\end{array}$$
So some R-code to do this would simply be (with the matrices $\bf{X}$ and $\boldsymbol{\beta}$ defined as above). The main parts are reading in the data, and padding it with 1s and 0s for $\bf{X}$ and $\boldsymbol{\beta}$:
beta<-cbind(as.matrix( structure(c(3.68021133487111,
-0.903496528862169, -1.56339830041814,
-1.13238307296064, -1.67706243532044, -0.177585202845615, 0.0611115470557421,
-0.0458373863009504, 0.0881133593132653, -0.0686190052488972,
0.0163917121907627, 0.0165232098847022, 0.0373815294869855, -0.0353209839724262,
-0.00698911507852077), .Names = c("(Intercept):1", "(Intercept):2",
"(Intercept):3", "(Intercept):4",
"(Intercept):5", "Deviance:1",
"Deviance:2", "Deviance:3",
"Deviance:4", "Deviance:5", "Votes:1",
"Votes:2", "Votes:3", "Votes:4",
"Votes:5")
, .Dim=c(3L,5L)) ),0)
X<-cbind(1,as.matrix(
structure(c(0.71847390030784,
1.01838748408701, 1.01838748408701,
1.20499277373001, 0.71847390030784, 1.20499277373001, 0.56393315893118,
1.20499277373001, 0.71847390030784, 0.56393315893118, 2, 5, 4, 7, 10, 27, 5, 9, 17, 16), .Dim = c(10L, 2L)) ))
P<-diag(as.vector(exp(X %*% beta) %*%
as.matrix(rep(1,ncol(beta))))^-1) %*%
exp(X %*% beta)
? cor.test() says:
For Spearman's test, p-values are computed using algorithm AS 89 for n < 1290 and exact = TRUE, otherwise via the asymptotic t approximation.
(and default value is exact = TRUE). And its Warning message:
In cor.test.default(x2, y2, method = "spearman") :
Cannot compute exact p-value with ties
They mean that if you don't give an argument exact, cor.test() uses algorithm AS 89 when n < 1290 and no tie. rcorr() always calculates p-values via the asymptotic t approximation, so cor.test() and rcorr() return the same p-value when any ties are present.
getS3method("cor.test", "default") gives you the code. D. J. Best and D. E. Roberts, Algorithm AS 89: The Upper Tail Probabilities of Spearman's Rho, Applied Statistics, Vol.24, No.3 1975, pp.377-379 gives the details about AS 89.
[Edited]
Algorithm AS 89 has two methods, the "exact" method and the "semi-exact" method. cor.test() uses the former when n < 10 and the latter when 9 < n < 1290. In your case, p.value was calculated by the "exact" method.
If you do it manually;
library(e1071) # to get permutations
x2 <- c(1,2,3,4,5,6,7)
y2 <- c(1,3,6,2,7,4,9) # second y2
permutation <- permutations(7)
n <- length(x2)
rho <- cor.test(x2, y2, method="spearman")$estimate # 0.75
# the function to calculate rho between an argument and y2
f.rho <- function(a) 1 - 6 * sum((rank(a) - rank(y2))^2) / (n^3 - n)
x2.all.permutation <- matrix(x2[permutation], ncol=7)
x2.all.permutation.rho <- apply(x2.all.permutation, 1, function(a) f.rho(a))
sum(x2.all.permutation.rho > rho) / factorial(7) * 2 # [1] 0.06626984
cor.test(x2, y2, method="spearman")$p.value # [1] 0.06626984 # the same
Best Answer
Your are looking for
assign().gives
and
Update
I agree that using loops is (very often) bad R coding style (see discussion above). Using
list2env()(thanks to @mbq for mentioning it), this is another solution to @Han Lin Shang's question: