[Rsiena-commits] r23 - in pkg/RSiena: R src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sun Nov 8 18:11:12 CET 2009
Author: ripleyrm
Date: 2009-11-08 18:11:12 +0100 (Sun, 08 Nov 2009)
New Revision: 23
Modified:
pkg/RSiena/R/Sienatest.r
pkg/RSiena/R/effects.r
pkg/RSiena/R/phase1.r
pkg/RSiena/R/phase3.r
pkg/RSiena/R/siena07.r
pkg/RSiena/R/sienaprint.r
pkg/RSiena/R/simstatsc.r
pkg/RSiena/src/siena07.cpp
Log:
bug fixes, internal changes to score test and derivative calculation, display effect number on siena07 screen.
Modified: pkg/RSiena/R/Sienatest.r
===================================================================
--- pkg/RSiena/R/Sienatest.r 2009-11-08 17:08:39 UTC (rev 22)
+++ pkg/RSiena/R/Sienatest.r 2009-11-08 17:11:12 UTC (rev 23)
@@ -67,9 +67,9 @@
tmp
}
##@TestOutput siena07 Print report
-TestOutput <- function(z,x)
+TestOutput <- function(z, x)
{
- testn<- sum(z$test)
+ testn <- sum(z$test)
# browser()
if (testn)
{
@@ -95,7 +95,7 @@
Report(' \n',outf)
if (testn > 1)
Report('Joint test:\n-----------\n',outf)
- Report(c(' c = ',sprintf("%8.4f",z$testresOverall),
+ Report(c(' c = ',sprintf("%8.4f", z$testresOverall),
' d.f. = ',j,' p-value '),sep='',outf)
pvalue <- 1-pchisq(z$testresOverall,j)
if (pvalue < 0.0001)
@@ -141,45 +141,50 @@
}
}
##@ScoreTest siena07 Do score tests
-ScoreTest<- function(z,x)
+ScoreTest<- function(pp, dfra, msf, fra, test, maxlike)
{
- z$testresult<- rep(NA,z$pp) ##for chisq per parameter
- z$testresulto <- rep(NA,z$pp) ##for one-sided tests per parameter
+ testresult<- rep(NA, pp) ##for chisq per parameter
+ testresulto <- rep(NA, pp) ##for one-sided tests per parameter
##first the general one
- ans<-EvaluateTestStatistic(x,z$test,z$dfra,z$msf,z$fra)
- z$testresOverall<- ans$cvalue
- if (sum(z$test)==1)
- z$testresulto[1]<- ans$oneSided
+ ans <- EvaluateTestStatistic(maxlike, test, dfra, msf, fra)
+ testresOverall <- ans$cvalue
+ covMatrix <- ans$covMatrix
+ if (sum(test) == 1)
+ {
+ testresulto[1] <- ans$oneSided
+ }
else
{
## single df tests
- use<- !z$test
- k<- 0
- for (i in 1:z$pp)
+ use <- !test
+ k <- 0
+ for (i in 1:pp)
{
- if (z$test[i])
+ if (test[i])
{
- k<- k+1
- use[i]<- TRUE
- ans<-EvaluateTestStatistic(x,z$test[use],z$dfra[use,use],
- z$msf[use,use],z$fra[use])
- z$testresult[k]<- ans$cvalue
- z$testresulto[k]<- ans$oneSided
- use[i]<- FALSE
+ k <- k+1
+ use[i] <- TRUE
+ ans <- EvaluateTestStatistic(maxlike, test[use], dfra[use, use],
+ msf[use, use], fra[use])
+ testresult[k] <- ans$cvalue
+ testresulto[k] <- ans$oneSided
+ use[i] <- FALSE
}
}
}
##onestep estimator
- if (x$maxlike)
- dfra2<- z$dfra+ z$msf
+ if (maxlike)
+ dfra2 <- dfra + msf
else
- dfra2<- z$dfra
- dinv2<- solve(dfra2)
- z$oneStep<- -dinv2%*%z$fra
- z
+ dfra2 <- dfra
+ dinv2 <- solve(dfra2)
+ oneStep<- -dinv2 %*% fra
+ list(testresult=testresult, testresulto=testresulto,
+ testresOverall=testresOverall, covMatrix=covMatrix,
+ oneStep=oneStep, dinv2= dinv2, dfra2=dfra2)
}
##@EvaluateTestStatistic siena07 Calculate score test statistics
-EvaluateTestStatistic<- function(x,test,dfra,msf,fra)
+EvaluateTestStatistic<- function(maxlike, test, dfra, msf, fra)
{
##uses local arrays set up in the calling procedure
d11 <- dfra[!test,!test,drop=FALSE]
@@ -194,7 +199,7 @@
z2 <- fra[test]
id11 <- solve(d11)
rg<- d21%*%id11
- if (!x$maxlike)
+ if (!maxlike)
{
##orthogonalise deviation vector
ov<- z2-rg%*%z1
@@ -219,10 +224,10 @@
oneSided <- ov * sqrt(vav)
else
oneSided <- 0
- if (!x$maxlike) oneSided<- - oneSided
+ if (maxlike) oneSided<- - oneSided
## change the sign for intuition for users
}
else
oneSided <- 0
- list(cvalue=cvalue,oneSided=oneSided)
+ list(cvalue=cvalue, oneSided=oneSided, covMatrix=v9)
}
Modified: pkg/RSiena/R/effects.r
===================================================================
--- pkg/RSiena/R/effects.r 2009-11-08 17:08:39 UTC (rev 22)
+++ pkg/RSiena/R/effects.r 2009-11-08 17:11:12 UTC (rev 23)
@@ -111,7 +111,7 @@
}
for (j in seq(along = xx$dycCovars))
{
- if (attr(x$dycCovars[[j]], 'nodeSet')[1] == nodeSet)
+ if (attr(xx$dycCovars[[j]], 'nodeSet')[1] == nodeSet)
{
objEffects <- rbind(objEffects,
createEffects("dyadObjective",
@@ -120,7 +120,7 @@
}
for (j in seq(along = xx$dyvCovars))
{
- if (attr(x$dvvCovars[[j]], 'nodeSet')[1] == nodeSet)
+ if (attr(xx$dvvCovars[[j]], 'nodeSet')[1] == nodeSet)
{
objEffects <- rbind(objEffects,
createEffects("dyadObjective",
@@ -129,7 +129,7 @@
}
for (j in seq(along = xx$cCovars))
{
- if (attr(x$cCovars[[j]], 'nodeSet') == nodeSet)
+ if (attr(xx$cCovars[[j]], 'nodeSet') == nodeSet)
{
tmp <- covarOneModeEff(names(xx$cCovars)[j],
attr(xx$cCovars[[j]], 'poszvar'),
@@ -152,9 +152,9 @@
rateEffects <- rbind(rateEffects, tmp$rateEff)
}
}
- for (j in seq(along=x$vCovars))
+ for (j in seq(along=xx$vCovars))
{
- if (attr(x$cCovars[[j]], 'nodeSet') == nodeSet)
+ if (attr(xx$vCovars[[j]], 'nodeSet') == nodeSet)
{
tmp <- covarOneModeEff(names(xx$vCovars)[j],
attr(xx$vCovars[[j]], 'poszvar'),
@@ -491,7 +491,7 @@
}
for (j in seq(along=xx$vCovars))
{
- covNodeset <- match(attr(xx$cCovars[[j]], "nodeSet"),
+ covNodeset <- match(attr(xx$vCovars[[j]], "nodeSet"),
nodeSets)
if (covNodeset > 0)
{
Modified: pkg/RSiena/R/phase1.r
===================================================================
--- pkg/RSiena/R/phase1.r 2009-11-08 17:08:39 UTC (rev 22)
+++ pkg/RSiena/R/phase1.r 2009-11-08 17:11:12 UTC (rev 23)
@@ -417,22 +417,9 @@
##note that warnings is never set as this piece of code is not executed
##if force_fin_diff_phase_1 is true so warning is zero on entry
##browser()
- dfra<- matrix(0, nrow = z$pp, ncol = z$pp)
- # browser()
- for (i in 1 : (f$observations-1))
- {
- tmp <- sapply(1 : z$n1,function(j, y, s)
- outer(y[j, ], s[j, ]),
- y = matrix(z$sf2[, i, ], ncol=z$pp),
- s = matrix(z$ssc[, i, ], ncol=z$pp))
- tmp <- array(tmp, dim = c(z$pp, z$pp, z$n1))
- dfra <- dfra + apply(tmp, c(1, 2), sum)
- }
- dfra <- dfra / z$n1
- tmp <- matrix(sapply(1 : (f$observations - 1), function(i)
- outer(colMeans(z$sf2)[i,],
- colMeans(z$ssc)[i,])), ncol=f$observations-1)
- dfra <- dfra - matrix(rowSums(tmp), nrow = z$pp)
+
+ dfra <-derivativeFromScoresAndDeviations(z$ssc, z$sf2)
+
z$jacobianwarn1 <- rep(FALSE, z$pp)
if (any(diag(dfra) <= 0))
{
@@ -545,3 +532,23 @@
# browser()
z
}
+##@derivativeFromScoresAndDeviations siena07 create dfra from scores and deviations
+derivativeFromScoresAndDeviations <- function(scores, deviations)
+{
+ nIterations <- dim(scores)[1]
+ nWaves <- dim(scores)[2]
+ nParameters <- dim(scores)[3]
+ dfra <- rowSums(sapply(1:nWaves, function(i)
+ sapply(1:nIterations, function(j, y, s)
+ outer(y[j, ], s[j, ]),
+ y=matrix(deviations[, i, ], ncol=nParameters),
+ s=matrix(scores[, i, ], ncol=nParameters))))
+ dim(dfra)<- c(nParameters, nParameters, nIterations)
+ dfra<- apply(dfra, c(1, 2), sum)
+ dfra<- dfra / nIterations
+ tmp <- matrix(sapply(1 : nWaves, function(i)
+ outer(colMeans(deviations)[i,],
+ colMeans(scores)[i,])), ncol=nWaves)
+
+ dfra - matrix(rowSums(tmp), nrow=nParameters)
+}
Modified: pkg/RSiena/R/phase3.r
===================================================================
--- pkg/RSiena/R/phase3.r 2009-11-08 17:08:39 UTC (rev 22)
+++ pkg/RSiena/R/phase3.r 2009-11-08 17:11:12 UTC (rev 23)
@@ -443,20 +443,7 @@
}
else
{
- dfra <-rowSums(sapply(1:(f$observations-1), function(i,y,s)
- sapply(1:z$Phase3nits, function(j,y,s)
- outer(y[j,], s[j,]),
- y=matrix(z$sf2[,i,], ncol=z$pp),
- s=matrix(z$ssc[,i,], ncol=z$pp))))
- dim(dfra)<- c(z$pp, z$pp, z$Phase3nits)
- dfra<- apply(dfra,c(1,2),sum)
- dfra<- dfra/z$Phase3nits
- tmp <- matrix(sapply(1 : (f$observations - 1), function(i)
- outer(colMeans(z$sf2)[i,],
- colMeans(z$ssc)[i,])), ncol=f$observations-1)
-
- dfra<- dfra - matrix(rowSums(tmp), nrow=z$pp)
-
+ dfra <- derivativeFromScoresAndDeviations(z$ssc, z$sf2)
if (any(diag(dfra) < 0))
{
sub <- which(diag(dfra) < 0)
Modified: pkg/RSiena/R/siena07.r
===================================================================
--- pkg/RSiena/R/siena07.r 2009-11-08 17:08:39 UTC (rev 22)
+++ pkg/RSiena/R/siena07.r 2009-11-08 17:11:12 UTC (rev 23)
@@ -145,45 +145,7 @@
Report(sprintf("Current random number seed is %d.\n", seed), outf)
}
}
-##@WriteOutTheta siena07 Progress reporting
-WriteOutTheta <- function(z)
-{
- if (!is.batch())
- {
- tkdelete(z$tkvars$current, "1.0", "end")
- tmp <- paste(c("", rep("\n", z$pp - 1)),
- format(round(z$theta,4), width=12, sep=""),
- collapse="")
- tkinsert(z$tkvars$current, "1.0", tmp)
- }
- else
- {
- Report(c("theta:", format(z$theta, digits=3), "\n"))
- }
- Report("Current parameter values:\n", cf)
- Report(format(z$theta), cf, fill=80)
-}
-##@DisplayTheta siena07 Progress reporting
-DisplayTheta<- function(z)
-{
- if ((z$Phase == 2 || z$nit == 1 ) && (z$nit <= 30))
- {
- if (!is.batch())
- {
- tkdelete(z$tkvars$current, "1.0", "end")
- tmp<- paste(c("", rep("\n", z$pp - 1)),
- format(z$theta, width=12, sep=""),
- collapse="")
- tkinsert(z$tkvars$current, "1.0", tmp)
- }
- else
- {
- Report(c("theta:", format(z$theta, digits=3), "\n"))
- }
- }
-
-}
##@AnnouncePhase siena07 Progress reporting
AnnouncePhase <- function(z, x, subphase=NULL)
{
@@ -211,9 +173,9 @@
{
if (!is.batch())
{
- tkconfigure(z$tkvars$current, height=z$pp)
- tkconfigure(z$tkvars$deviation, height=z$pp)
- tkconfigure(z$tkvars$quasi, height=z$pp)
+ tkconfigure(z$tkvars$current, height=min(z$pp, 30))
+ tkconfigure(z$tkvars$deviation, height=min(z$pp, 30))
+ tkconfigure(z$tkvars$quasi, height=min(z$pp, 30))
}
n1pos <- z$n1 * (z$pp + 1)
z$n2min0 <- 7 + z$pp
@@ -276,21 +238,29 @@
w
}
+##@WriteOutTheta siena07 Progress reporting
+WriteOutTheta <- function(z)
+{
+ if (!is.batch())
+ {
+ DisplayTheta(z)
+ }
+ else
+ {
+ Report(c("theta:", format(z$theta, digits=3), "\n"))
+ }
+ Report("Current parameter values:\n", cf)
+ Report(format(z$theta), cf, fill=80)
+}
+
##@DisplayThetaAutocor siena07 Progress reporting
DisplayThetaAutocor <- function(z)
{
if (!is.batch())
{
- tkdelete(z$tkvars$current, "1.0", "end")
- tmp<- paste(c("", rep("\n", z$pp - 1)),
- format(round(z$theta, 4), width=12, sep=""),
- collapse="")
- tkinsert(z$tkvars$current, "1.0", tmp)
+ DisplayTheta(z)
tkdelete(z$tkvars$quasi, "1.0", "end")
- tmp<- paste(c("", rep("\n", z$pp - 1)),
- format(round(z$ac, 4), width=12, sep=""),
- collapse="")
- tkinsert(z$tkvars$quasi, "1.0", tmp)
+ tkinsert(z$tkvars$quasi, "1.0", FormatString(z$pp, z$ac))
}
else
{
@@ -304,11 +274,7 @@
{
if (!is.batch())
{
- tkdelete(z$tkvars$current, "1.0", "end")
- tmp<- paste(c("", rep("\n", z$pp - 1)),
- format(round(z$theta, 4), width=12, nsmall=4, sep=""),
- collapse="")
- tkinsert(z$tkvars$current, "1.0", tmp)
+ DisplayTheta(z)
}
else
{
@@ -318,27 +284,31 @@
##@DisplayTheta siena07 Progress reporting
DisplayTheta <- function(z)
{
- if (!is.batch())
- {
- tkdelete(z$tkvars$current, "1.0", "end")
- tmp<- paste(c("", rep("\n", z$pp - 1)),
- format(round(z$theta, 4), width=12, sep="", nsmall=4),
- collapse="")
- tkinsert(z$tkvars$current, "1.0", tmp)
- }
+ if (!is.batch())
+ {
+ tkdelete(z$tkvars$current, "1.0", "end")
+ tkinsert(z$tkvars$current, "1.0", FormatString(z$pp, z$theta))
+ }
}
+##@FormatString siena07 Progress Reporting
+FormatString <- function(pp, value)
+{
+ ppuse <- min(30, pp)
+ nbrs <- format(1:ppuse)
+ nch <- nchar(nbrs[1])
+ formatstr <- paste("%", nch, "d.%", (13 - nch), ".4f\n", sep="",
+ collapse="")
+ paste(sprintf(formatstr, 1:ppuse, value[1:ppuse]), collapse="")
+}
##@DisplayDeviations siena07 Progress reporting
DisplayDeviations <- function(z, fra)
{
- if (!is.batch())
- {
- tkdelete(z$tkvars$deviations, "1.0", "end")
- tmp<- paste(c("", rep("\n", z$pp - 1)),
- format(round(fra, 4), width=12, sep="", nsmall=4),
- collapse="")
- tkinsert(z$tkvars$deviations, "1.0", tmp)
- }
+ if (!is.batch())
+ {
+ tkdelete(z$tkvars$deviations, "1.0", "end")
+ tkinsert(z$tkvars$deviations, "1.0", FormatString(z$pp, fra))
+ }
}
##@DisplayIteration siena07 Progress reporting
DisplayIteration <- function(z)
Modified: pkg/RSiena/R/sienaprint.r
===================================================================
--- pkg/RSiena/R/sienaprint.r 2009-11-08 17:08:39 UTC (rev 22)
+++ pkg/RSiena/R/sienaprint.r 2009-11-08 17:11:12 UTC (rev 23)
@@ -24,6 +24,7 @@
sapply(x$nodesets,length)))
print(tmp)
}
+ invisible(x)
}
##@print.sienaGroup Methods
print.sienaGroup <- function(x, ...)
@@ -35,6 +36,7 @@
cat(paste(att$netnames, ":", att$types),'\n')
cat('Total number of periods:', att$observations)
cat("\nmore to be added!\n")
+ invisible(x)
}
##@print.sienafit Methods
@@ -88,6 +90,7 @@
cat(" \n*** Warning ***",
"Estimation terminated early at user request.\n")
}
+ invisible(x)
}
##@summary.sienaFit Methods
@@ -120,6 +123,7 @@
covcor[lower.tri(covcor)] <- correl[lower.tri(correl)]
printMatrix(format(round(t(covcor),digits=3),width=12))
}
+ invisible(x)
}
##@printMatrix Miscellaneous
@@ -157,6 +161,7 @@
if (x$condvarno > 0)
cat('conditioned on First variable')
}
+ invisible(x)
}
##@sienaFitThetaTable Miscellaneous
Modified: pkg/RSiena/R/simstatsc.r
===================================================================
--- pkg/RSiena/R/simstatsc.r 2009-11-08 17:08:39 UTC (rev 22)
+++ pkg/RSiena/R/simstatsc.r 2009-11-08 17:11:12 UTC (rev 23)
@@ -343,7 +343,8 @@
if (sum(z$test))
{
z$fra <- colMeans(z$sf, na.rm=TRUE)
- z <- ScoreTest(z, x)
+ ans <- ScoreTest(z$pp, z$dfra, z$msf, z$fra, z$test, x$maxlike)
+ z <- c(z, ans)
TestOutput(z, x)
}
dimnames(z$dfra)[[1]] <- as.list(z$effects$shortName)
@@ -353,6 +354,22 @@
f <- FRANstore()
# browser()
# cat(f$randomseed2, f$storedseed, '\n')
+ if (fromFiniteDiff)
+ {
+ returnDeps <- FALSE
+ }
+ else
+ {
+ returnDeps <- f$returnDeps
+ }
+ if (is.null(f$seeds))
+ {
+ seeds <- NULL
+ }
+ else
+ {
+ seeds <- f$seeds
+ }
if (is.null(f$randomseed2))
{
randomseed2 <- NULL
@@ -371,13 +388,14 @@
## cat(randomseed2, '\n')
}
ans <- .Call('model', PACKAGE="RSiena",
- z$Deriv, f$pData, f$seeds,
+ z$Deriv, f$pData, seeds,
fromFiniteDiff, f$pModel, f$myeffects, z$theta,
- randomseed2, f$returnDeps, z$FinDiff.method)
+ randomseed2, returnDeps, z$FinDiff.method)
# browser()
if (!fromFiniteDiff)
{
- f$seeds <- ans[[3]]
+ if (z$FinDiff.method)
+ f$seeds <- ans[[3]]
}
if (z$Deriv)
{
@@ -391,8 +409,11 @@
fra <- t(ans[[1]])
f$randomseed2 <- ans[[5]]#[c(1,4,3,2)]
FRANstore(f)
- sims <- ans[[6]]
- if (returnDeps)
+ if (f$returnDeps)
+ sims <- ans[[6]]
+ else
+ sims <- NULL
+ if (f$returnDeps)
{
## attach the names
names(sims) <- names(f)[1:length(sims)]
Modified: pkg/RSiena/src/siena07.cpp
===================================================================
--- pkg/RSiena/src/siena07.cpp 2009-11-08 17:08:39 UTC (rev 22)
+++ pkg/RSiena/src/siena07.cpp 2009-11-08 17:11:12 UTC (rev 23)
@@ -553,7 +553,7 @@
{
error("cannot find effect3");
}
-// Rprintf("%d parmcol\n", *parmCol);
+//Rprintf("%d parmcol\n", *parmCol);
}
@@ -754,7 +754,7 @@
setupOneModeNetwork(VECTOR_ELT(ONEMODES, period),
pOneModeNetworkLongitudinalData,
period);
- }
+ }
UNPROTECT(2);
}
/**
@@ -1006,6 +1006,7 @@
actorSet, 0)));
BehaviorLongitudinalData * pBehaviorData =
pData->createBehaviorData(CHAR(STRING_ELT(name, 0)), myActorSet);
+// Rprintf("%x\n", pBehaviorData);
setupBehavior(VECTOR_ELT(BEHGROUP, behavior), pBehaviorData);
UNPROTECT(2);
}
@@ -2662,9 +2663,9 @@
{
SET_VECTOR_ELT(ans, 1, scores);
}
- if (deriv || (!fromFiniteDiff))
+ if (returnDependents)
{
- SET_VECTOR_ELT(ans, 5, sims);/* not done in phase 2 */
+ SET_VECTOR_ELT(ans, 5, sims);/* not done in phase 2 !!!!test properly*/
}
SET_VECTOR_ELT(ans, 0, fra);
SET_VECTOR_ELT(ans, 3, ntim);
More information about the Rsiena-commits
mailing list