Nothing
#/******************************************************************************
# * SIENA: Simulation Investigation for Empirical Network Analysis
# *
# * Web: http://www.stats.ox.ac.uk/~snidjers/siena
# *
# * File: Sienatest.r
# *
# * Description: This module contains the function for instability analysis and
# * score tests.
# *
# *****************************************************************************/
##@InstabilityAnalysis siena07 Not currently used
InstabilityAnalysis<- function(z)
{
##I think this is not correct, because of scaling. cond number of var matrix of X
## can be obtained via svd(data) (which is stored in z$sf). Square of ratio
## of smallest to largest singular value.
Report('Instability Analysis\n')
constant <- z$diver
test <- z$test
covtheta <- z$covtheta
covZ <- z$msf
covth <- covtheta[!(test|constant), !(test|constant)]
covth <- MatrixNorm(covth)
eigenv <- eigen(covth,symmetric=TRUE)$values
ma <- max(eigenv)
mi <- min(eigenv)
if (mi!=0)
{
cond.n <- ma / mi
}
Report('Instability analysis\n', lf)
Report('--------------------\n\n', lf)
Report('Variance-covariance matrix of parameter estimates', lf)
##if (global boolean1 )
## Report(' (without coordinates that are kept constant):\n',lf)
##else
Report(c(':\n\nCondition number = ',format(cond.n, width=4, nsmall=4,
digits=1),
' \n\n'),sep='',lf)
Report(c('Eigen Values ', format(eigenv, width=6, nsmall=6, digits=1)), lf)
Report('\n\n',lf)
covZ <- MatrixNorm(covZ)
eigenvZ <-eigen(covZ, symmetric=TRUE)$values
ma <- max(eigenvZ)
mi <- min(eigenvZ)
if (mi!=0)
{
cond.n <- ma / mi
}
Report('Variance-covariance matrix of X',lf)
Report(c(':\n\nCondition number = ', format(cond.n, width=4, nsmall=4,
digits=1),
' \n\n'), sep='', lf)
Report(c('Eigen Values ', format(eigenvZ, width=6, nsmall=6, digits=1)), lf)
Report(c('\n\n',date(),'\n'),sep='',lf)
mysvd <- svd(z$sf)$d
ma <- max(mysvd)
mi <- min(mysvd)
cond.n <- (ma/mi)^2
Report(c(':\n\nCondition number2 = ',
format(cond.n, width=4, nsmall=4, digits=1),
' \n\n'), sep='', lf)
Report(c('Singular Values ', format(mysvd, width=6, nsmall=6, digits=1)), lf)
Report(c('\n\n', date(), '\n'), sep='', lf)
}
##@MatrixNorm siena07 Not currently used. May be incorrect.
MatrixNorm <- function(mat)
{
tmp <- apply(mat, 2, function(x) x / sqrt(crossprod(x)))
##or sweep(mat,2,apply(mat,2,function(x)x/sqrt(crossprod(x))
tmp
}
##@TestOutput siena07 Print report
TestOutput <- function(z, x)
{
testn <- sum(z$test)
## browser()
if (testn)
{
if (x$maxlike)
{
Heading(2, outf,'Score test <c>')
}
else
{
Heading(2, outf, 'Generalised score test <c>')
}
Report('Testing the goodness-of-fit of the model restricted by\n', outf)
j <- 0
for (k in 1:z$pp)
{
if (z$test[k])
{
j <- j + 1
Report(c(" (", j, ") ",
format(paste(z$requestedEffects$type[k], ": ",
z$requestedEffects$effectName[k],
sep=""),
width=50), " = ",
sprintf("%8.4f", z$theta[k]), "\n"),
sep = "", outf)
}
}
Report('_________________________________________________\n',outf)
Report(' ',outf)
Report(' \n',outf)
if (testn > 1)
Report('Joint test:\n-----------\n',outf)
Report(c(' c = ',sprintf("%8.4f", z$testresOverall),
' d.f. = ',j,' p-value '),sep='',outf)
pvalue <- 1-pchisq(z$testresOverall,j)
if (!is.na(pvalue))
{
if (pvalue < 0.0001)
{
Report('< 0.0001',outf)
}
else
{
Report(c('= ',sprintf("%8.4f",pvalue)), sep = '', outf)
}
}
else
{
Report(' NA ',outf)
}
if (testn==1)
Report(c('\n one-sided (normal variate): ',
sprintf("%8.4f",z$testresulto[1])), sep = '', outf)
if (testn> 1)
{
Report('\n\n',outf)
for (k in 1:j)
{
Report(c('(',k,') tested separately:\n'),sep='',outf)
Report('-----------------------\n',outf)
Report(' - two-sided:\n',outf)
Report(c(' c = ', sprintf("%8.4f", z$testresult[k]),
' d.f. = 1 p-value '), sep = '', outf)
pvalue<- 1-pchisq(z$testresult[k],1)
if (!is.na(pvalue))
{
if (pvalue < 0.0001)
{
Report('< 0.0001\n',outf)
}
else
{
Report(c('= ', sprintf("%8.4f", pvalue), '\n'), sep = '',
outf)
}
}
else
{
Report(' NA ',outf)
}
Report(c(' - one-sided (normal variate): ',
sprintf("%8.4f", z$testresulto[k])), sep = '', outf)
if (k < j)
{
Report('\n\n',outf)
}
}
}
Report(' \n_________________________________________________\n\n',
outf)
Report('One-step estimates: \n\n', outf)
for (i in 1 : z$pp)
{
onestepest <- z$oneStep[i] + z$theta[i]
Report(c(format(paste(z$requestedEffects$type[i], ': ',
z$requestedEffects$effectName[i], sep = ''),
width=50),
sprintf("%8.4f", onestepest), '\n'), sep = '', outf)
}
Report('\n',outf)
}
}
##@ScoreTest siena07 Do score tests
ScoreTest<- function(pp, dfra, msf, fra, test, redundant, maxlike)
{
testresult <- rep(NA, pp) ##for chisq per parameter
testresulto <- rep(NA, pp) ##for one-sided tests per parameter
##first the general one
ans <- EvaluateTestStatistic(maxlike, test, redundant, dfra, msf, fra)
testresOverall <- ans$cvalue
covMatrix <- ans$covMatrix
if (sum(test) == 1)
{
testresulto[1] <- ans$oneSided
}
else
{
## single df tests
use <- !test
k <- 0
for (i in 1:pp)
{
if (test[i])
{
k <- k + 1
use[i] <- TRUE
ans <- EvaluateTestStatistic(maxlike, test[use], redundant[use],
dfra[use, use], msf[use, use], fra[use])
testresult[k] <- ans$cvalue
testresulto[k] <- ans$oneSided
use[i] <- FALSE
}
}
}
##onestep estimator
if (maxlike)
dfra2 <- dfra + msf
else
dfra2 <- dfra
if (inherits(try(dinv2 <- solve(dfra2), silent=TRUE), "try-error"))
{
Report("Error message for inversion to get onestep estimator: \n", cf)
dinv2 <- dfra2
dinv2[] <- NA
oneStep <- rep(NA, nrow(dfra2))
}
else
{
oneStep<- -dinv2 %*% fra
}
list(testresult=testresult, testresulto=testresulto,
testresOverall=testresOverall, covMatrix=covMatrix,
oneStep=oneStep, dinv2= dinv2, dfra2=dfra2)
}
##@EvaluateTestStatistic siena07 Calculate score test statistics
EvaluateTestStatistic<- function(maxlike, test, redundant, dfra, msf, fra)
{
##uses local arrays set up in the calling procedure
d11 <- dfra[!(test|redundant), !(test|redundant), drop=FALSE]
d22 <- dfra[test, test, drop=FALSE]
d21 <- dfra[test, !(test|redundant), drop=FALSE]
d12 <- t(d21)
sigma11 <- msf[!(test|redundant), !(test|redundant), drop=FALSE]
sigma22<- msf[test, test,drop=FALSE]
sigma12 <- msf[!(test|redundant), test, drop=FALSE]
sigma21<- t(sigma12)
z1 <- fra[!(test|redundant)]
z2 <- fra[test]
if (inherits(try(id11 <- solve(d11), silent=TRUE), "try-error"))
{
warning('Score test: Error for inversion of d11 \n')
oneSided <- NA
v9 <- d22
v9[] <- NA
cvalue <- matrix(NA, 1, 1)
}
else
{
rg <- d21 %*% id11
if (!maxlike)
{
##orthogonalise deviation vector
ov <- z2 - rg %*% z1
##compute var(ov) = sigma22 - (d21 %*% id11) %*% sigma12 -
## sigma21 %*% t(id11)%*% t(d21) +
## d21%*%id11 %*% sigma11 %*% t(id11) %*% t(d21)
v2 <- sigma21 - rg %*% sigma11
v6 <- v2 %*% t(id11) %*% t(d21)
v9 <- sigma22 - rg %*% sigma12 - v6
}
else
{
ov <- -z2
v9 <- d22 - rg %*% d12
}
if (inherits(try(vav <- solve(v9), silent=TRUE), "try-error"))
## vav is the inverse variance matrix of ov
{
warning('Score test: Error for inversion of v9\n')
vav <- v9
vav[] <- NA
cvalue <- NA
oneSided <- NA
}
else
{
cvalue <- t(ov) %*% vav %*% ov
if (cvalue < 0) cvalue <- 0
if (sum(test) == 1)
{
if (vav > 0)
oneSided <- ov * sqrt(vav)
else
oneSided <- 0
if (!maxlike) oneSided <- - oneSided
## change the sign for intuition for users
}
else
{
oneSided <- 0
}
}
}
list(cvalue=cvalue, oneSided=oneSided, covMatrix=v9)
}
##@scoreTest Calculate score test
score.Test <- function(ans, test=ans$test)
# use: ans must be a sienaFit object;
# test must be a boolean vector with length equal to the number of parameters of ans,
# or a vector of integer numbers between 1 and ans$pp.
{
if ((is.numeric(test)) || (is.integer(test)))
{
if (max(test) > ans$pp)
{
stop(paste('The maximum requested coordinate is too high:',
max(test)))
}
test <- (1:ans$pp) %in% test
}
if (sum(test) <= 0) stop(paste('Something should be tested, but the total requested is',
sum(test)))
if (length(test) != ans$pp) stop('Dimensions of test must agree')
if (any(test & (!ans$fix))) warning('Warning: some tested parameters were not fixed; do you know what you are doing??? \n')
fra <- colMeans(ans$sf, na.rm=TRUE)
redundant <- (ans$fix & (!test))
tests <- EvaluateTestStatistic(ans$maxlike, test, redundant, ans$dfra, ans$msf, fra)
teststat <- tests$cvalue
df <- sum(test)
if (df == 1)
{
onesided <- tests$oneSided
}
else
{
onesided <- NULL
}
pval <- 1 - pchisq(teststat, df)
efnames <- paste(ans$effects$name[test], ans$effects$effectName[test], sep=': ')
if (any(ans$effects$type != "eval"))
{
efnames <- paste(efnames, ans$effects$type[test], sep=' ')
}
t.ans <- list(chisquare=teststat, df=df, pvalue=pval, onesided=onesided, efnames=efnames)
class(t.ans) <- "sienaTest"
t.ans
}
##@Wald.RSiena Calculate Wald test statistics
Wald.RSiena <- function(A, ans)
{
if (is.vector(A)){A <- matrix(A, nrow=1)}
if (dim(A)[2] != ans$pp){stop(paste('A must have', ans$pp, 'columns.'))}
sigma <- ans$covtheta
if (any(is.na(sigma))) {
# happens when some coordinates were fixed
# in the call of siena07 leading to ans;
# then the non-used part of sigma,
# which partially consists of NA,
# is replaced by the identity matrix.
zero.cols <- apply(A, 2, function(colum){all(colum==0)})
sigma[zero.cols, ] <- 0
sigma[, zero.cols] <- 0
diag(sigma)[zero.cols] <- 1
}
th <- A %*% ans$theta
covmat <- A %*% sigma %*% t(A)
chisq <- drop(t(th) %*% solve(covmat) %*% th)
df <- nrow(A)
pval <- 1 - pchisq(chisq, df)
if (df == 1)
{
onesided <- sign(th) * sqrt(chisq)
}
else
{
onesided <- NULL
}
t.ans <- list(chisquare=chisq, df=df, pvalue=pval, onesided=onesided)
class(t.ans) <- "sienaTest"
t.ans
}
##@Multipar.RSiena Calculate Wald test statistic for hypothesis that subvector = 0.
Multipar.RSiena <- function(ans, ...)
{
p <- length(ans$theta)
tested <- c(...)
efnames <- paste(ans$effects$name[tested], ans$effects$effectName[tested], sep=': ')
if (any(ans$effects$type != "eval"))
{
efnames <- paste(efnames, ans$effects$type[tested], sep=' ')
}
k <- length(tested)
A <- matrix(0, nrow=k, ncol=p)
A[cbind(1:k,tested)] <- 1
t.ans <- Wald.RSiena(A, ans)
t.ans$efnames <- efnames
t.ans
}
##@print.sienaTest Methods
print.sienaTest <- function(x, ...)
{
if (!inherits(x, "sienaTest"))
{
stop("not a legitimate sienaTest object")
}
if (any(is.na(unlist(x))))
{
stop("some estimates or standard errors are NA")
}
if (!is.null(x$efnames))
{
cat('Tested effects:\n ')
cat(paste(x$efnames,'\n'))
}
cat(paste('chi-squared = ',
format(round(x$chisquare, digits=2), nsmall = 2),
', d.f. = ', x$df, '; ', sep=''))
if ((x$df == 1) & (!is.null(x$onesided)))
{
cat(paste('one-sided Z = ',
format(round(x$onesided, digits=2), nsmall = 2), '; ', sep=''))
cat('two-sided')
}
if (x$pvalue < 0.001)
{
cat(' p < 0.001. \n')
}
else
{
cat(paste(' p = ',
format(round(x$pvalue, digits=3), nsmall = 2), '. \n', sep=''))
}
invisible(x)
}
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.