inst/doc/gelS4.R

In momentfit: Methods of Moments

## ----echo=FALSE---------------------------------------------------------------
library(knitr)
opts_chunk$set(size='footnotesize')

## -----------------------------------------------------------------------------
library(momentfit)
data(simData)
lin <- momentModel(y~x1+x2, ~x2+z1+z2, data=simData, vcov="iid")  

## -----------------------------------------------------------------------------
theta0=c(mu=1,sig=1)
x <- simData$x1
dat <- data.frame(x=x, x2=x^2, x3=x^3, x4=x^4)
gform <- list(x~mu,
              x2~mu^2+sig, 
              x3~mu^3+3*mu*sig,
              x4~mu^4+6*mu^2*sig+3*sig^2)
form <- momentModel(gform, NULL, theta0=theta0, vcov="MDS", data=dat)
form

## -----------------------------------------------------------------------------
fct <- function(theta, x)
    cbind(x-theta[1], (x-theta[1])^2-theta[2],
          (x-theta[1])^3, (x-theta[1])^4-3*theta[2]^2)
dfct <- function(theta, x)
    {
        m1 <- mean(x-theta[1])
        m2 <- mean((x-theta[1])^2)
        m3 <- mean((x-theta[1])^3)
        matrix(c(-1, -2*m1, -3*m2, -4*m3, 
                 0, -1, 0, -6*theta[2]), 4, 2)
    }
theta0=c(mu=1,sig2=1)
func <- momentModel(fct, simData$x3, theta0=theta0, grad=dfct, vcov="iid")
func

## -----------------------------------------------------------------------------
theta0=c(b0=1, b1=1, b2=1)
gform <- y~exp(b0+b1*x1+b2*x2)
nlin <- momentModel(gform, ~z1+z2+z3+x2, theta0=theta0, vcov="MDS", data=simData)
nlin

## -----------------------------------------------------------------------------
rhoEL

## -----------------------------------------------------------------------------
args(getLambda)

## -----------------------------------------------------------------------------
X <- simData[c("x3","x4","z5")]
(res <- getLambda(X, gelType="EL"))$lambda
res$convergence$convergence

## -----------------------------------------------------------------------------
(res <- getLambda(X, gelType="ET",control=list(maxit=2000)))$lambda
res$convergence$convergence

## -----------------------------------------------------------------------------
(res <- getLambda(X, rhoFct=rhoEEL))$lambda
res$convergence$convergence

## -----------------------------------------------------------------------------
showMethods("solveGel")

## -----------------------------------------------------------------------------
mylSolve <- function(gmat, lambda0, gelType=NULL, rhoFct=NULL, k=1, ...) 
    {
        lambda <-  rep(0,ncol(gmat))
        obj <-  sum(colMeans(gmat)^2)
        list(lambda=lambda, convergence=0, obj=obj)
    }

## -----------------------------------------------------------------------------
solveGel(lin,theta0=c(0,0,0), lamSlv=mylSolve)

## -----------------------------------------------------------------------------
mylSolve <- function(gmat, lambda0=NULL, gelType=NULL, rhoFct=NULL, k=1, ...) 
    {
        gmat <- as.matrix(gmat)
        res  <-  getLambda(gmat, lambda0, gelType="ET", k=k)
        gml <- c(gmat%*%res$lambda)
        w <- exp(gml)
        w <- w/sum(w)
        n <- nrow(gmat)
        res$obj <- mean(-log(w*n))
        res
    }
etelFit <- solveGel(lin,theta0=c(1,1,0), lamSlv=mylSolve)
etelFit$theta

## -----------------------------------------------------------------------------
solveGel(update(lin, gelType="ETEL"), theta0=c(1,1,0))$theta

## ----warning=FALSE------------------------------------------------------------
solveGel(lin, theta0=c(1,1,0), lControl=list(algo="nlminb"))$theta

## -----------------------------------------------------------------------------
res <- solveGel(lin, theta0=c(1,1,0), lControl=list(restrictedLam=c(2L,3L)))
res$lambda

## -----------------------------------------------------------------------------
solveGel(lin, theta0=c(1,1,0),
         tControl=list(method="BFGS", control=list(maxit=2000, reltol=1e-9)))$theta

## ----warning=FALSE------------------------------------------------------------
fit <- gelFit(lin)

## -----------------------------------------------------------------------------
showClass("gelfit")

## -----------------------------------------------------------------------------
fit <- gelFit(lin)
print(fit, lambda=FALSE, model=FALSE)

## -----------------------------------------------------------------------------
pt <- getImpProb(fit)
pt$convMom

## -----------------------------------------------------------------------------
summary(fit)

## -----------------------------------------------------------------------------
specTest(fit)

## -----------------------------------------------------------------------------
confint(fit, parm=1:2, level=0.90)
confint(fit, level=0.90, lambda=TRUE)

## -----------------------------------------------------------------------------
confint(fit, 1:2, type="invLR", level=0.90)

## -----------------------------------------------------------------------------
muModel <- momentModel(x1~1, ~1, data=simData)

## -----------------------------------------------------------------------------
muFit <- gelFit(muModel, gelType="EL", tControl=list(method="Brent", lower=-10, upper=10))
muFit@theta

## -----------------------------------------------------------------------------
mean(simData$x1)

## -----------------------------------------------------------------------------
confint(muFit, type="invLR")

## -----------------------------------------------------------------------------
fit <- gelFit(lin)
fit@theta

## -----------------------------------------------------------------------------
fit2 <- update(fit, coefSlv="nlminb")
fit2@theta

## -----------------------------------------------------------------------------
fit3 <- update(fit2, gelType="ET", lControl=list())
fit3@theta

## -----------------------------------------------------------------------------
rlin <- restModel(fit3@model, "x1=1")
update(fit3, newModel=rlin)

## -----------------------------------------------------------------------------
mod1 <- momentModel(y~x1+x2, ~x2+z1+z2+z3, data=simData)
fit1 <- gelFit(mod1)

## -----------------------------------------------------------------------------
eta <- c(coef(fit1), fit1@lambda)
names(eta) <- NULL

## -----------------------------------------------------------------------------
mod2 <-  momentModel(momFct, fit1, theta0=eta, vcov="MDS")

## -----------------------------------------------------------------------------
fit2 <- evalGmm(mod2, eta)
v <- vcov(fit2)[1:3,1:3]
sqrt(diag(v))

## -----------------------------------------------------------------------------
sqrt(diag(vcov(fit1)$vcov_par))

## -----------------------------------------------------------------------------
sqrt(diag(vcov(fit1, robToMiss=TRUE)$vcov_par))

## -----------------------------------------------------------------------------
summary(fit1, robToMiss=TRUE)@coef

## -----------------------------------------------------------------------------
print(evalGel(lin, theta=c(1,2,3), lambda=rep(.1,4)), model=FALSE)

## -----------------------------------------------------------------------------
specTest(evalGel(muModel, theta=4), type="LR")

## -----------------------------------------------------------------------------
rmuModel <- restModel(muModel, R=matrix(1,1,1), rhs=4)
specTest(gelFit(rmuModel))

## -----------------------------------------------------------------------------
fit <- gel4(y~x1+x2, ~x2+z1+z2, data=simData, vcov="iid", gelType="EL",
             lControl=list(algo="Wu"))  
print(fit, model=FALSE)

## -----------------------------------------------------------------------------
fit <- gel4(y~x1+x2, ~x2+z1+z2, data=simData, vcov="iid", gelType="EL",
             lControl=list(algo="Wu"), theta0=c(1,1,0), initTheta="theta0")  
coef(fit)

## -----------------------------------------------------------------------------
fit <- gel4(y~x1+x2, ~x2+z1+z2, data=simData, vcov="iid", gelType="EL",
             lControl=list(algo="Wu"), cstLHS="x2=x1")
print(fit, model=FALSE)

## -----------------------------------------------------------------------------
theta0=c(mu=1,sig=1)
x <- simData$x1
dat <- data.frame(x=x, x2=x^2, x3=x^3, x4=x^4)
gform <- list(x~mu,
              x2~mu^2+sig, 
              x3~mu^3+3*mu*sig,
              x4~mu^4+6*mu^2*sig+3*sig^2)
fit <- gel4(g=gform, gelType="EEL", theta0=theta0, vcov="MDS", data=dat)
print(fit, model=FALSE)

## -----------------------------------------------------------------------------
fct <- function(theta, x) 
   cbind(x-theta[1], (x-theta[1])^2-theta[2],
          (x-theta[1])^3, (x-theta[1])^4-3*theta[2]^2)
dfct <- function(theta, x)
    {
        m1 <- mean(x-theta[1])
        m2 <- mean((x-theta[1])^2)
        m3 <- mean((x-theta[1])^3)
        matrix(c(-1, -2*m1, -3*m2, -4*m3, 
                 0, -1, 0, -6*theta[2]), 4, 2)
    }
fit <- gel4(g=fct, x=simData$x3, theta0=c(1,1), grad=dfct, vcov="iid", gelType="ET")
print(fit, model=FALSE)

## -----------------------------------------------------------------------------
theta0=c(b0=1, b1=0, b2=0)
gform <- y~exp(b0+b1*x1+b2*x2)
fit <- gel4(gform, ~z1+z2+z3+x2, theta0=theta0, vcov="MDS", data=simData,
                 gelType="HD")
print(fit, model=FALSE)

## -----------------------------------------------------------------------------
x2 <- simData$x2
confint(x2, gelType="EL")
print(confint(x2, gelType="EL", type="invLR"), digits=5)

## -----------------------------------------------------------------------------
confint(simData, "x2", gelType="EL")

## -----------------------------------------------------------------------------
res <- confint(simData, c("x2","y"), npoints=20)
res

## ----fig.height=5, eval=FALSE-------------------------------------------------
#  res2 <- confint(simData, c("x2","y"), type="invLR", npoints=20)
#  c1 <- col2rgb("darkorange2")/255
#  c1 <- rgb(c1[1],c1[2],c1[3],.5)
#  c2 <- col2rgb("lightblue2")/255
#  c2 <- rgb(c2[1],c2[2],c2[3],.5)
#  plot(res, pch=20, bg=1, Pcol=c1, col=c1, density=10, ylim=c(3.5,6.5),
#       xlim=c(4.8,7.5))
#  plot(res2, pch=20, bg=2, Pcol=c2, col=c2, density=10, add=TRUE)
#  legend("topright", c("Wald","LR"), col=c(c1,c2), pch=20 , bty="n")

## ----extract, message=FALSE, warning=FALSE------------------------------------
library(texreg)
setMethod("extract", "gelfit", 
          function(model, includeSpecTest=TRUE, 
                   specTest=c("LR","LM","J"), include.nobs=TRUE, 
                   include.obj.fcn=TRUE, ...)
              {
                  specTest <- match.arg(specTest)
                  s <- summary(model, ...)
                  wspecTest <- grep(specTest, rownames(s@specTest@test))
                  spec <- modelDims(model@model)
                  coefs <- s@coef
                  names <- rownames(coefs)
                  coef <- coefs[, 1]
                  se <- coefs[, 2]
                  pval <- coefs[, 4]
                  n <- model@model@n
                  gof <- numeric()
                  gof.names <- character()
                  gof.decimal <- logical()
                  if (includeSpecTest) {
                      if (spec$k == spec$q)
                          {
                              obj.fcn <- NA
                              obj.pv <- NA
                          } else {
                              obj.fcn <- s@specTest@test[wspecTest,1]
                              obj.pv <- s@specTest@test[wspecTest,3]
                          }
                      gof <- c(gof, obj.fcn, obj.pv)                      
                      gof.names <- c(gof.names, 
                                     paste(specTest,"-test Statistics", sep=""),
                                     paste(specTest,"-test p-value", sep=""))
                      gof.decimal <- c(gof.decimal, TRUE, TRUE)
                  }
                  if (include.nobs == TRUE) {
                      gof <- c(gof, n)
                      gof.names <- c(gof.names, "Num.\\ obs.")
                      gof.decimal <- c(gof.decimal, FALSE)
                  }
                  tr <- createTexreg(coef.names = names, coef = coef, se = se, 
                                     pvalues = pval, gof.names = gof.names, gof = gof, 
                                     gof.decimal = gof.decimal)
                  return(tr)
              })

## ----results='asis'-----------------------------------------------------------
fit1 <- gel4(y~x1, ~x2+x3+x4, data=simData)
fit2 <- gel4(y~x1+x2, ~x2+x3+x4, data=simData)
texreg(list(fit1,fit2), digits=4)