iwlsm: Function to fit an iterated weighted least squares model.

In RSiena: Siena - Simulation Investigation for Empirical Network Analysis

Function to fit an iterated weighted least squares model.

Description

Fits an iterated weighted least squares model.

Usage

iwlsm(x, ...)

## S3 method for class 'formula'
iwlsm(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)

## Default S3 method:
iwlsm(x, y, weights, ses, ..., w = rep(1/nrow(x), 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)

psi.iwlsm(u, k, deriv = 0, w, sj2, hh)

Arguments

formula	a formula of the form y ~ x1 + x2 + ....
data	data frame from which variables specified in formula are preferentially to be taken.
weights	a vector of prior weights for each case.
subset	An index vector specifying the cases to be used in fitting.
ses	Estimated variance of the responses. Will be paseed to psi as sj2
na.action	A function to specify the action to be taken if NAs are found. The ‘factory-fresh’ default action in R is na.omit, and can be changed by options(na.action=).
x	a matrix or data frame containing the explanatory variables.
y	the response: a vector of length the number of rows of x.
method	Must be "M". (argument not used here).
wt.method	are the weights case weights (giving the relative importance of case, so a weight of 2 means there are two of these) or the inverse of the variances, so a weight of two means this error is half as variable? This will not work at present.
model	should the model frame be returned in the object?
x.ret	should the model matrix be returned in the object?
y.ret	should the response be returned in the object?
contrasts	optional contrast specifications: se lm.
w	(optional) initial down-weighting for each case. Will not work at present.
init	(optional) initial values for the coefficients OR a method to find initial values OR the result of a fit with a coef component. Known methods are "ls" (the default) for an initial least-squares fit using weights w*weights, and "lts" for an unweighted least-trimmed squares fit with 200 samples. Probably not functioning.
psi	the psi function is specified by this argument. It must give (possibly by name) a function g(x, ..., deriv, w) that for deriv=0 returns psi(x)/x and for deriv=1 returns some value. Extra arguments may be passed in via ....
scale.est	method of scale estimation: re-scaled MAD of the residuals (default) or Huber"s proposal 2 (which can be selected by either "Huber" or "proposal 2").
k2	tuning constant used for Huber proposal 2 scale estimation.
maxit	the limit on the number of IWLS iterations.
acc	the accuracy for the stopping criterion.
test.vec	the stopping criterion is based on changes in this vector.
...	additional arguments to be passed to iwlsm.default or to the psi function.
lqs.control	An optional list of control values for lqs.
u	numeric vector of evaluation points.
k	tuning constant. Not used.
deriv	0 or 1: compute values of the psi function or of its first derivative. (Latter not used).
sj2	Estimated variance of the responses
hh	Diagonal values of the hat matrix

Details

This function is very slightly adapted from rlm in packages MASS. It alternates between weighted least squares and estimation of variance on the basis of a common variance. The function psi.iwlsm calculates the weights for the next iteration. Used by siena08 to combine estimates from different sienaFits.

Value

An object of class "iwlsm" inheriting from "lm". Note that the df.residual component is deliberately set to NA to avoid inappropriate estimation of the residual scale from the residual mean square by "lm" methods.

The additional components not in an lm object are

s	the robust scale estimate used
w	the weights used in the IWLS process
psi	the psi function with parameters substituted
conv	the convergence criteria at each iteration
converged	did the IWLS converge?
wresid	a working residual, weighted for "inv.var" weights only.

Note

The function has been changed as little as possible, but has only been used with default arguments. The other options have been retained just in case they may prove useful.

Author(s)

Ruth Ripley

References

Venables, W. N. and Ripley, B. D. (2002) Modern Applied Statistics with S. Fourth edition. Springer. See also https://www.stats.ox.ac.uk/~snijders/siena/

See Also

siena08, sienaMeta, sienaFit

Examples

## Not run: 
##not enough data here for a sensible example, but shows the idea.
myalgorithm <- sienaAlgorithmCreate(nsub=2, n3=100)
mynet1 <- sienaDependent(array(c(s501, s502), dim=c(50, 50, 2)))
mynet2 <- sienaDependent(array(c(s502, s503), dim=c(50, 50, 2)))
mydata1 <- sienaDataCreate(mynet1)
mydata2 <- sienaDataCreate(mynet2)
myeff1 <- getEffects(mydata1)
myeff2 <- getEffects(mydata2)
myeff1 <- setEffect(myeff1, transTrip, fix=TRUE, test=TRUE)
myeff2 <- setEffect(myeff2, transTrip, fix=TRUE, test=TRUE)
myeff1 <- setEffect(myeff1, cycle3, fix=TRUE, test=TRUE)
myeff2 <- setEffect(myeff2, cycle3, fix=TRUE, test=TRUE)
ans1 <- siena07(myalgorithm, data=mydata1, effects=myeff1, batch=TRUE)
ans2 <- siena07(myalgorithm, data=mydata2, effects=myeff2, batch=TRUE)
meta <- siena08(ans1, ans2)
metadf <- split(meta$thetadf, meta$thetadf$effects)[[1]]
metalm <- iwlsm(theta ~ tconv, metadf, ses=se^2)

## End(Not run)

RSiena documentation built on Nov. 2, 2023, 5:19 p.m.