Nothing
starting.values <- function(family, link, vlink, rn, model, Kw, X, y, id, cl, start, other){
# the relevant length for the start vector provided by the user is
# either : 0 (nothing), K (the length of the beta vector or the
# number of the additional parameter (1 most of the time, but 2 for
# the negbin model)
ls <- length(start)
K <- ncol(X)
if (family == "binomial"){
if (model == "pooling"){
# use Phi^-1(ybar) as starting value for the intercept and 0 for
# the other coefficients
if (!ls %in% c(0, K)) stop("irrelevant length for the start vector")
if (ls == 0){
interc <- ifelse(link == "probit", qnorm(mean(y)), - log(1 / mean(y) - 1))
start <- c(interc, rep(0, ncol(X) - 1))
names(start) <- colnames(X)
}
}
if (model == "random"){
if (!ls %in% c(0, 1, K + 1)) stop("irrelevant length for the start vector")
if (ls <= 1){
# the case where start is null or of length one (the sigma
# parameter). In these cases, compute the pooling model
startcl <- cl
startcl[c('model', 'method', 'print.level')] <- c('pooling', 'nr', 0)
glmest <- eval(startcl, parent.frame())
if (ls == 0){
# when start is null compute the starting value of sigma
sigma <- lnl.binomial(coef(glmest), y = y, X = X, id = id,
link = link, model = "pooling", start.sigma = TRUE)
start <- c(coef(glmest), sigma = sigma)
}
else start <- c(coef(glmest), sigma = start)
}
}
}
if (family == "ordinal"){
J <- length(unique(y))
if (model == "pooling"){
if (!ls %in% c(0, K + J - 2)) stop("irrelevant length for the start vector")
if (ls == 0){
fy <- prop.table(table(y))
fyc <- qnorm(cumsum(fy))
beta0 <- - fyc[1]
sup.coef <- fyc[2:(J-1)] + beta0
start <- c(beta0, rep(0, K-1), sup.coef)
}
names(start) <- c(colnames(X), paste("mu", 1:(J - 2), sep = "_"))
}
if (model == "random"){
sup.coef.names <- "sigma"
if (!ls %in% c(0, 1, K + J - 1)) stop("irrelevant length for the start vector")
if (ls <= 1){
# the case where start is null or of length one (the sigma
# parameter). In these cases, compute the pooling model
startcl <- cl
startcl[c('model', 'method', 'print.level')] <- c('pooling', 'nr', 0)
glmest <- eval(startcl, parent.frame())
if (ls == 0){
# when start is null compute the starting value of sigma
sigma <- lnl.ordinal(coef(glmest), y = y, X = X, id = id,
link = link, model = "pooling", start.sigma = TRUE)
start <- c(coef(glmest), sigma = sigma)
}
else start <- c(coef(glmest), sigma = start)
}
}
}
if (family == "poisson"){
if (model == "pooling"){
# use Phi^-1(ybar) as starting value for the intercept
if (!ls %in% c(0, K)) stop("irrelevant length for the start vector")
if (ls == 0){
start <- c(log(mean(y)), rep(0, K - 1))
names(start) <- colnames(X)
}
}
else{
if (model == "within") relevant.values <- c(0, Kw)
else relevant.values <- c(0, 1, K + 1)
if (!ls %in% relevant.values) stop("irrelevant length for the start vector")
if (ls <= 1){
# the case where start is null or of length one (the sigma
# parameter). In these cases, compute the pooling model
startcl <- cl
startcl$start <- NULL
startcl[c('model', 'method', 'print.level')] <- c('pooling', 'nr', 0)
glmest <- eval(startcl, parent.frame())
if (ls == 0){
if (model == "within") start <- coef(glmest)[Kw]
else{
if (FALSE){
big <- 100
init.val <- ifelse(other == "inv", big, 1/big)
haty <- lnl.poisson(coef(glmest), y = y, X = X, id = id,
link = link, model = "pooling")
haty <- attr(haty, "fitted.values")
print(head(haty))
alphait <- y / haty
alphai <- tapply(alphait, id, mean)
print(head(alphai))
print(var(alphai));stop()
sigma <- lnl.poisson(c(coef(glmest), init.val), y = y, X = X, id = id,
link = link, model = model)
g <- apply(attr(sigma, "gradient"), 2, sum)
g <- g[length(g)]
h <- attr(sigma, "hessian")
h <- h[nrow(h), nrow(h)]
sigma <- init.val - g / h
start <- c(coef(glmest), sigma = sigma)
}
start <- c(coef(glmest), sigma = 1.2)
}
}
else start <- c(coef(glmest), sigma = start)
}
}
}
if (family == "negbin"){
if (model == "pooling"){
# Use the Poisson model for starting values
startcl <- cl
startcl[c('model', 'method', 'print.level', 'family')] <- c('pooling', 'nr', 0, poisson)
glmest <- eval(startcl, parent.frame())
haty <- lnl.poisson(coef(glmest), y = y, X = X, id = id,
link = link, model = "pooling")
haty <- attr(haty, "fitted.values")
res <- y - haty
V <- var(res)
E <- mean(y)
if (vlink == 'nb1'){
alpha <- V / E - 1
}
if (vlink == 'nb2'){
alpha <- V / E^2 - 1 / E
## print(alpha)
## print(summary(lm(I(res^2/haty) ~ haty )))
## print(summary(lm(I(res^2/haty) ~ haty -1)))
## alpha <- coef(lm(I(res^2/haty)~ haty - 1)) - 1
## print(alpha)
}
start <- c(coef(glmest), sigma = alpha)
}
else{
if (model != "within") relevant.values <- c(0, K + 1L)
else relevant.values <- c(0, 2, K + 3L)
if (!ls %in% relevant.values) stop("irrelevant length for the start vector")
if (ls <= 1){
# the case where start is null or of length one (the sigma
# parameter). In these cases, compute the pooling model
startcl <- cl
startcl[c('model', 'method', 'print.level', 'family')] <- c('pooling', 'nr', 0, poisson)
glmest <- eval(startcl, parent.frame())
if (ls == 0){
if (model == "within") start <- coef(glmest)#[within.var]
else start <- c(coef(glmest), a = 2, b = 1)
}
else start <- c(coef(glmest), a = 2, b = 1)
}
}
}
if (family == "gaussian"){
if (is.null(other)) other <- "sd"
if (is.null(start)){
if (FALSE){
startcl <- cl
startcl[[1]] <- as.name("plm")
startcl$family<- NULL
startcl$model <- "random"
glmest <- eval(startcl, parent.frame())
theta <- glmest$ercomp$theta
sigma2 <- glmest$ercomp$sigma2$idios
sigmu <- glmest$ercomp$sigma2$id^0.5
sigeps <- glmest$ercomp$sigma2$idiod^0.5
}
else{
m <- match(c("formula", "data", "subset", "na.action"),names(cl),0)
lmcl <- cl[c(1,m)]
lmcl[[1]] <- as.name("lm")
lmcl <- eval(lmcl, parent.frame())
eb <- tapply(resid(lmcl), id, mean)[as.character(id)]
sig2mu <- var(eb)
ew <- resid(lmcl) - eb
sig2eta <- var(ew)
gamma <- sig2mu / sig2eta
}
if (other == "sd"){
start <- c(coef(lmcl), sd.id = sqrt(sig2mu), sd.idios = sqrt(sig2eta))
}
else{
start <- c(coef(lmcl), gamma = gamma, sig2eta = sig2eta)
}
}
}
if (family == "tobit"){
if (model == "pooling"){
if (!ls %in% c(0, K + 1)) stop("irrelevant length for the start vector")
if (ls == 0){
m <- match(c("formula", "data", "subset", "na.action"),names(cl),0)
lmcl <- cl[c(1,m)]
lmcl[[1]] <- as.name("lm")
lmcl <- eval(lmcl, parent.frame())
sig2 <- deviance(lmcl) / df.residual(lmcl)
if (is.null(other)) other <- "sd"
if (other == "var") sigma <- sig2
if (other == "sd") sigma <- sqrt(sig2)
if (other == "lsd") sigma <- log(sqrt(sig2))
start <- c(coef(lmcl), sd.idios = sigma)
}
}
else{
if (ls <= 1){
startcl <- cl
startcl$model <- "pooling"
pglmest <- eval(startcl, parent.frame())
thestart <- coef(pglmest)
if (ls == 1){
start <- c(thestart, start)
}
else{
sigma <- lnl.tobit(coef(pglmest), y = y, X = X, id = id,
link = link, model = "pooling", start.sigma = TRUE)
start <- c(thestart, sd.id = sigma)
}
}
}
}
start
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.