[Rsiena-commits] r69 - in pkg/RSiena: . R data inst/examples man src src/data src/model src/model/effects src/model/filters src/model/ml src/model/variables src/utils
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Mar 24 19:12:18 CET 2010
Author: ripleyrm
Date: 2010-03-24 19:12:18 +0100 (Wed, 24 Mar 2010)
New Revision: 69
Added:
pkg/RSiena/R/iwlsm.R
pkg/RSiena/R/siena08.r
pkg/RSiena/man/iwlsm.Rd
pkg/RSiena/man/print.sienaMeta.Rd
pkg/RSiena/man/siena08.Rd
pkg/RSiena/man/summary.iwlsm.Rd
pkg/RSiena/src/model/ml/MLSimulation.cpp
pkg/RSiena/src/model/ml/MLSimulation.h
pkg/RSiena/src/model/ml/Option.cpp
pkg/RSiena/src/model/ml/Option.h
Modified:
pkg/RSiena/NAMESPACE
pkg/RSiena/R/RSienaRDocumentation.r
pkg/RSiena/R/getTargets.r
pkg/RSiena/R/globals.r
pkg/RSiena/R/phase1.r
pkg/RSiena/R/phase2.r
pkg/RSiena/R/phase3.r
pkg/RSiena/R/print01Report.r
pkg/RSiena/R/printInitialDescription.r
pkg/RSiena/R/robmon.r
pkg/RSiena/R/siena01.r
pkg/RSiena/R/siena07.r
pkg/RSiena/R/sienaDataCreate.r
pkg/RSiena/R/sienaDataCreateFromSession.r
pkg/RSiena/R/sienaModelCreate.r
pkg/RSiena/R/sienaprint.r
pkg/RSiena/R/simstatsc.r
pkg/RSiena/R/zzz.R
pkg/RSiena/changeLog
pkg/RSiena/data/allEffects.csv
pkg/RSiena/inst/examples/baerveldt3.csv
pkg/RSiena/inst/examples/baerveldt4.csv
pkg/RSiena/man/RSiena-package.Rd
pkg/RSiena/man/siena07.Rd
pkg/RSiena/man/sienaModelCreate.Rd
pkg/RSiena/src/data/BehaviorLongitudinalData.cpp
pkg/RSiena/src/data/BehaviorLongitudinalData.h
pkg/RSiena/src/data/Data.cpp
pkg/RSiena/src/data/LongitudinalData.cpp
pkg/RSiena/src/data/LongitudinalData.h
pkg/RSiena/src/data/NetworkLongitudinalData.cpp
pkg/RSiena/src/data/NetworkLongitudinalData.h
pkg/RSiena/src/data/OneModeNetworkLongitudinalData.cpp
pkg/RSiena/src/data/OneModeNetworkLongitudinalData.h
pkg/RSiena/src/model/EpochSimulation.cpp
pkg/RSiena/src/model/EpochSimulation.h
pkg/RSiena/src/model/effects/EffectFactory.cpp
pkg/RSiena/src/model/filters/AtLeastOneFilter.cpp
pkg/RSiena/src/model/filters/AtLeastOneFilter.h
pkg/RSiena/src/model/filters/DisjointFilter.cpp
pkg/RSiena/src/model/filters/DisjointFilter.h
pkg/RSiena/src/model/filters/HigherFilter.cpp
pkg/RSiena/src/model/filters/HigherFilter.h
pkg/RSiena/src/model/filters/LowerFilter.cpp
pkg/RSiena/src/model/filters/LowerFilter.h
pkg/RSiena/src/model/filters/PermittedChangeFilter.h
pkg/RSiena/src/model/ml/BehaviorChange.cpp
pkg/RSiena/src/model/ml/BehaviorChange.h
pkg/RSiena/src/model/ml/Chain.cpp
pkg/RSiena/src/model/ml/Chain.h
pkg/RSiena/src/model/ml/MiniStep.cpp
pkg/RSiena/src/model/ml/MiniStep.h
pkg/RSiena/src/model/ml/NetworkChange.cpp
pkg/RSiena/src/model/ml/NetworkChange.h
pkg/RSiena/src/model/variables/BehaviorVariable.cpp
pkg/RSiena/src/model/variables/BehaviorVariable.h
pkg/RSiena/src/model/variables/DependentVariable.cpp
pkg/RSiena/src/model/variables/DependentVariable.h
pkg/RSiena/src/model/variables/NetworkVariable.cpp
pkg/RSiena/src/model/variables/NetworkVariable.h
pkg/RSiena/src/siena07.cpp
pkg/RSiena/src/utils/Random.cpp
pkg/RSiena/src/utils/Random.h
Log:
Siena08 meta analysis, multiple processors by wave, fixes
Modified: pkg/RSiena/NAMESPACE
===================================================================
--- pkg/RSiena/NAMESPACE 2010-03-21 15:37:21 UTC (rev 68)
+++ pkg/RSiena/NAMESPACE 2010-03-24 18:12:18 UTC (rev 69)
@@ -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,
- installGui)#, sienaTimeTest)
+ effectsDocumentation, sienaDataConstraint, #maxlikefn,
+ installGui, siena08, iwlsm)#, sienaTimeTest)
import(Matrix)
import(xtable)
@@ -18,6 +18,22 @@
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, coCovar)
+S3method(addAttributes, varCovar)
+S3method(addAttributes, coDyadCovar)
+S3method(addAttributes, varDyadCovar)
#S3method(print, sienaTimeTest)
#S3method(summary, sienaTimeTest)
#S3method(print, summary.sienaTimeTest)
Modified: pkg/RSiena/R/RSienaRDocumentation.r
===================================================================
--- pkg/RSiena/R/RSienaRDocumentation.r 2010-03-21 15:37:21 UTC (rev 68)
+++ pkg/RSiena/R/RSienaRDocumentation.r 2010-03-24 18:12:18 UTC (rev 69)
@@ -80,7 +80,7 @@
## get the calls (global)
codet <- lapply(comms3[,2], function(x)
{
- x <- try(getFromNamespace(x, "RSiena"), silent=TRUE)
+ x <- try(getFromNamespace(x, pkgname), silent=TRUE)
if (is.function(x))
{
tmp1 <- findGlobals(x, merge=FALSE)[[1]]
@@ -108,7 +108,7 @@
{
yy <- y[x]
zz <- z[[x]]
- yy <- getFromNamespace(yy, "RSiena")
+ yy <- getFromNamespace(yy, pkgname)
targs <- formals(yy)
n <- length(targs)
myargs <- targs
Modified: pkg/RSiena/R/getTargets.r
===================================================================
--- pkg/RSiena/R/getTargets.r 2010-03-21 15:37:21 UTC (rev 68)
+++ pkg/RSiena/R/getTargets.r 2010-03-24 18:12:18 UTC (rev 69)
@@ -3,32 +3,30 @@
{
f <- unpackData(data)
effects <- effects[effects$include,]
- ##dyn.load(dllpath)
##
- ## dyn.load('d:/sienasvn/siena/src/RSiena.dll')
- pData <- .Call('setupData', PACKAGE="RSiena",
+ pData <- .Call('setupData', PACKAGE=pkgname,
list(as.integer(f$observations)),
list(f$nodeSets))
## register a finalizer
ans <- reg.finalizer(pData, clearData, onexit = FALSE)
- ans<- .Call('OneMode', PACKAGE="RSiena",
+ ans<- .Call('OneMode', PACKAGE=pkgname,
pData, list(f$nets))
- ans<- .Call('Behavior', PACKAGE="RSiena", pData,
+ ans<- .Call('Behavior', PACKAGE=pkgname, pData,
list(f$behavs))
- ans<-.Call('ConstantCovariates', PACKAGE="RSiena",
+ ans<-.Call('ConstantCovariates', PACKAGE=pkgname,
pData, list(f$cCovars))
- ans<-.Call('ChangingCovariates',PACKAGE="RSiena",
+ ans<-.Call('ChangingCovariates',PACKAGE=pkgname,
pData,list(f$vCovars))
- ans<-.Call('DyadicCovariates',PACKAGE="RSiena",
+ ans<-.Call('DyadicCovariates',PACKAGE=pkgname,
pData,list(f$dycCovars))
- ans<-.Call('ChangingDyadicCovariates',PACKAGE="RSiena",
+ ans<-.Call('ChangingDyadicCovariates',PACKAGE=pkgname,
pData, list(f$dyvCovars))
storage.mode(effects$parm) <- 'integer'
storage.mode(effects$group) <- 'integer'
storage.mode(effects$period) <- 'integer'
effects$effectPtr <- NA
myeffects <- split(effects, effects$name)
- ans<- .Call('effects', PACKAGE="RSiena",
+ ans<- .Call('effects', PACKAGE=pkgname,
pData, myeffects)
pModel <- ans[[1]][[1]]
for (i in 1:length(ans[[2]])) ## ans[[2]] is a list of lists of
@@ -38,7 +36,7 @@
effectPtr <- ans[[2]][[i]]
myeffects[[i]]$effectPtr <- effectPtr
}
- ans <- .Call('getTargets', PACKAGE="RSiena",
+ ans <- .Call('getTargets', PACKAGE=pkgname,
pData, pModel, myeffects)
ans
}
Modified: pkg/RSiena/R/globals.r
===================================================================
--- pkg/RSiena/R/globals.r 2010-03-21 15:37:21 UTC (rev 68)
+++ pkg/RSiena/R/globals.r 2010-03-24 18:12:18 UTC (rev 69)
@@ -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/RSiena/R/iwlsm.R
===================================================================
--- pkg/RSiena/R/iwlsm.R (rev 0)
+++ pkg/RSiena/R/iwlsm.R 2010-03-24 18:12:18 UTC (rev 69)
@@ -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))
+ }
+
+}
Property changes on: pkg/RSiena/R/iwlsm.R
___________________________________________________________________
Name: svn:eol-style
+ native
Modified: pkg/RSiena/R/phase1.r
===================================================================
--- pkg/RSiena/R/phase1.r 2010-03-21 15:37:21 UTC (rev 68)
+++ pkg/RSiena/R/phase1.r 2010-03-24 18:12:18 UTC (rev 69)
@@ -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/RSiena/R/phase2.r
===================================================================
--- pkg/RSiena/R/phase2.r 2010-03-21 15:37:21 UTC (rev 68)
+++ pkg/RSiena/R/phase2.r 2010-03-24 18:12:18 UTC (rev 69)
@@ -26,6 +26,11 @@
phase2.1<- function(z, x, ...)
{
#initialise phase2
+ 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
@@ -81,6 +86,7 @@
z$ctime <- proc.time()[3]
z$time1 <- proc.time()[3]
z$thav <- z$theta
+ z$thavn <- 1
## cat(z$thav, z$theta, '\n')
z$prod0 <- rep(0, z$pp)
z$prod1 <- rep(0, z$pp)
@@ -149,7 +155,7 @@
' = ', format(subphaseTime, nsmall=4, digits=4),
'\n', sep=''), lf)
}
- z$theta <- z$thav / (z$nit + 1)
+ z$theta <- z$thav / z$thavn #(z$nit + 1)
DisplayThetaAutocor(z)
## cat('it',z$nit,'\n')
##recalculate autocor using -1 instead of -2 as error
@@ -188,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
@@ -243,9 +251,11 @@
}
}
}
+ zsmall$nit <- z$nit
if (z$int == 1)
{
zz <- x$FRAN(zsmall, xsmall)
+ ## browser()
fra <- colSums(zz$fra) - z$targets
if (!zz$OK)
{
@@ -264,6 +274,11 @@
break
}
}
+ if (x$maxlike)
+ {
+ z$phase2fras[subphase, ,z$nit] <- fra
+ z$rejectprops[subphase, , z$nit] <- zz$rejectprop
+ }
if (z$nit %% 2 == 1)
{
prev.fra <- fra
@@ -281,9 +296,7 @@
DisplayThetaAutocor(z)
}
}
- ## limit change. not sure what to do here sd is not set up
- ## unless finite differences are used or ML and
- ## ML is specifically excluded here. Reporting is delayed to
+ ## limit change. Reporting is delayed to
## end of phase.
## browser()
if (x$diag)## !maxlike at present
@@ -300,7 +313,6 @@
}
else
{
- browser()
fchange <- as.vector(z$gain * fra %*% z$dinv)
}
## browser()
@@ -311,6 +323,12 @@
zsmall$theta <- zsmall$theta - fchange
z$theta <- zsmall$theta
z$thav <- z$thav + zsmall$theta
+ z$thavn <- z$thavn + 1
+ ## if (x$maxlike && !is.null(x$moreUpdates) && x$moreUpdates > 0)
+ ## {
+ ## z <- doMoreUpdates(z, x, x$moreUpdates * subphase)
+ ## zsmall$theta <- z$theta
+ ## }
##check for user interrupt
## browser()
CheckBreaks()
Modified: pkg/RSiena/R/phase3.r
===================================================================
--- pkg/RSiena/R/phase3.r 2010-03-21 15:37:21 UTC (rev 68)
+++ pkg/RSiena/R/phase3.r 2010-03-24 18:12:18 UTC (rev 69)
@@ -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,7 @@
}
}
if (nit <= 5 || nit == 10 || (int==1 && nit %% z$writefreq == 0 ) ||
- (int > 1 && (z$writefreq - 1) < z$n3%/%int &&
+ (int > 1 && (z$writefreq + 1) < z$n3%/%int &&
nit %in% nits[seq(z$writefreq + 1, x$n3 %/% int,
z$writefreq)]))
{
@@ -157,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())
@@ -209,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)
@@ -225,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()
}
@@ -248,19 +250,19 @@
z$sims[[nit + (i - 1)]] <- zz[[i]]$sims
}
}
- if (z$cconditional)
+ 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
@@ -309,7 +311,7 @@
Report(c('Total of', z$n,'iterations.\n'), outf)
Report(c('Parameter estimates based on', z$n - z$Phase3nits,
'iterations,\n'), outf)
- if (z$cconditional)
+ if (!x$maxlike && z$cconditional)
Report(c('basic rate parameter',
c('', 's')[as.integer(z$observations > 2) + 1],
' as well as \n'), outf)
@@ -319,10 +321,10 @@
Report(c('Averages, standard deviations, ',
'and t-ratios for deviations from targets:\n'), sep='', outf)
# Report(c(date(),'\n'),bof)
- if (z$cconditional)
+ if (x$maxlike)
+ Report('\nMaximum Likelihood estimation.', bof)
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/rsiena -r 69
More information about the Rsiena-commits
mailing list