[Rsiena-commits] r63 - in pkg/RSienaTest: . R man src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Mon Mar 15 16:02:22 CET 2010
Author: ripleyrm
Date: 2010-03-15 16:02:21 +0100 (Mon, 15 Mar 2010)
New Revision: 63
Added:
pkg/RSienaTest/R/iwlsm.R
pkg/RSienaTest/R/siena08.r
pkg/RSienaTest/man/iwlsm.Rd
pkg/RSienaTest/man/print.sienaMeta.Rd
pkg/RSienaTest/man/siena08.Rd
pkg/RSienaTest/man/summary.iwlsm.Rd
Modified:
pkg/RSienaTest/DESCRIPTION
pkg/RSienaTest/NAMESPACE
pkg/RSienaTest/R/globals.r
pkg/RSienaTest/R/phase1.r
pkg/RSienaTest/R/phase2.r
pkg/RSienaTest/R/phase3.r
pkg/RSienaTest/R/print01Report.r
pkg/RSienaTest/R/robmon.r
pkg/RSienaTest/R/siena01.r
pkg/RSienaTest/R/siena07.r
pkg/RSienaTest/R/sienaDataCreate.r
pkg/RSienaTest/R/sienaModelCreate.r
pkg/RSienaTest/R/simstatsc.r
pkg/RSienaTest/man/siena07.Rd
pkg/RSienaTest/man/sienaModelCreate.Rd
pkg/RSienaTest/src/siena07.cpp
Log:
1. Added siena08 and iwlsm functions. 2. Fixed rlecuyer random numbers to go on to new substream with each iteration. 3. missing values in dyadic covariates (incomplete). 4. New feature in Report: any output file optional. 5. Can now use multiple processors by period rather than by iteration. 6 Few other minor fixes.
Modified: pkg/RSienaTest/DESCRIPTION
===================================================================
--- pkg/RSienaTest/DESCRIPTION 2010-03-15 14:56:31 UTC (rev 62)
+++ pkg/RSienaTest/DESCRIPTION 2010-03-15 15:02:21 UTC (rev 63)
@@ -6,7 +6,7 @@
Author: Various
Depends: R (>= 2.9.0), xtable
Imports: Matrix
-Suggests: tcltk, rlecuyer, snow, network, codetools
+Suggests: tcltk, snow, rlecuyer, network, codetools, lattice
SystemRequirements: GNU make, tcl/tk 8.5, Tktable
Maintainer: Ruth Ripley <ruth at stats.ox.ac.uk>
Description: Fits models to longitudinal networks
Modified: pkg/RSienaTest/NAMESPACE
===================================================================
--- pkg/RSienaTest/NAMESPACE 2010-03-15 14:56:31 UTC (rev 62)
+++ pkg/RSienaTest/NAMESPACE 2010-03-15 15:02:21 UTC (rev 63)
@@ -2,10 +2,10 @@
export(coCovar, coDyadCovar, getEffects, model.create, print01Report,
siena01Gui, siena07, sienaCompositionChange,
sienaCompositionChangeFromFile, sienaDataCreate, sienaDataCreateFromSession,
-sienaGroupCreate, sienaModelCreate, sienaNet, sienaNodeSet, simstats0c,
+sienaGroupCreate, sienaModelCreate, sienaNet, sienaNodeSet, #simstats0c,
varCovar, varDyadCovar, setEffect, includeEffects, includeInteraction,
- effectsDocumentation, sienaDataConstraint, maxlikefn,
- installGui)#, sienaTimeTest)
+ effectsDocumentation, sienaDataConstraint, #maxlikefn,
+ installGui, siena08, iwlsm)#, sienaTimeTest)
import(Matrix)
import(xtable)
@@ -18,6 +18,23 @@
S3method(summary, sienaFit)
S3method(xtable, sienaFit)
S3method(print, xtable.sienaFit)
+S3method(print, sienaMeta)
+S3method(print, summary.sienaMeta)
+S3method(summary, sienaMeta)
+S3method(plot, sienaMeta)
+S3method(predict, iwlsm)
+S3method(print, iwlsm)
+S3method(print, summary.iwlsm)
+S3method(summary, iwlsm)
+S3method(iwlsm, formula)
+S3method(iwlsm, default)
+S3method(se.contrast, iwlsm)
+S3method(vcov, iwlsm)
+S3method(addAttributes, default)
+S3method(addAttributes, coCovar)
+S3method(addAttributes, varCovar)
+S3method(addAttributes, coDyadCovar)
+S3method(addAttributes, varDyadCovar)
#S3method(print, sienaTimeTest)
#S3method(summary, sienaTimeTest)
#S3method(print, summary.sienaTimeTest)
Modified: pkg/RSienaTest/R/globals.r
===================================================================
--- pkg/RSienaTest/R/globals.r 2010-03-15 14:56:31 UTC (rev 62)
+++ pkg/RSienaTest/R/globals.r 2010-03-15 15:02:21 UTC (rev 63)
@@ -23,11 +23,11 @@
x <- x
beverbose <- verbose
besilent <- silent
- function(txt, dest, fill=FALSE, sep=" ", hdest,
- open=FALSE, close=FALSE,
- type=c("a", "w"), projname="Siena" , verbose=FALSE, silent=FALSE)
+ function(txt, dest, fill=FALSE, sep=" ", hdest, openfiles=FALSE,
+ closefiles=FALSE, type=c("a", "w", "n"), projname="Siena" ,
+ verbose=FALSE, silent=FALSE)
{
- if (open)
+ if (openfiles)
{
type <- match.arg(type)
beverbose <<- verbose
@@ -36,15 +36,16 @@
{
x$outf <<- file(paste(projname, ".out", sep=""), open="w")
}
- else
+ else if (type =="a")
{
x$outf <<- file(paste(projname, ".out", sep=""), open="a")
}
}
- else if (close)
+ else if (closefiles)
{
close(x[["outf"]])
+ x$outf <<- NULL
}
else
{
@@ -82,15 +83,23 @@
}
else
{
- cat(txt, file=x[[deparse(substitute(dest))]],
- fill=fill, sep=sep)
+ if (is.null(x[[deparse(substitute(dest))]]))
+ {
+ cat(txt, fill=fill, sep=sep)
+ }
+ else
+ {
+ cat(txt, file=x[[deparse(substitute(dest))]],
+ fill=fill, sep=sep)
+ }
}
}
}
- }
+ }
}
}
+
##@Report Globals
Report <- local({verbose <- NULL; silent <- NULL;
Reportfun(list(outf=outf, lf=lf, cf=cf, bof=bof), verbose,
Added: pkg/RSienaTest/R/iwlsm.R
===================================================================
--- pkg/RSienaTest/R/iwlsm.R (rev 0)
+++ pkg/RSienaTest/R/iwlsm.R 2010-03-15 15:02:21 UTC (rev 63)
@@ -0,0 +1,445 @@
+# file iwlsm.R derived from MASS/R/rlm.R which is
+# copyright (C) 1994-2005 W. N. Venables and B. D. Ripley
+#
+# 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/
+#
+#/******************************************************************************
+# * SIENA: Simulation Investigation for Empirical Network Analysis
+# *
+# * Web: http://www.stats.ox.ac.uk/~snidjers/siena
+# *
+# * File: iwlsm.r
+# *
+# * Description: This module contains code to fit iterated weighted least
+# * squares model
+# *
+# *****************************************************************************/
+##@iwlsm iwlsm
+iwlsm <- function(x, ...) UseMethod("iwlsm")
+##@iwlsm.formula iwlsm
+iwlsm.formula <-
+ function(formula, data, weights, ses, ..., subset, na.action,
+ method = c("M", "MM", "model.frame"),
+ wt.method = c("inv.var", "case"),
+ model = TRUE, x.ret = TRUE, y.ret = FALSE, contrasts = NULL)
+{
+ mf <- match.call(expand.dots = FALSE)
+ mf$method <- mf$wt.method <- mf$model <- mf$x.ret <- mf$y.ret <- mf$contrasts <- mf$... <- NULL
+ mf[[1L]] <- as.name("model.frame")
+ mf <- eval.parent(mf)
+ method <- match.arg(method)
+ wt.method <- match.arg(wt.method)
+ if(method == "model.frame") return(mf)
+ mt <- attr(mf, "terms")
+ y <- model.response(mf)
+ offset <- model.offset(mf)
+ if(!is.null(offset)) y <- y - offset
+ x <- model.matrix(mt, mf, contrasts)
+ xvars <- as.character(attr(mt, "variables"))[-1L]
+ if ((yvar <- attr(mt, "response")) > 0L)
+ xvars <- xvars[-yvar]
+ xlev <- if (length(xvars) > 0L) {
+ xlev <- lapply(mf[xvars], levels)
+ xlev[!sapply(xlev, is.null)]
+ }
+ weights <- model.weights(mf)
+ if(!length(weights)) weights <- rep(1, nrow(x))
+ ses <- mf[["(ses)"]]
+ fit <- iwlsm.default(x, y, weights, method = method,
+ wt.method = wt.method, ses=ses, ...)
+ fit$terms <- mt
+ ## fix up call to refer to the generic, but leave arg name as `formula'
+ cl <- match.call()
+ cl[[1L]] <- as.name("iwlsm")
+ fit$call <- cl
+ fit$contrasts <- attr(x, "contrasts")
+ fit$xlevels <- .getXlevels(mt, mf)
+ fit$na.action <- attr(mf, "na.action")
+ if(model) fit$model <- mf
+ if(!x.ret) fit$x <- NULL
+ if(y.ret) fit$y <- y
+ fit
+}
+##@iwlsm.default iwlsm
+iwlsm.default <-
+ function(x, y, weights, ses, ..., w = rep(1, nrow(x)),
+ init = "ls", psi = psi.iwlsm,
+ scale.est = c("MAD", "Huber", "proposal 2"), k2 = 1.345,
+ method = c("M", "MM"), wt.method = c("inv.var", "case"),
+ maxit = 20, acc = 1e-4, test.vec = "resid", lqs.control=NULL)
+{
+ irls.delta <- function(old, new)
+ sqrt(sum((old - new)^2)/max(1e-20, sum(old^2)))
+ irls.rrxwr <- function(x, w, r)
+ {
+ w <- sqrt(w)
+ max(abs((matrix(r * w, 1L, length(r)) %*% x)/
+ sqrt(matrix(w, 1L, length(r)) %*% (x^2))))/sqrt(sum(w * r^2))
+ }
+ wmad <- function(x, w)
+ {
+ o <- sort.list(abs(x)); x <- abs(x)[o]; w <- w[o]
+ p <- cumsum(w)/sum(w)
+ n <- sum(p < 0.5)
+ if (p[n + 1L] > 0.5) x[n + 1L]/0.6745 else (x[n + 1L] + x[n + 2L])/(2*0.6745)
+ }
+
+ method <- match.arg(method)
+ wt.method <- match.arg(wt.method)
+ nmx <- deparse(substitute(x))
+ if(is.null(dim(x))) {
+ x <- as.matrix(x)
+ colnames(x) <- nmx
+ } else x <- as.matrix(x)
+ if(is.null(colnames(x)))
+ colnames(x) <- paste("X", seq(ncol(x)), sep="")
+ if(qr(x)$rank < ncol(x))
+ stop("'x' is singular: singular fits are not implemented in iwlsm")
+
+ if(!(any(test.vec == c("resid", "coef", "w", "NULL"))
+ || is.null(test.vec))) stop("invalid 'test.vec'")
+ ## deal with weights
+ xx <- x
+ yy <- y
+ if(!missing(weights)) {
+ if(length(weights) != nrow(x))
+ stop("length of 'weights' must equal number of observations")
+ if(any(weights < 0)) stop("negative 'weights' value")
+ if(wt.method == "inv.var") {
+ fac <- sqrt(weights)
+ y <- y*fac; x <- x* fac
+ wt <- NULL
+ } else {
+ w <- w * weights
+ wt <- weights
+ }
+ } else wt <- NULL
+
+ if(method == "M") {
+ scale.est <- match.arg(scale.est)
+ if(!is.function(psi)) psi <- get(psi, mode="function")
+ ## match any ... args to those of psi.
+ arguments <- list(...)
+ if(length(arguments)) {
+ pm <- pmatch(names(arguments), names(formals(psi)), nomatch = 0L)
+ if(any(pm == 0L)) warning("some of ... do not match")
+ pm <- names(arguments)[pm> 0L]
+ ## formals(psi)[pm] <- unlist(arguments[pm])
+ formals(psi)[pm] <- arguments[pm]
+ }
+ if(is.character(init)) {
+ temp <- if(init == "ls") lm.wfit(x, y, w, method="qr")
+ else if(init == "lts") {
+ if(is.null(lqs.control)) lqs.control <- list(nsamp=200L)
+ do.call("lqs", c(list(x, y, intercept = FALSE), lqs.control))
+ } else stop("'init' method is unknown")
+ coef <- temp$coefficient
+ resid <- temp$residuals
+ } else {
+ if(is.list(init)) coef <- init$coef
+ else coef <- init
+ resid <- drop(y - x %*% coef)
+ }
+ } else if(method == "MM") {
+ scale.est <- "MM"
+ temp <- do.call("lqs",
+ c(list(x, y, intercept = FALSE, method = "S",
+ k0 = 1.548), lqs.control))
+ coef <- temp$coefficients
+ resid <- temp$residuals
+ # psi <- psi.bisquare
+ if(length(arguments <- list(...)))
+ if(match("c", names(arguments), nomatch = 0L)) {
+ c0 <- arguments$c
+ if (c0 > 1.548) formals(psi)$c <- c0
+ else
+ warning("'c' must be at least 1.548 and has been ignored")
+ }
+ scale <- temp$scale
+ } else stop("'method' is unknown")
+ environment(psi) <- .GlobalEnv
+ done <- FALSE
+ conv <- NULL
+ n1 <- (if(is.null(wt)) nrow(x) else sum(wt)) - ncol(x)
+ theta <- 2*pnorm(k2) - 1
+ gamma <- theta + k2^2 * (1 - theta) - 2 * k2 * dnorm(k2)
+ ## At this point the residuals are weighted for inv.var and
+ ## unweighted for case weights. Only Huber handles case weights
+ ## correctly.
+ if(scale.est != "MM")
+ scale <- if(is.null(wt)) mad(resid, 0) else wmad(resid, wt)
+ for(iiter in 1L:maxit) {
+ if(!is.null(test.vec)) testpv <- get(test.vec)
+ if(scale.est != "MM") {
+ scale <- if(scale.est == "MAD")
+ if(is.null(wt)) median(abs(resid))/0.6745 else wmad(resid, wt)
+ else if(is.null(wt))
+ sqrt(sum(pmin(resid^2, (k2 * scale)^2))/(n1*gamma))
+ else sqrt(sum(wt*pmin(resid^2, (k2 * scale)^2))/(n1*gamma))
+ if(scale == 0) {
+ done <- TRUE
+ break
+ }
+ }
+ ## w <- psi(resid/scale)###
+ w <- psi(resid, w=w, sj2=ses)
+ if(!is.null(wt)) w <- w * weights
+ temp <- lm.wfit(x, y, w, method="qr")
+ coef <- temp$coefficients
+ resid <- temp$residuals ## w * res* res
+ if(!is.null(test.vec)) convi <- irls.delta(testpv, get(test.vec))
+ else convi <- irls.rrxwr(x, w, resid)
+ conv <- c(conv, convi)
+ done <- (convi <= acc)
+ if(done) break
+ }
+ if(!done)
+ warning(gettextf("iwlsm failed to converge in %d steps", maxit),
+ domain = NA)
+ fitted <- drop(xx %*% coef)
+ ## fix up call to refer to the generic, but leave arg name as `formula'
+ cl <- match.call()
+ cl[[1L]] <- as.name("iwlsm")
+ fit <- list(coefficients = coef, residuals = yy - fitted, wresid = resid,
+ effects = temp$effects,
+ rank = temp$rank, fitted.values = fitted,
+ assign = temp$assign, qr = temp$qr, df.residual = NA, w = w,
+ s = scale, psi = psi, k2 = k2,
+ weights = if(!missing(weights)) weights,
+ conv = conv, converged = done, x = xx, call = cl)
+ class(fit) <- c("iwlsm", "lm")
+ fit
+}
+##@print.iwlsm Methods
+print.iwlsm <- function(x, ...)
+{
+ if(!is.null(cl <- x$call)) {
+ cat("Call:\n")
+ dput(cl, control=NULL)
+ }
+ if(x$converged)
+ cat("Converged in", length(x$conv), "iterations\n")
+ else cat("Ran", length(x$conv), "iterations without convergence\n")
+ coef <- x$coefficients
+ cat("\nCoefficients:\n")
+ print(coef, ...)
+ nobs <- length(x$residuals)
+ rdf <- nobs - length(coef)
+ cat("\nDegrees of freedom:", nobs, "total;", rdf, "residual\n")
+ if(nzchar(mess <- naprint(x$na.action))) cat(" (", mess, ")\n", sep="")
+ cat("Scale estimate:", format(signif(x$s,3)), "\n")
+ invisible(x)
+}
+##@summary.iwlsm Methods
+summary.iwlsm <- function(object, method = c("XtX", "XtWX"),
+ correlation = FALSE, ...)
+{
+ method <- match.arg(method)
+ s <- object$s
+ coef <- object$coefficients
+ ptotal <- length(coef)
+ wresid <- object$wresid
+ res <- object$residuals
+ n <- length(wresid)
+ if(any(na <- is.na(coef))) coef <- coef[!na]
+ cnames <- names(coef)
+ p <- length(coef)
+ rinv <- diag(p)
+ dimnames(rinv) <- list(cnames, cnames)
+ wts <- if(length(object$weights)) object$weights else rep(1, n)
+ if(length(object$call$wt.method) && object$call$wt.method == "case") {
+ rdf <- sum(wts) - p
+ w <- object$psi(wresid/s)
+ S <- sum(wts * (wresid*w)^2)/rdf
+ psiprime <- object$psi(wresid/s, deriv=1)
+ m1 <- sum(wts*psiprime)
+ m2 <- sum(wts*psiprime^2)
+ nn <- sum(wts)
+ mn <- m1/nn
+ kappa <- 1 + p*(m2 - m1^2/nn)/(nn-1)/(nn*mn^2)
+ stddev <- sqrt(S)*(kappa/mn)
+ } else {
+ res <- res * sqrt(wts)
+ rdf <- n - p
+ ## w <- object$psi(wresid/s)
+ w <- object$w
+ ## S <- sum((wresid*w)^2)/rdf
+ ## psiprime <- object$psi(wresid, deriv=1, w=w)
+ ## mn <- mean(psiprime)
+ ## kappa <- 1 + p*var(psiprime)/(n*mn^2)
+ hh <- hatvalues(lm(object$fitted.values ~ object$x,
+ weights=w))
+ S <- sum(w * (wresid^2)) / (1 - sum(w * hh))
+ # browser()
+ stddev <- sqrt(S)#*(kappa/mn)
+ }
+ X <- if(length(object$weights)) object$x * sqrt(object$weights) else object$x
+ if(method == "XtWX") {
+ mn <- sum(wts*w)/sum(wts)
+ X <- X * sqrt(w/mn)
+ }
+ R <- qr(X)$qr
+ R <- R[1L:p, 1L:p, drop = FALSE]
+ R[lower.tri(R)] <- 0
+ rinv <- solve(R, rinv)
+ dimnames(rinv) <- list(cnames, cnames)
+ rowlen <- (rinv^2 %*% rep(1, p))^0.5
+ names(rowlen) <- cnames
+ if(correlation) {
+ correl <- rinv * array(1/rowlen, c(p, p))
+ correl <- correl %*% t(correl)
+ } else correl <- NULL
+ coef <- array(coef, c(p, 3L))
+ dimnames(coef) <- list(cnames, c("Value", "Std. Error", "t value"))
+ coef[, 2] <- rowlen %o% stddev
+ coef[, 3] <- coef[, 1]/coef[, 2]
+ object <- object[c("call", "na.action")]
+ object$residuals <- res
+ object$coefficients <- coef
+ object$sigma <- s
+ object$stddev <- stddev
+ object$df <- c(p, rdf, ptotal)
+ object$r.squared <- NA
+ object$cov.unscaled <- rinv %*% t(rinv)
+ object$correlation <- correl
+ object$terms <- NA
+ class(object) <- "summary.iwlsm"
+ object
+}
+##@print.summary.iwlsm Methods
+print.summary.iwlsm <-
+function(x, digits = max(3, .Options$digits - 3), ...)
+{
+ cat("\nCall: ")
+ dput(x$call, control=NULL)
+ resid <- x$residuals
+ df <- x$df
+ rdf <- df[2L]
+ cat(if(!is.null(x$weights) && diff(range(x$weights))) "Weighted ",
+ "Residuals:\n", sep="")
+ if(rdf > 5L) {
+ if(length(dim(resid)) == 2L) {
+ rq <- apply(t(resid), 1L, quantile)
+ dimnames(rq) <- list(c("Min", "1Q", "Median", "3Q", "Max"),
+ colnames(resid))
+ } else {
+ rq <- quantile(resid)
+ names(rq) <- c("Min", "1Q", "Median", "3Q", "Max")
+ }
+ print(rq, digits = digits, ...)
+ } else if(rdf > 0L) {
+ print(resid, digits = digits, ...)
+ }
+ if(nsingular <- df[3L] - df[1L])
+ cat("\nCoefficients: (", nsingular,
+ " not defined because of singularities)\n", sep = "")
+ else cat("\nCoefficients:\n")
+ print(format(round(x$coefficients, digits = digits)), quote = FALSE, ...)
+ cat("\nResidual standard error:", format(signif(x$stddev, digits)),
+ "on", rdf, "degrees of freedom\n")
+ if(nzchar(mess <- naprint(x$na.action))) cat(" (", mess, ")\n", sep="")
+ if(!is.null(correl <- x$correlation)) {
+ p <- dim(correl)[2L]
+ if(p > 1L) {
+ cat("\nCorrelation of Coefficients:\n")
+ ll <- lower.tri(correl)
+ correl[ll] <- format(round(correl[ll], digits))
+ correl[!ll] <- ""
+ print(correl[-1L, -p, drop = FALSE], quote = FALSE, digits = digits, ...)
+ }
+ }
+ invisible(x)
+}
+
+##@se.contrast.iwlsm Methods
+se.contrast.iwlsm <-
+ function(object, contrast.obj, coef = contr.helmert(ncol(contrast))[, 1],
+ data = NULL, ...)
+{
+ contrast.weight.aov <- function(object, contrast)
+ {
+ asgn <- object$assign[object$qr$pivot[1L:object$rank]]
+ uasgn <- unique(asgn)
+ nterms <- length(uasgn)
+ nmeffect <- c("(Intercept)",
+ attr(object$terms, "term.labels"))[1L + uasgn]
+ effects <- as.matrix(qr.qty(object$qr, contrast))
+ res <- matrix(0, nrow = nterms, ncol = ncol(effects),
+ dimnames = list(nmeffect, colnames(contrast)))
+ for(i in seq(nterms)) {
+ select <- (asgn == uasgn[i])
+ res[i,] <- colSums(effects[seq_along(asgn)[select], , drop = FALSE]^2)
+ }
+ res
+ }
+ if(is.null(data)) contrast.obj <- eval(contrast.obj)
+ else contrast.obj <- eval(substitute(contrast.obj), data, parent.frame())
+ if(!is.matrix(contrast.obj)) { # so a list
+ if(sum(coef) != 0)
+ stop("'coef' must define a contrast, i.e., sum to 0")
+ if(length(coef) != length(contrast.obj))
+ stop("'coef' must have same length as 'contrast.obj'")
+ contrast <-
+ sapply(contrast.obj, function(x)
+ {
+ if(!is.logical(x))
+ stop(gettextf("each element of %s must be logical",
+ substitute(contrasts.list)), domain = NA)
+ x/sum(x)
+ })
+ contrast <- contrast %*% coef
+ if(!any(contrast) || all(is.na(contrast)))
+ stop("the contrast defined is empty (has no TRUE elements)")
+ } else {
+ contrast <- contrast.obj
+ if(any(abs(colSums(contrast)) > 1e-8))
+ stop("columns of 'contrast.obj' must define a contrast (sum to zero)")
+ if(length(colnames(contrast)) == 0L)
+ colnames(contrast) <- paste("Contrast", seq(ncol(contrast)))
+ }
+ weights <- contrast.weight.aov(object, contrast)
+ object$stddev * if(!is.matrix(contrast.obj)) sqrt(sum(weights)) else sqrt(colSums(weights))
+}
+
+##@predict.iwlsm Methods
+predict.iwlsm <- function (object, newdata = NULL, scale = NULL, ...)
+{
+ ## problems with using predict.lm are the scale and
+ ## the QR decomp which has been done on down-weighted values.
+ object$qr <- qr(sqrt(object$weights) * object$x)
+ predict.lm(object, newdata = newdata, scale = object$s, ...)
+}
+
+##@vcov.iwlsm Methods
+vcov.iwlsm <- function (object, ...)
+{
+ so <- summary(object, corr = FALSE)
+ so$stddev^2 * so$cov.unscaled
+}
+##@psi.iwlsm iwlsm
+psi.iwlsm <- function(u, k, deriv=0, w, sj2)
+{
+ if (!deriv)
+ {
+ v <- sum(w * u^2) / (1 - sum(w * w))
+ v1 <- max(0, v - weighted.mean(sj2, w))
+ ww <- 1 /(v1 + sj2)
+ ww / sum(ww)
+ }
+ else ## dummy: I have removed the call with deriv = 1
+ {
+ rep(1, length(u))
+ }
+
+}
Modified: pkg/RSienaTest/R/phase1.r
===================================================================
--- pkg/RSienaTest/R/phase1.r 2010-03-15 14:56:31 UTC (rev 62)
+++ pkg/RSienaTest/R/phase1.r 2010-03-15 15:02:21 UTC (rev 63)
@@ -56,6 +56,8 @@
zsmall$Phase <- z$Phase
zsmall$nit <- z$nit
zsmall$FinDiff.method <- z$FinDiff.method
+ zsmall$int2 <- z$int2
+ zsmall$cl <- z$cl
xsmall<- NULL
zsmall$cconditional <- z$cconditional
zsmall$condvar <- z$condvar
@@ -97,8 +99,7 @@
## #######################################
## # do iteration
## #######################################
- z <- doPhase1it(z, x, cl=z$cl, int=z$int, zsmall=zsmall,
- xsmall=xsmall, ...)
+ z <- doPhase1it(z, x, zsmall=zsmall, xsmall=xsmall, ...)
## #######################################
## #
## #######################################
@@ -164,8 +165,9 @@
}
##@doPhase1it siena07 does 1 iteration in Phase 1
-doPhase1it<- function(z, x, cl, int, zsmall, xsmall, ...)
+doPhase1it<- function(z, x, zsmall, xsmall, ...)
{
+ int <- z$int
DisplayIteration(z)
if (int == 1)
{
@@ -181,9 +183,10 @@
}
else
{
- zz <- clusterCall(cl, usesim, zsmall, xsmall)
+ zz <- clusterCall(z$cl, usesim, zsmall, xsmall)
z$n <- z$n + z$int
z$phase1Its <- z$phase1Its + int
+ # browser()
}
if (int == 1)
{
@@ -207,7 +210,7 @@
}
if (z$FinDiff.method)
{
- z <- FiniteDifferences(z, x, fra + z$targets, z$cl, z$int, ...)
+ z <- FiniteDifferences(z, x, fra + z$targets, ...)
z$sdf[z$nit:(z$nit + (z$int - 1)), , ] <- z$sdf0
}
else if (x$maxlike)
@@ -273,6 +276,8 @@
zsmall$Deriv <- z$Deriv
zsmall$Phase <- z$Phase
zsmall$nit <- z$nit
+ zsmall$int2 <- z$int2
+ zsmall$cl <- z$cl
xsmall<- NULL
zsmall$cconditional <- z$cconditional
zsmall$condvar <- z$condvar
@@ -292,8 +297,7 @@
}
}
z$nit <- nit
- z <- doPhase1it(z, x, cl=z$cl, int=int, zsmall=zsmall,
- xsmall=xsmall, ...)
+ z <- doPhase1it(z, x, zsmall=zsmall, xsmall=xsmall, ...)
if (!z$OK || UserInterruptFlag() || UserRestartFlag())
{
return(z)
@@ -477,14 +481,16 @@
}
##@FiniteDifferences siena07 Does the extra iterations for finite differences
-FiniteDifferences <- function(z, x, fra, cl, int=1, ...)
+FiniteDifferences <- function(z, x, fra, ...)
{
+ int <- z$int
fras <- array(0, dim = c(int, z$pp, z$pp))
xsmall<- NULL
##browser()
for (i in 1 : z$pp)
{
- zdummy <- z[c('theta', 'Deriv', 'cconditional', 'FinDiff.method')]
+ zdummy <- z[c('theta', 'Deriv', 'cconditional', 'FinDiff.method',
+ 'int2', 'cl')]
if (!z$fixed[i])
{
zdummy$theta[i] <- z$theta[i] + z$epsilon[i]
@@ -502,8 +508,9 @@
}
else
{
- zz <- clusterCall(cl, usesim, zdummy, xsmall,
+ zz <- clusterCall(z$cl, usesim, zdummy, xsmall,
INIT=FALSE, fromFiniteDiff=TRUE)
+ #browser()
}
if (int == 1)
{
Modified: pkg/RSienaTest/R/phase2.r
===================================================================
--- pkg/RSienaTest/R/phase2.r 2010-03-15 14:56:31 UTC (rev 62)
+++ pkg/RSienaTest/R/phase2.r 2010-03-15 15:02:21 UTC (rev 63)
@@ -26,8 +26,11 @@
phase2.1<- function(z, x, ...)
{
#initialise phase2
- z$phase2fras <- array(0, dim=c(4, z$pp, 1000))
- z$rejectprops <- array(0, dim=c(4, 4, 1000))
+ if (x$maxlike)
+ {
+ z$phase2fras <- array(0, dim=c(4, z$pp, 1000))
+ z$rejectprops <- array(0, dim=c(4, 4, 1000))
+ }
int <- 1
f <- FRANstore()
z$Phase <- 2
@@ -191,7 +194,9 @@
zsmall$theta <- z$theta
zsmall$Deriv <- z$Deriv
zsmall$Phase <- z$Phase
+ zsmall$int2 <- z$int2
zsmall$FinDiff.method <- z$FinDiff.method
+ zsmall$cl <- z$cl
xsmall <- NULL
zsmall$cconditional <- z$cconditional
zsmall$condvar <- z$condvar
@@ -269,9 +274,9 @@
break
}
}
- z$phase2fras[subphase, ,z$nit] <- fra
if (x$maxlike)
{
+ z$phase2fras[subphase, ,z$nit] <- fra
z$rejectprops[subphase, , z$nit] <- zz$rejectprop
}
if (z$nit %% 2 == 1)
Modified: pkg/RSienaTest/R/phase3.r
===================================================================
--- pkg/RSienaTest/R/phase3.r 2010-03-15 14:56:31 UTC (rev 62)
+++ pkg/RSienaTest/R/phase3.r 2010-03-15 15:02:21 UTC (rev 63)
@@ -66,6 +66,8 @@
zsmall$Deriv <- z$Deriv
zsmall$Phase <- z$Phase
zsmall$nit <- z$nit
+ zsmall$int2 <- z$int2
+ zsmall$cl <- z$cl
xsmall<- NULL
zsmall$cconditional <- z$cconditional
zsmall$condvar <- z$condvar
@@ -124,7 +126,8 @@
}
}
if (nit <= 5 || nit == 10 || (int==1 && nit %% z$writefreq == 0 ) ||
- (int > 1 && nit %in% nits[seq(z$writefreq + 1, x$n3 %/% int,
+ (int > 1 && (z$writefreq - 1) < z$n3%/%int &&
+ nit %in% nits[seq(z$writefreq + 1, x$n3 %/% int,
z$writefreq)]))
{
if (!is.batch())
@@ -156,8 +159,7 @@
}
}
###############################
- z <- doPhase3it(z, x, nit, cl=z$cl, int=int, zsmall=zsmall,
- xsmall=xsmall, ...)
+ z <- doPhase3it(z, x, nit, zsmall=zsmall, xsmall=xsmall, ...)
##############################
## browser()
if (!is.batch())
@@ -208,8 +210,9 @@
}
##@doPhase3it siena07 Does one iteration in phase 3
-doPhase3it<- function(z, x, nit, cl, int, zsmall, xsmall, ...)
+doPhase3it<- function(z, x, nit, zsmall, xsmall, ...)
{
+ int <- z$int
if (int == 1)
{
zz <- x$FRAN(zsmall, xsmall)
@@ -224,7 +227,7 @@
else
{
##zz <- clusterCall(cl, simstats0c, zsmall, xsmall)
- zz <- clusterCall(cl, usesim, zsmall, xsmall)
+ zz <- clusterCall(z$cl, usesim, zsmall, xsmall)
z$n <- z$n + z$int
# browser()
}
@@ -250,16 +253,16 @@
if ((!x$maxlike) && z$cconditional)
{
if (int==1)
- z$ntim[nit,]<- zz$ntim0
+ z$ntim[nit,] <- zz$ntim0
else
{
for (i in 1:int)
- z$ntim[nit+(i-1),]<- zz[[i]]$ntim0
+ z$ntim[nit+(i-1),] <- zz[[i]]$ntim0
}
}
if (z$FinDiff.method)
{
- z <- FiniteDifferences(z, x, fra + z$targets, cl, int, ...)
+ z <- FiniteDifferences(z, x, fra + z$targets, ...)
z$sdf[nit:(nit + (int - 1)), , ] <- z$sdf0
}
else if (x$maxlike) ## as far as I can see
Modified: pkg/RSienaTest/R/print01Report.r
===================================================================
--- pkg/RSienaTest/R/print01Report.r 2010-03-15 14:56:31 UTC (rev 62)
+++ pkg/RSienaTest/R/print01Report.r 2010-03-15 15:02:21 UTC (rev 63)
@@ -292,8 +292,8 @@
Report(c("Total number of missing data: ",
sum(missrow),
", corresponding to a fraction of ",
- round(sum(missrow)/atts$netdims[1] /
- (atts$netdims[1] - 1), 3),
+ format(round(sum(missrow)/atts$netdims[1] /
+ (atts$netdims[1] - 1), 3), nsmall=3),
".\n"), sep="", outf)
if (k > 1)
Report(c("In reported in- and outdegrees,",
@@ -649,13 +649,31 @@
use <- ! covars %in% names(x$dycCovars) ## need an attributes to say
nCovars <- length(x$dyvCovars[use])
Heading(2, outf, "Reading exogenous dyadic covariates.")
- Report(c("Note that no missing values are considered yet for",
- "changing dyadic covariates.\n"), outf)
+ ## Report(c("Note that no missing values are considered yet for",
+ ## "changing dyadic covariates.\n"), outf)
for (i in seq(along=covars))
{
Report(c("Exogenous dyadic covariate named ", covars[i], '.\n'),
sep="", outf)
}
+ Report("Number of tie variables with missing data per period:\n", outf)
+ Report(c(" period ", format(1:(x$observations - 1) +
+ periodFromStart, width=9),
+ " overall\n"), sep="", outf)
+ for (i in seq(along=covars))
+ {
+ if (use[i])
+ {
+ thiscovar <- x$dyvCovars[[i]] ## array
+ missvals <- colSums(is.na(thiscovar), dims=2)
+ Report(c(format(covars[i], width=10),
+ format(missvals, width=8),
+ format(sum(missvals), width=9), " (",
+ format(round(sum(missvals)/nrow(thiscovar)/
+ ncol(thiscovar), 1), nsmall=1,
+ width=3), '%)\n'), outf)
+ }
+ }
Report("\nMeans of covariates:\n", outf)
for (i in seq(along=covars))
{
@@ -747,7 +765,7 @@
tt <- getInternals()
return(tt)
}
- Report(open=TRUE, type="w", projname=modelname)
+ Report(openfiles=TRUE, type="w", projname=modelname)
Report(" ************************\n", outf)
Report(c(" ", modelname, ".out\n"),
sep='', outf)
Modified: pkg/RSienaTest/R/robmon.r
===================================================================
--- pkg/RSienaTest/R/robmon.r 2010-03-15 14:56:31 UTC (rev 62)
+++ pkg/RSienaTest/R/robmon.r 2010-03-15 15:02:21 UTC (rev 63)
@@ -12,7 +12,8 @@
## z: model fitting object
## returns updated z
##@robmon siena07 Controls MOM process
-robmon <- function(z, x, useCluster, nbrNodes, initC, clusterString, ...)
+robmon <- function(z, x, useCluster, nbrNodes, initC, clusterString,
+ clusterIter, ...)
{
z$FinDiff.method<- x$FinDiff.method
z$n <- 0
@@ -29,6 +30,10 @@
#######################################################
##do initial setup call of FRAN
#######################################################
+ if (!is.function(x$FRAN))
+ {
+ x$FRAN <- getFromNamespace(x$FRANname, pos=grep("RSiena", search())[1])
+ }
z <- x$FRAN(z, x, INIT=TRUE, ...)
##
##if conditional, FRAN changes z$theta etc
@@ -39,10 +44,14 @@
{
stop("Multiple processors only for simstats0c at present")
}
+ if (!clusterIter && nbrNodes >= z$f$observations)
+ {
+ stop("Not enough observations to use the nodes")
+ }
cl <- makeCluster(clusterString, type = "SOCK",
- outfile = 'cluster.out')
- clusterSetupRNG(cl, seed = rep(1, 6))
+ outfile = "cluster.out")
clusterCall(cl, library, pkgname, character.only = TRUE)
+ clusterSetupRNG(cl, seed = as.integer(runif(6, max=.Machine$integer.max)))
clusterCall(cl, storeinFRANstore, FRANstore())
if (initC)
{
@@ -50,11 +59,10 @@
INIT=FALSE, initC = initC)
}
z$cl <- cl
- z$int <- nbrNodes
}
else
{
- z$int <- 1
+ z$cl <- NULL
}
z$newFixed <- rep(FALSE, z$pp)
z$AllNowFixed <- FALSE
Modified: pkg/RSienaTest/R/siena01.r
===================================================================
--- pkg/RSienaTest/R/siena01.r 2010-03-15 14:56:31 UTC (rev 62)
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/rsiena -r 63
More information about the Rsiena-commits
mailing list