[Rsiena-commits] r217 - in pkg/RSienaTest: . R inst/doc man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Mon Jun 11 03:11:39 CEST 2012
Author: jalospinoso
Date: 2012-06-11 03:11:37 +0200 (Mon, 11 Jun 2012)
New Revision: 217
Modified:
pkg/RSienaTest/DESCRIPTION
pkg/RSienaTest/NAMESPACE
pkg/RSienaTest/R/sienaGOF.r
pkg/RSienaTest/inst/doc/RSiena_Manual.tex
pkg/RSienaTest/man/sienaGOF.Rd
Log:
Adding flexibility to sienaGOF for multiple dependent terms.
Modified: pkg/RSienaTest/DESCRIPTION
===================================================================
--- pkg/RSienaTest/DESCRIPTION 2012-06-07 08:17:43 UTC (rev 216)
+++ pkg/RSienaTest/DESCRIPTION 2012-06-11 01:11:37 UTC (rev 217)
@@ -1,7 +1,7 @@
Package: RSienaTest
Type: Package
Title: Siena - Simulation Investigation for Empirical Network Analysis
-Version: 1.1-216
+Version: 1.1-217
Date: 2012-06-07
Author: Various
Depends: R (>= 2.10.0)
Modified: pkg/RSienaTest/NAMESPACE
===================================================================
--- pkg/RSienaTest/NAMESPACE 2012-06-07 08:17:43 UTC (rev 216)
+++ pkg/RSienaTest/NAMESPACE 2012-06-11 01:11:37 UTC (rev 217)
@@ -6,9 +6,7 @@
varCovar, varDyadCovar, setEffect, includeEffects, includeInteraction,
effectsDocumentation, sienaDataConstraint, print.xtable.sienaFit,
installGui, siena08, iwlsm, sienaTimeTest, includeTimeDummy, sienaGOF,
- sparseMatrixExtraction, snaEdgelistExtraction, snaSociomatrixExtraction,
- igraphEdgelistExtraction, OutdegreeDistribution, IndegreeDistribution,
- GeodesicDistribution, TriadCensus, KnnDistribution, xtable, algorithms,
+ xtable, algorithms,
profileLikelihoods, siena.table)
import(Matrix)
Modified: pkg/RSienaTest/R/sienaGOF.r
===================================================================
--- pkg/RSienaTest/R/sienaGOF.r 2012-06-07 08:17:43 UTC (rev 216)
+++ pkg/RSienaTest/R/sienaGOF.r 2012-06-11 01:11:37 UTC (rev 217)
@@ -1,768 +1,618 @@
-## /*****************************************************************************
-## * SIENA: Simulation Investigation for Empirical Network Analysis
-## *
-## * Web: http://www.stats.ox.ac.uk/~snidjers/siena
-## *
-## * File: gof.r
-## *
-## * Description: This file contains the code to assess goodness of fit
-## *
-## ****************************************************************************/
-##@sienaGOF siena07 Does test for goodness of fit
-sienaGOF <- function(
- sienaFitObject,
- auxiliaryFunction,
- groupName=NULL,
- varName=NULL,
- wave=NULL,
- verbose=FALSE, join=TRUE,
- twoTailed=FALSE,
- cluster=NULL, robust=FALSE, ...)
- {
- require(MASS)
- ##require(Matrix)
- ## Check input
- if (! sienaFitObject$returnDeps)
- {
- stop("You must instruct siena07 to return the simulated networks")
- }
- iterations <- length(sienaFitObject$sims)
- if (iterations < 1)
- {
- stop("You need at least one iteration.")
- }
- groups <- length(sienaFitObject$f$groupNames)
- if (verbose) cat("Detected", iterations, "iterations and", groups,
- "groups.\n")
- if (! is.null(groupName) ) {
- groupNumber <- match( groupName, sienaFitObject$f$groupNames)
- }
- else
- {
- groupName <- sienaFitObject$f$groupNames[1]
- groupNumber <- 1
- }
- if (is.na(groupNumber))
- {
- stop("Invalid group name.")
- }
- if (verbose) cat("Group", groupName, "corresponds to index",
- groupNumber,".\n")
- if(! is.null(varName)) {
- varNumber <- match( varName, sienaFitObject$f$depNames )
- }
- else
- {
- varNumber<-1
- varName <- sienaFitObject$f$depNames[1]
- }
-
- if (varNumber < 1 || varNumber > length(sienaFitObject$f$depNames))
- {
- stop("Invalid variable number -- out of bounds.")
- }
- if (verbose)
- {
- cat("Variable", varName,
- "corresponds to index",varNumber,".\n")
- }
- if (is.null(wave) )
- {
- wave <- 1:(attr(sienaFitObject$f[[groupName]]$depvars[[varName]],
- "netdims")[3] - 1)
- }
- if (varNumber < 1 || varNumber >
- length(sienaFitObject$sims[[1]][[groupNumber]]))
- {
- stop("Invalid variable number -- out of bounds.")
- }
- if (min(wave) < 1 || max(wave) >
- attr(sienaFitObject$f[[groupName]]$depvars[[varName]],
- "netdims")[3] - 1)
- {
- stop("Invalid wave index -- out of bounds")
- }
-
- obsStatsByWave <- lapply(wave, function (j) {
- matrix(
- auxiliaryFunction(NULL,
- sienaFitObject$f, sienaFitObject$sims,
- groupName, varName, j)
- , nrow=1)
- })
- if (join)
- {
- obsStats <- Reduce("+", obsStatsByWave)
- obsStats <- list(Joint=obsStats)
- }
- else
- {
- obsStats <- obsStatsByWave
- names(obsStats) <- paste("Wave",wave)
- }
- class(obsStats) <- "observedAuxiliaryStatistics"
- attr(obsStats,"auxiliaryStatisticName") <-
- deparse(substitute(auxiliaryFunction))
- attr(obsStats,"joint") <- join
-
- ## Calculate the simulated auxiliary statistics
- if (verbose) cat("Calculating auxiliary statistics for waves",wave,".\n")
-
- if (!is.null(cluster)) {
- ttcSimulation <- system.time(simStatsByWave <- lapply(wave,
- function (j) {
- simStatsByWave <- parSapply(cluster, 1:iterations,
- function (i)
- { auxiliaryFunction(i,
- sienaFitObject$f,
- sienaFitObject$sims,
- groupName, varName, j)
- if (verbose && (i %% 100 == 0) )
- cat(" > Completed ", i,
- " calculations\n")
- })
- simStatsByWave <- matrix(simStatsByWave,
- ncol=iterations)
- dimnames(simStatsByWave)[[2]] <- 1:iterations
- t(simStatsByWave)
- }))
- }
- else
- {
- ttcSimulation <- system.time(simStatsByWave <- lapply(wave,
- function (j) {
- simStatsByWave <- sapply(1:iterations, function (i)
- {
- if (verbose && (i %% 100 == 0) )
- {
- cat(" > Completed ", i,
- " calculations\n")
- }
- auxiliaryFunction(i,
- sienaFitObject$f,
- sienaFitObject$sims,
- groupName, varName, j)
- })
- simStatsByWave <-
- matrix(simStatsByWave, ncol=iterations)
- dimnames(simStatsByWave)[[2]] <- 1:iterations
- t(simStatsByWave)
- }))
- }
-
- ## Aggregate by wave if necessary to produce simStats
- if (join)
- {
- simStats <- Reduce("+", simStatsByWave)
- simStats <- list(Joint=simStats)
- }
- else
- {
- simStats <- simStatsByWave
- names(simStats) <- paste("Wave",wave)
- }
- class(simStats) <- "simulatedAuxiliaryStatistics"
- attr(simStats,"auxiliaryStatisticName") <-
- deparse(substitute(auxiliaryFunction))
- attr(simStats,"joint") <- join
- attr(simStats,"time") <- ttcSimulation
-
- applyTest <- function (observed, simulated)
- {
- if (class(simulated) != "matrix")
- {
- stop("Invalid input.")
- }
- if (class(observed) != "matrix")
- {
- observed <- matrix(observed,nrow=1)
- }
- if (class(observed) != "matrix")
- {
- stop("Observation must be a matrix.")
- }
- if (ncol(observed) != ncol(simulated))
- {
- stop("Dimensionality of function parameters do not match.")
- }
- observations <- nrow(observed)
- # simulations<-nrow(simulated)
- variates<-ncol(simulated)
- if (robust) {
- a <- cov.rob(simulated)$cov
- }
- else
- {
- a <- cov(simulated)
- }
- ainv <- ginv(a)
- arank <- rankMatrix(a)
- expectation <- colMeans(simulated);
- centeredSimulations <- scale(simulated, scale=FALSE)
- if (variates==1)
- {
- centeredSimulations <- t(centeredSimulations)
- }
- mhd <- function(x)
- {
- x %*% ainv %*% x
- }
- simTestStat <- apply(centeredSimulations, 1, mhd)
- centeredObservations <- observed - expectation
- obsTestStat <- apply(centeredObservations, 1, mhd)
- if (twoTailed)
- {
- p <- sapply(1:observations, function (i)
- 1 - abs(1 - 2 * sum(obsTestStat[i] <=
- simTestStat)/length(simTestStat)) )
- }
- else
- {
- p <- sapply(1:observations, function (i)
- sum(obsTestStat[i] <= simTestStat) /length(simTestStat))
- }
-
- ret <- list( p = p,
- SimulatedTestStat=simTestStat,
- ObservedTestStat=obsTestStat,
- TwoTailed=twoTailed,
- Simulations=simulated,
- Observations=observed,
- InvCovSimStats=a,
- Rank=arank)
- class(ret) <- "sienaGofTest"
- attr(ret,"auxiliaryStatisticName") <-
- attr(obsStats,"auxiliaryStatisticName")
- ret
- }
-
- res <- lapply(1:length(simStats),
- function (i) {
- applyTest(obsStats[[i]], simStats[[i]]) })
- mhdTemplate <- rep(0, sum(sienaFitObject$test))
- names(mhdTemplate) <- rep(0, sum(sienaFitObject$test))
-
- JoinedOneStepMHD <- mhdTemplate
- JoinedPartialOneStepMHD <- mhdTemplate
- JoinedGmmMhdValue <- mhdTemplate
-
- OneStepMHD <- lapply(wave, function(i) (mhdTemplate))
- PartialOneStepMHD <- lapply(wave, function(i) (mhdTemplate))
- GmmMhdValue <- lapply(wave, function(i) (mhdTemplate))
-
- obsMhd <- NULL
-
- ExpStat <- lapply(wave, function(i) {
- colMeans(simStatsByWave[[i]])
- })
- OneStepSpecs <- matrix(0, ncol=sum(sienaFitObject$test),
- nrow=length(sienaFitObject$theta))
- PartialOneStepSpecs <- matrix(0, ncol=sum(sienaFitObject$test),
- nrow=length(sienaFitObject$theta))
- GmmOneStepSpecs <- matrix(0, ncol=sum(sienaFitObject$test),
- nrow=length(sienaFitObject$theta))
- if (robust) {
- covInvByWave <- lapply(wave, function(i) ginv(
- cov.rob(simStatsByWave[[i]]) ))
- }
- else
- {
- covInvByWave <- lapply(wave, function(i) ginv(
- cov(simStatsByWave[[i]]) ))
- }
-
- obsMhd <- sapply(wave, function (i) {
- (obsStatsByWave[[i]] - ExpStat[[i]]) %*%
- covInvByWave[[i]] %*%
- t(obsStatsByWave[[i]] - ExpStat[[i]] )
- })
-
- if (sum(sienaFitObject$test) > 0) {
- effectsObject <- sienaFitObject$requestedEffects
- nSims <- sienaFitObject$Phase3nits
- for (i in wave) {
- names(OneStepMHD[[i]]) <-
- effectsObject$effectName[sienaFitObject$test]
- names(PartialOneStepMHD[[i]]) <-
- effectsObject$effectName[sienaFitObject$test]
- names(GmmMhdValue[[i]]) <-
- effectsObject$effectName[sienaFitObject$test]
- }
- names(JoinedOneStepMHD) <-
- effectsObject$effectName[sienaFitObject$test]
- names(JoinedPartialOneStepMHD) <-
- effectsObject$effectName[sienaFitObject$test]
- names(JoinedGmmMhdValue) <-
- effectsObject$effectName[sienaFitObject$test]
-
- rownames(OneStepSpecs) <- effectsObject$effectName
- colnames(OneStepSpecs) <- effectsObject$effectName[sienaFitObject$test]
- rownames(PartialOneStepSpecs) <- effectsObject$effectName
- colnames(PartialOneStepSpecs) <-
- effectsObject$effectName[sienaFitObject$test]
- rownames(GmmOneStepSpecs) <- effectsObject$effectName
- colnames(GmmOneStepSpecs) <-
- effectsObject$effectName[sienaFitObject$test]
-
- counterTestEffects <- 0
- for(index in which(sienaFitObject$test)) {
- if (verbose) {
- cat("Estimating test statistic for model including ",
- effectsObject$effectName[index], "\n")
- }
- counterTestEffects <- counterTestEffects + 1
- effectsToInclude <- !sienaFitObject$test
- effectsToInclude[index] <- TRUE
- theta0 <- sienaFitObject$theta
- names(theta0) <- effectsObject$effectName
- theta0 <- theta0[effectsToInclude]
- obsSuffStats <-
- t(sienaFitObject$targets2[effectsToInclude, , drop=FALSE])
- G <- sienaFitObject$sf2[, , effectsToInclude, drop=FALSE] -
- rep(obsSuffStats, each=nSims)
- sigma <- cov(apply(G, c(1, 3), sum))
- SF <- sienaFitObject$ssc[ , , effectsToInclude, drop=FALSE]
- dimnames(SF)[[3]] <- effectsObject$effectName[effectsToInclude]
- dimnames(G) <- dimnames(SF)
- if (!(sienaFitObject$maxlike || sienaFitObject$FinDiff.method))
- {
- D <- RSienaTest:::derivativeFromScoresAndDeviations(SF, G)
- }
- else
- {
- DF <- sienaFitObject$
- sdf2[ , , effectsToInclude, effectsToInclude,
- drop=FALSE]
- D <- t(apply(DF, c(3, 4), mean))
- }
- fra <- apply(G, 3, sum) / nSims
- doTests <- rep(FALSE, sum(effectsToInclude))
- names(doTests) <- effectsObject$effectName[effectsToInclude]
- doTests[effectsObject$effectName[index]] <- TRUE
- mmThetaDelta <- as.numeric(
- RSienaTest:::ScoreTest(length(doTests), D,
- sigma, fra, doTests,
- maxlike=sienaFitObject$maxlike)$oneStep )
- mmPartialThetaDelta <- rep(0,length(theta0))
- mmPartialThetaDelta[length(theta0)] <- mmThetaDelta[length(theta0)]
- JacobianExpStat <- lapply(wave, function (i) {
- t(SF[,i,]) %*% simStatsByWave[[i]]/ nSims })
- Gradient <- lapply(wave, function(i) {
- -2 * JacobianExpStat[[i]] %*%
- covInvByWave[[i]] %*%
- t( obsStatsByWave[[i]] - ExpStat[[i]] ) })
- Hessian <- lapply(wave, function (i) {
- 2 *
- JacobianExpStat[[i]] %*%
- covInvByWave[[i]] %*%
- t(JacobianExpStat[[i]])
- })
- gmmThetaDelta <- -1 * as.numeric( ginv(Reduce("+", Hessian)) %*%
- Reduce("+", Gradient) )
- OneStepSpecs[effectsToInclude,counterTestEffects] <- theta0 +
- mmThetaDelta
- PartialOneStepSpecs[effectsToInclude,counterTestEffects] <-
- theta0 + mmPartialThetaDelta
- GmmOneStepSpecs[effectsToInclude,counterTestEffects] <- theta0 +
- gmmThetaDelta
- for (i in 1:length(obsMhd)) {
- OneStepMHD[[i]][counterTestEffects] <- as.numeric(
- obsMhd[i] +
- mmThetaDelta %*% Gradient[[i]] + 0.5 *
- mmThetaDelta %*% Hessian[[i]] %*% mmThetaDelta)
- GmmMhdValue[[i]][counterTestEffects] <-
- as.numeric( obsMhd[i] +
- gmmThetaDelta %*%
- Gradient[[i]] + 0.5 *
- gmmThetaDelta %*%
- Hessian[[i]] %*%
- gmmThetaDelta )
- PartialOneStepMHD[[i]][counterTestEffects] <-
- as.numeric(
- obsMhd[i] +
- mmPartialThetaDelta %*%
- Gradient[[i]] +
- 0.5 *
- mmPartialThetaDelta %*%
- Hessian[[i]] %*%
- mmPartialThetaDelta)
- }
- JoinedOneStepMHD[counterTestEffects] <-
- Reduce("+",OneStepMHD)[counterTestEffects]
- JoinedPartialOneStepMHD[counterTestEffects] <-
- Reduce("+",PartialOneStepMHD)[counterTestEffects]
- JoinedGmmMhdValue[counterTestEffects] <-
- Reduce("+",GmmMhdValue)[counterTestEffects]
- }
- }
-
- names(res) <- names(obsStats)
- class(res) <- "sienaGOF"
- attr(res, "scoreTest") <- (sum(sienaFitObject$test) > 0)
- attr(res, "originalMahalanobisDistances") <- obsMhd
- attr(res, "joinedOneStepMahalanobisDistances") <-
- JoinedOneStepMHD
- attr(res, "oneStepSpecs") <- OneStepSpecs
- attr(res, "partialOneStepMahalanobisDistances") <- PartialOneStepMHD
- attr(res, "joinedPartialOneStepMahalanobisDistances") <-
- JoinedPartialOneStepMHD
- attr(res, "partialOneStepSpecs") <- PartialOneStepSpecs
- attr(res, "gmmOneStepSpecs") <- GmmOneStepSpecs
- attr(res, "gmmOneStepMahalanobisDistances") <- GmmMhdValue
- attr(res, "joinedGmmOneStepMahalanobisDistances") <- JoinedGmmMhdValue
- attr(res,"auxiliaryStatisticName") <-
- attr(obsStats,"auxiliaryStatisticName")
- attr(res, "simTime") <- attr(simStats,"time")
- attr(res, "twoTailed") <- twoTailed
- attr(res, "joined") <- join
- res
-}
-##@print.sienaGOF siena07 Print method for sienaGOF
-print.sienaGOF <- function (x, ...) {
- ## require(Matrix)
- levels <- 1:length(x)
- pVals <- sapply(levels, function(i) x[[i]]$p)
- titleStr <- "Monte Carlo Mahalanobis distance test P-value: "
- cat("Siena Goodness of Fit (",
- attr(x,"auxiliaryStatisticName") ,")\n=====\n")
- if (! attr(x,"joined"))
- {
- cat(" >",titleStr, "\n")
- for (i in 1:length(pVals))
- {
- cat(" > Wave ", i, ": ", pVals[i], "\n")
- }
- for (i in 1:length(pVals))
- {
- cat(" * Note for wave ",i,
- ": Only ", x[[i]]$Rank, " statistics are ",
- "necessary in the auxiliary function.\n")
- }
- }
- else
- {
- cat(titleStr,pVals[1], "\n")
- cat("**Note: Only ", x[[1]]$Rank, " statistics are ",
- "necessary in the auxiliary function.\n")
- }
-
- if ( attr(x, "twoTailed") )
- {
- cat("-----\nTwo tailed test used.")
- }
- else
- {
- cat("-----\nOne tailed test used ",
- "(i.e. area under curve for greater distance than observation).")
- }
- originalMhd <- attr(x, "originalMahalanobisDistances")
- if (attr(x, "joined")) {
- cat("-----\nCalculated joint MHD = (",
- sum(originalMhd),") for current model.\n")
- }
- else
- {
- for (j in 1:length(originalMhd)) {
- cat("-----\nCalculated wave ", j, " MHD = (",
- originalMhd[j],") for current model.")
- }
- }
-}
-##@summary.sienaGOF siena07 Summary method for sienaGOF
-summary.sienaGOF <- function(object, ...) {
- x <- object
- print(x)
- if (attr(x, "scoreTest")) {
- oneStepSpecs <- attr(x, "oneStepSpecs")
- oneStepMhd <- attr(x, "oneStepMahalanobisDistances")
- gmmMhd <- attr(x, "gmmOneStepMahalanobisDistances")
- gmmOneStepSpecs <- attr(x, "gmmOneStepSpecs")
- partialOneStepSpecs <- attr(x, "partialOneStepSpecs")
- partialOneStepMhd <- attr(x, "partialOneStepMahalanobisDistances")
- joinedPartialOneStepMhd <-
- attr(x, "joinedPartialOneStepMahalanobisDistances")
- joinedOneStepMhd <- attr(x, "joinedOneStepMahalanobisDistances")
- joinedGmmMhd <- attr(x, "joinedGmmOneStepMahalanobisDistances")
- if (attr(x, "joined")) {
- for (i in 1:ncol(oneStepSpecs)) {
- a <- cbind(oneStepSpecs[,i, drop=FALSE],
- partialOneStepSpecs[,i, drop=FALSE],
- gmmOneStepSpecs[,i, drop=FALSE] )
- b <- matrix( c(joinedOneStepMhd[i],
- joinedPartialOneStepMhd[i],
- joinedGmmMhd[i]), ncol=3)
- rownames(b) <- c("MHD")
- a <- rbind(a, b)
- a <- round(a, 3)
- cat("\n**Model", colnames(a)[1], "\n")
- colnames(a) <- c("MM(Full)","MM(Par.)", "GMM")
- print(a)
- }
- }
- else
- {
- for (j in 1:length(oneStepMhd)) {
- for (i in 1:ncol(oneStepSpecs)) {
- a <- cbind(oneStepSpecs[,i, drop=FALSE],
- partialOneStepSpecs[,i, drop=FALSE],
- gmmOneStepSpecs[,i, drop=FALSE] )
- b <- matrix( c(oneStepMhd[[j]][i],
- partialOneStepMhd[[j]][i],
- gmmMhd[[j]][i]), ncol=3)
- rownames(b) <- c("MHD")
- a <- rbind(a, b)
- a <- round(a, 3)
- cat("\n**Model", colnames(a)[1], "\n")
- colnames(a) <- c("MM(Full)","MM(Par.)", "GMM")
- print(a)
- }
- }
- }
- cat("\n-----")
- }
- cat("\nComputation time for auxiliary statistic calculations on simulations: ",
- attr(x, "simTime")["elapsed"] , "seconds.")
-}
-##@plot.sienaGOF siena07 Plot method for sienaGOF
-plot.sienaGOF <- function (x, center=FALSE, scale=FALSE, violin=TRUE,
- key=NULL, perc=.05, wave=1, main=main, ylab=ylab, ...)
-{
- require(lattice)
- args <- list(...)
- if (is.null(args$main))
- {
- main=paste("Goodness of Fit of",
- attr(x,"auxiliaryStatisticName"))
- if (!attr(x,"joined"))
- {
- main = paste(main, "Wave", wave)
- }
- }
- else
- {
- main=args
- }
-
- x <- x[[wave]]
- sims <- x$Simulations
- obs <- x$Observations
- itns <- nrow(sims)
-# vars <- ncol(sims)
- ## Need to check for useless statistics here:
- n.obs <- nrow(obs)
-
- if (center)
- {
- sims.median <- apply(sims, 2, median)
- sims <- sapply(1:ncol(sims), function(i)
- (sims[,i] - sims.median[i]) )
- obs <- matrix(sapply(1:ncol(sims), function(i)
- (obs[,i] - sims.median[i])), nrow=n.obs )
- }
- if (scale)
- {
- sims.min <- apply(sims, 2, min)
- sims.max <- apply(sims, 2, max)
- sims <- sapply(1:ncol(sims), function(i) sims[,i]/(sims.max[i] -
- sims.min[i] ) )
- obs <- matrix(sapply(1:ncol(sims), function(i) obs[,i] /(sims.max[i] -
- sims.min[i] )
- ), nrow=n.obs )
- }
-
- if (is.null(args$ylab))
- {
- ylabel = "Statistic"
- if (center && scale) {
- ylabel = "Statistic (centered and scaled)"
- }
- else if (scale)
- {
- ylabel = "Statistic (scaled)"
- }
- else if (center)
- {
- ylabel = "Statistic (center)"
- }
- else
- {
- ylabel = "Statistic"
- }
- }
- else
- {
- ylabel = args$ylab
- }
-
- screen <- sapply(1:ncol(obs),function(i){
- (sum(is.nan(rbind(sims,obs)[,i])) == 0) }) &
- (diag(var(rbind(sims,obs)))!=0)
- sims <- sims[,screen, drop=FALSE]
- obs <- obs[,screen, drop=FALSE]
- obsLabels <- round(x$Observations[,screen, drop=FALSE],3)
- key <- key[screen]
-
- if (is.null(args$xlab))
- {
- xlabel = paste( paste("p:", round(x$p, 3),
- collapse = " "), collapse = "\n")
- }
- else
- {
- xlabel = args$xlab
- }
-
- xAxis <- (1:sum(screen))
-
- if (!is.null(key))
- {
- if (length(key) != ncol(obs))
- {
- stop("Key length does not match the number of variates.")
- }
- }
- else
- {
- key=xAxis
- }
-
- br <- trellis.par.get("box.rectangle")
- br$col <- 1
- trellis.par.set("box.rectangle", br)
- bu <- trellis.par.get("box.umbrella")
- bu$col <- 1
- trellis.par.set("box.umbrella", bu)
- plot.symbol <- trellis.par.get("plot.symbol")
- plot.symbol$col <- "black"
- plot.symbol$pch <- 4
- plot.symbol$cex <- 1
- trellis.par.set("plot.symbol", plot.symbol)
-
- panelFunction <- function(..., x=x, y=y, box.ratio){
- ind.lower = max( round(itns * perc/2), 1)
- ind.upper = round(itns * (1-perc/2))
- yperc.lower = sapply(1:ncol(sims), function(i)
- sort(sims[,i])[ind.lower] )
- yperc.upper = sapply(1:ncol(sims), function(i)
- sort(sims[,i])[ind.upper] )
- if (violin) {
- panel.violin(x, y, box.ratio=box.ratio, col = "transparent", ...)
- }
- panel.bwplot(x, y, box.ratio=.1, fill = "gray", ...)
- panel.xyplot(xAxis, yperc.lower, lty=3, col = "gray", lwd=3, type="l",
- ...)
- panel.xyplot(xAxis, yperc.upper, lty=3, col = "gray", lwd=3, type="l",
- ...)
- for(i in 1:nrow(obs))
- {
- panel.xyplot(xAxis, obs[i,], col="red", type="l", lwd=1, ...)
- panel.xyplot(xAxis, obs[i,], col="red", type="p", lwd=3, pch=19,
- ...)
- panel.text(xAxis, obs[i,], labels=obsLabels[i,], pos=4)
- }
- }
-
- bwplot(as.numeric(sims)~rep(xAxis, each=itns), horizontal=FALSE,
- panel = panelFunction, xlab=xlabel, ylab=ylabel,
- scales=list(x=list(labels=key), y=list(draw=FALSE)),
- main=main, ...)
-}
-
-sparseMatrixExtraction <- function (i, data, sims, groupName, varName, wave) {
- #require(Matrix)
- dimsOfDepVar<-
- attr(data[[groupName]]$depvars[[varName]],
- "netdims")
- missing <- Matrix(is.na(data[[groupName]]$depvars[[varName]][,,wave+1])*1)
- if (is.null(i)) {
- # sienaGOF wants the observation:
- returnValue <- Matrix(data[[groupName]]$depvars[[varName]][,,wave+1])
- returnValue[is.na(returnValue)] <- 0
- }
- else
- {
- #sienaGOF wants the i-th simulation:
- returnValue <- sparseMatrix(
- sims[[i]][[groupName]][[varName]][[wave]][,1],
- sims[[i]][[groupName]][[varName]][[wave]][,2],
- x=sims[[i]][[groupName]][[varName]][[wave]][,3],
- dims=dimsOfDepVar )
- }
- ## Zero missings:
- 1*((returnValue - missing) > 0)
-}
-
-snaEdgelistExtraction <- function (i, data, sims, groupName, varName, wave) {
- require(sna)
- returnValue <- snaSociomatrixExtraction(i, data, sims, groupName, varName,
- wave)
- as.edgelist.sna(returnValue)
-}
-
-snaSociomatrixExtraction <- function (i, data, sims, groupName, varName, wave) {
- require(sna)
- actors <- attr(data[[groupName]]$nets[[varName]][[wave+1]]$mat1,
- "nActors")
- missing <- t(data[[groupName]]$nets[[varName]][[wave+1]]$mat1)
- attr(missing, "n") <- actors
- missing <- 1*is.na( as.sociomatrix.sna( missing ) )
-
- if (is.null(i)) {
- # sienaGOF wants the observation:
- returnValue <- t(data[[groupName]]$nets[[varName]][[wave+1]]$mat1)
- }
- else
- {
- #sienaGOF wants the i-th simulation:
- returnValue <- sims[[i]][[groupName]][[varName]][[wave]]
- }
- attr(returnValue, "n") <- actors
- returnValue <- as.sociomatrix.sna( returnValue )
- returnValue = 1*((returnValue - missing) > 0)
- returnValue
-}
-
-igraphEdgelistExtraction <- function (i, data, sims, groupName, varName, wave) {
- require(igraph)
- returnValue <- snaSociomatrixExtraction(i, data, sims, groupName, varName,
- wave)
- graph.adjacency(returnValue)
-}
-
-OutdegreeDistribution <- function (i, data, sims, groupName, varName, wave,
- levels=0:8, extractor=snaSociomatrixExtraction) {
- x <- extractor(i, data, sims, groupName, varName, wave)
- a <- apply(x, 1, sum)
- sapply(levels, function(i){ sum(a<=i) })
-}
-
-IndegreeDistribution <- function (i, data, sims, groupName, varName, wave,
- levels=0:8, extractor=snaSociomatrixExtraction) {
- x <- extractor(i, data, sims, groupName, varName, wave)
- a <- apply(x, 2, sum)
- sapply(levels, function(i){ sum(a<=i) })
-}
-
-GeodesicDistribution <- function (i, data, sims, groupName, varName, wave,
- extractor=snaEdgelistExtraction, levels=1:8) {
- require(sna)
- x <- extractor(i, data, sims, groupName, varName, wave)
- a <- geodist(x)$gdist
- sapply(levels, function(i){ sum(a<=i) })
-}
-
-TriadCensus <- function (i, data, sims, groupName, varName, wave,
- extractor=snaEdgelistExtraction) {
- require(sna)
- x <- extractor(i, data, sims, groupName, varName, wave)
- triad.census(x)
-}
-
-KnnDistribution <- function (i, data, sims, groupName, varName, wave,
- extractor=igraphEdgelistExtraction, levels=0:25) {
- require(igraph)
- x <- extractor(i, data, sims, groupName, varName, wave)
- a <- graph.knn(x)$knn
- a[is.nan(a)] <- 0
- sapply(levels, function(i){ sum(a<=i) })
-}
+## /*****************************************************************************
+## * SIENA: Simulation Investigation for Empirical Network Analysis
+## *
+## * Web: http://www.stats.ox.ac.uk/~snidjers/siena
+## *
+## * File: gof.r
+## *
+## * Description: This file contains the code to assess goodness of fit
+## *
+## ****************************************************************************/
+
+##@sienaGOF siena07 Does test for goodness of fit
+sienaGOF <- function(
+ sienaFitObject,
+ auxiliaryFunction,
+ wave=NULL,
+ verbose=FALSE, join=TRUE,
+ twoTailed=FALSE,
+ cluster=NULL, robust=FALSE, ...)
+ {
+ require(MASS)
+ ##require(Matrix)
+ ## Check input
+ if (! sienaFitObject$returnDeps)
+ {
+ stop("You must instruct siena07 to return the simulated networks")
+ }
+ iterations <- length(sienaFitObject$sims)
+ if (iterations < 1)
+ {
+ stop("You need at least one iteration.")
+ }
+ groups <- length(sienaFitObject$f$groupNames)
+ if (verbose) cat("Detected", iterations, "iterations and", groups,
+ "groups.\n")
+ if (is.null(wave) )
+ {
+ wave <- 1:(attr(sienaFitObject$f[[1]]$depvars[[1]], "netdims")[3] - 1)
+ }
+
+ obsStatsByWave <- lapply(wave, function (j) {
+ matrix(
+ auxiliaryFunction(NULL,
+ sienaFitObject$f,
+ sienaFitObject$sims, j)
+ , nrow=1)
+ })
+ if (join)
+ {
+ obsStats <- Reduce("+", obsStatsByWave)
+ obsStats <- list(Joint=obsStats)
+ }
+ else
+ {
+ obsStats <- obsStatsByWave
+ names(obsStats) <- paste("Wave", wave)
+ }
+ class(obsStats) <- "observedAuxiliaryStatistics"
+ attr(obsStats,"auxiliaryStatisticName") <-
+ deparse(substitute(auxiliaryFunction))
+ attr(obsStats,"joint") <- join
+
+ ## Calculate the simulated auxiliary statistics
+ if (verbose) cat("Calculating auxiliary statistics for waves", wave, ".\n")
+
+ if (!is.null(cluster)) {
+ ttcSimulation <- system.time(simStatsByWave <- lapply(wave,
+ function (j) {
+ simStatsByWave <- parSapply(cluster, 1:iterations,
+ function (i)
+ { auxiliaryFunction(i,
+ sienaFitObject$f,
+ sienaFitObject$sims, j)
+ if (verbose && (i %% 100 == 0) )
+ cat(" > Completed ", i,
+ " calculations\n")
+ })
+ simStatsByWave <- matrix(simStatsByWave,
+ ncol=iterations)
+ dimnames(simStatsByWave)[[2]] <- 1:iterations
+ t(simStatsByWave)
+ }))
+ }
+ else
+ {
+ ttcSimulation <- system.time( simStatsByWave <- lapply(wave,
+ function (j) {
+ simStatsByWave <- sapply(1:iterations, function (i)
+ {
+ if (verbose && (i %% 100 == 0) )
+ {
+ cat(" > Completed ", i,
+ " calculations\n")
+ }
+ auxiliaryFunction(i,
+ sienaFitObject$f,
+ sienaFitObject$sims, j)
+ })
+ simStatsByWave <-
+ matrix(simStatsByWave, ncol=iterations)
+ dimnames(simStatsByWave)[[2]] <- 1:iterations
+ t(simStatsByWave)
+ })
+ )
+ }
+
+ ## Aggregate by wave if necessary to produce simStats
+ if (join)
+ {
+ simStats <- Reduce("+", simStatsByWave)
+ simStats <- list(Joint=simStats)
+ }
+ else
+ {
+ simStats <- simStatsByWave
+ names(simStats) <- paste("Wave",wave)
+ }
+ class(simStats) <- "simulatedAuxiliaryStatistics"
+ attr(simStats,"auxiliaryStatisticName") <-
+ deparse(substitute(auxiliaryFunction))
+ attr(simStats,"joint") <- join
+ attr(simStats,"time") <- ttcSimulation
+
+ applyTest <- function (observed, simulated)
+ {
+ if (class(simulated) != "matrix")
+ {
+ stop("Invalid input.")
+ }
+ if (class(observed) != "matrix")
+ {
+ observed <- matrix(observed,nrow=1)
+ }
+ if (class(observed) != "matrix")
+ {
+ stop("Observation must be a matrix.")
+ }
+ if (ncol(observed) != ncol(simulated))
+ {
+ stop("Dimensionality of function parameters do not match.")
+ }
+ observations <- nrow(observed)
+ # simulations<-nrow(simulated)
+ variates<-ncol(simulated)
+ if (robust) {
+ a <- cov.rob(simulated)$cov
+ }
+ else
+ {
+ a <- cov(simulated)
+ }
+ ainv <- ginv(a)
+ arank <- rankMatrix(a)
+ expectation <- colMeans(simulated);
+ centeredSimulations <- scale(simulated, scale=FALSE)
+ if (variates==1)
+ {
+ centeredSimulations <- t(centeredSimulations)
+ }
+ mhd <- function(x)
+ {
+ x %*% ainv %*% x
+ }
+ simTestStat <- apply(centeredSimulations, 1, mhd)
+ centeredObservations <- observed - expectation
+ obsTestStat <- apply(centeredObservations, 1, mhd)
+ if (twoTailed)
+ {
+ p <- sapply(1:observations, function (i)
+ 1 - abs(1 - 2 * sum(obsTestStat[i] <=
+ simTestStat)/length(simTestStat)) )
+ }
+ else
+ {
+ p <- sapply(1:observations, function (i)
+ sum(obsTestStat[i] <= simTestStat) /length(simTestStat))
+ }
+
+ ret <- list( p = p,
+ SimulatedTestStat=simTestStat,
+ ObservedTestStat=obsTestStat,
+ TwoTailed=twoTailed,
+ Simulations=simulated,
+ Observations=observed,
+ InvCovSimStats=a,
+ Rank=arank)
+ class(ret) <- "sienaGofTest"
+ attr(ret,"auxiliaryStatisticName") <-
+ attr(obsStats,"auxiliaryStatisticName")
+ ret
+ }
+
+ res <- lapply(1:length(simStats),
+ function (i) {
+ applyTest(obsStats[[i]], simStats[[i]]) })
+ mhdTemplate <- rep(0, sum(sienaFitObject$test))
+ names(mhdTemplate) <- rep(0, sum(sienaFitObject$test))
+
+ JoinedOneStepMHD <- mhdTemplate
+ JoinedPartialOneStepMHD <- mhdTemplate
+ JoinedGmmMhdValue <- mhdTemplate
+
+ OneStepMHD <- lapply(wave, function(i) (mhdTemplate))
+ PartialOneStepMHD <- lapply(wave, function(i) (mhdTemplate))
+ GmmMhdValue <- lapply(wave, function(i) (mhdTemplate))
+
+ obsMhd <- NULL
+
+ ExpStat <- lapply(wave, function(i) {
+ colMeans(simStatsByWave[[i]])
+ })
+ OneStepSpecs <- matrix(0, ncol=sum(sienaFitObject$test),
+ nrow=length(sienaFitObject$theta))
+ PartialOneStepSpecs <- matrix(0, ncol=sum(sienaFitObject$test),
+ nrow=length(sienaFitObject$theta))
+ GmmOneStepSpecs <- matrix(0, ncol=sum(sienaFitObject$test),
+ nrow=length(sienaFitObject$theta))
+ if (robust) {
+ covInvByWave <- lapply(wave, function(i) ginv(
+ cov.rob(simStatsByWave[[i]]) ))
+ }
+ else
+ {
+ covInvByWave <- lapply(wave, function(i) ginv(
+ cov(simStatsByWave[[i]]) ))
+ }
+
+ obsMhd <- sapply(wave, function (i) {
+ (obsStatsByWave[[i]] - ExpStat[[i]]) %*%
+ covInvByWave[[i]] %*%
+ t(obsStatsByWave[[i]] - ExpStat[[i]] )
+ })
+
+ if (sum(sienaFitObject$test) > 0) {
+ effectsObject <- sienaFitObject$requestedEffects
+ nSims <- sienaFitObject$Phase3nits
+ for (i in wave) {
+ names(OneStepMHD[[i]]) <-
+ effectsObject$effectName[sienaFitObject$test]
+ names(PartialOneStepMHD[[i]]) <-
+ effectsObject$effectName[sienaFitObject$test]
+ names(GmmMhdValue[[i]]) <-
+ effectsObject$effectName[sienaFitObject$test]
+ }
+ names(JoinedOneStepMHD) <-
+ effectsObject$effectName[sienaFitObject$test]
+ names(JoinedPartialOneStepMHD) <-
+ effectsObject$effectName[sienaFitObject$test]
+ names(JoinedGmmMhdValue) <-
+ effectsObject$effectName[sienaFitObject$test]
+
+ rownames(OneStepSpecs) <- effectsObject$effectName
+ colnames(OneStepSpecs) <- effectsObject$effectName[sienaFitObject$test]
+ rownames(PartialOneStepSpecs) <- effectsObject$effectName
+ colnames(PartialOneStepSpecs) <-
+ effectsObject$effectName[sienaFitObject$test]
+ rownames(GmmOneStepSpecs) <- effectsObject$effectName
+ colnames(GmmOneStepSpecs) <-
+ effectsObject$effectName[sienaFitObject$test]
+
+ counterTestEffects <- 0
+ for(index in which(sienaFitObject$test)) {
+ if (verbose) {
+ cat("Estimating test statistic for model including ",
+ effectsObject$effectName[index], "\n")
+ }
+ counterTestEffects <- counterTestEffects + 1
+ effectsToInclude <- !sienaFitObject$test
+ effectsToInclude[index] <- TRUE
+ theta0 <- sienaFitObject$theta
+ names(theta0) <- effectsObject$effectName
+ theta0 <- theta0[effectsToInclude]
+ obsSuffStats <-
+ t(sienaFitObject$targets2[effectsToInclude, , drop=FALSE])
+ G <- sienaFitObject$sf2[, , effectsToInclude, drop=FALSE] -
+ rep(obsSuffStats, each=nSims)
+ sigma <- cov(apply(G, c(1, 3), sum))
+ SF <- sienaFitObject$ssc[ , , effectsToInclude, drop=FALSE]
+ dimnames(SF)[[3]] <- effectsObject$effectName[effectsToInclude]
+ dimnames(G) <- dimnames(SF)
+ if (!(sienaFitObject$maxlike || sienaFitObject$FinDiff.method))
+ {
+ D <- RSienaTest:::derivativeFromScoresAndDeviations(SF, G)
+ }
+ else
+ {
+ DF <- sienaFitObject$
+ sdf2[ , , effectsToInclude, effectsToInclude,
+ drop=FALSE]
+ D <- t(apply(DF, c(3, 4), mean))
+ }
+ fra <- apply(G, 3, sum) / nSims
+ doTests <- rep(FALSE, sum(effectsToInclude))
+ names(doTests) <- effectsObject$effectName[effectsToInclude]
+ doTests[effectsObject$effectName[index]] <- TRUE
+ mmThetaDelta <- as.numeric(
+ RSienaTest:::ScoreTest(length(doTests), D,
+ sigma, fra, doTests,
+ maxlike=sienaFitObject$maxlike)$oneStep )
+ mmPartialThetaDelta <- rep(0,length(theta0))
+ mmPartialThetaDelta[length(theta0)] <- mmThetaDelta[length(theta0)]
+ JacobianExpStat <- lapply(wave, function (i) {
+ t(SF[,i,]) %*% simStatsByWave[[i]]/ nSims })
+ Gradient <- lapply(wave, function(i) {
+ -2 * JacobianExpStat[[i]] %*%
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/rsiena -r 217
More information about the Rsiena-commits
mailing list