R version 2.9.0 (2009-04-17)
Copyright (C) 2009 The R Foundation for Statistical Computing
ISBN 3-900051-07-0


>  # cinlls3.r                  November 2007, 2008, 2009
>  #
>  # simulation similar to Chapter 9, Exercise 24
>  #
>  # nonlinear least squares problem nllsqt1
>  #  test problem from Fuller(1976)
>  #
>  # compare different methods for se's in nonlinear regression
>  #  
>   n <- 12
>   p <- 3
>   x <- c(0, .4, .8, 1.6, 0, 1, 2, 4, 0, 1.5, 3, 5.9) 
>   x
 [1] 0.0 0.4 0.8 1.6 0.0 1.0 2.0 4.0 0.0 1.5 3.0 5.9
>  # no y's -- we'll be generating them
>  #
>  # function to compute everything: estimates, se's
>  #
>   methods <- function(y) {   # nonlinear least squares confidence int
+  # create data frame or list for nls
+   fulnls <- list(y,x)
+  #  find better starting values
+   t1 <- y[ x==max(x)]             # max value as asymptote
+   xs <- x[ x < 2 ]                # use smallest x's
+   ys <- y[ x < 2 ]
+   ok <- ( t1 > ys )               # second subsetting
+   xs <- xs[ok]
+   ys <- log( t1 - ys[ok] )        # transform reversed y's
+  #  simple linear regression
+   t3 <- sum( (xs - mean(xs) )*ys )/ sum( (xs-mean(xs))**2 ) # slope
+   t2 <- mean(ys)-t3*mean(xs)                                # intercept
+  #  transform back from log and flip signs
+   beta0 <- c(t1,-exp(t2),-t3)
+  # sum of squares function
+   Sbeta <- function(t) {      # sum of squares function
+     e <- y - t[1] - t[2]*exp(-t[3]*x)
+     Sbeta <- sum(e*e)  }
+  # call nls
+   soln <- nls( y ~ th1 + th2*exp(-th3*x),
+         start=list(th1=beta0[1],th2=beta0[2],th3=beta0[3]) )
+   beta <- coef(soln)
+  # different methods for standard errors for coefficient estimates
+     e <- y - beta[1] - beta[2]*exp(-beta[3]*x)    # residuals
+   sse <- sum(e*e)                                 # error sum of squares
+  # gradient matrix -- works as design matrix
+   G <- matrix( c(rep(1,n),exp(-beta[3]*x),-beta[2]*x*exp(-beta[3]*x)),
+           n,p)
+  # meat for sandwich
+   cs <- (n-1)*cov(G*e)
+  # using other pieces -- first 2nd partial of gi * residuali
+  #  start with zero and just do the nonzero terms
+   del2ge <- matrix(0,p,p)
+   del2ge[2,3] <- sum(e*(-x*exp(-beta[3])) )     # d2 g / db2 db3
+   del2ge[3,2] <- del2ge[2,3]                    # symmetry
+   del2ge[3,3] <- sum(e*(beta[2]*x*x*exp(-beta[3]*x))) # d2 g / db3 db3
+   GG <- t(G) %*% G                      # analogue to X'X
+   hdel2S <- GG - del2ge                 # half of hessian of S(beta)
+  #
+  # all of the different methods
+   GGinv <- chol2inv(chol(GG))
+   hdel2Sinv <- chol2inv(chol(hdel2S))
+   m1 <- diag((sse/(n-p))*GGinv)           # sigma2hat*inv(G'G)
+   m2 <- diag((sse/(n-p))*hdel2Sinv)       # using hessian
+   m3 <- diag(GGinv%*%cs%*%GGinv)          # easy sandwich
+   m4 <- diag(hdel2Sinv%*%cs%*%hdel2Sinv)  # harder sandwich
+  # finish function by putting all together
+   methods <- c(beta,sse,sqrt(c(m1,m2,m3,m4)))    } # end of function
>  #
>  # simulation experiment
>   set.seed(5151917)                # set seed
>   nrep <- 1000                     # number of replications
>   betatrue <- c(150, -80, 1.5)     # true value of beta
>   sigtrue <- 8                     # true noise level
>  # mean vector is constant for all replications
>   yhat <- betatrue[1] + betatrue[2]*exp(-betatrue[3]*x)
>  #
>  # generate Y's
>   Ys <- matrix(rep(yhat,nrep)+sigtrue*rnorm(nrep*n),n,nrep)
>  #
>  # *** main step *** compute estimates, etc for each sample ***
>   ests <- apply(Ys,2,methods)
>  #
>  # difference of estimates from true
>   diffs <- ests[1:3,]-betatrue
>  # sampling variation of estimates
>   bias <- apply(diffs,1,mean)
>   bias
        th1         th2         th3 
 0.35920277 -0.45479742  0.05170761 
>   cov <- var(t(diffs))
>  # is bias significant?
>   sqrt(diag(cov)/nrep)        # statistically
       th1        th2        th3 
0.13459991 0.18562926 0.01100742 
>   bias/betatrue               # practically
        th1         th2         th3 
0.002394685 0.005684968 0.034471741 
>  # how do the four methods compare -- average se
>   table <- matrix(c(sqrt(diag(cov)),apply(ests[5:16,],1,mean)),3,5)
>   table     
          [,1]      [,2]      [,3]      [,4]      [,5]
[1,] 4.2564228 4.2325788 4.2489003 3.4731201 3.5234068
[2,] 5.8701128 5.8755030 5.9127910 4.8915843 4.9787759
[3,] 0.3480853 0.3279204 0.3313207 0.2677572 0.2774447
>  #  sd(bhat) mean(se1) mean(se2) mean(se3) mean(se4)
>  #
>  #                  now look at coverage of confidence intervals
>  #
>  # method 1 first
>   cvr1 <- apply( (abs(diffs) < 1.96*ests[5:7]),1,sum)
>  # method 2
>   cvr2 <- apply( (abs(diffs) < 1.96*ests[8:10]),1,sum)
>  # method 3
>   cvr3 <- apply( (abs(diffs) < 1.96*ests[11:13]),1,sum)
>  # method 4 last
>   cvr4 <- apply( (abs(diffs) < 1.96*ests[14:16]),1,sum)
>  #
>  # put in table -- method as column, parameter as row
>   matrix(c(cvr1,cvr2,cvr3,cvr4),3,4)
     [,1] [,2] [,3] [,4]
[1,]  982  986  986  995
[2,]  968  972  945  960
[3,]  778  803  743  785
>  # done
>   rm(list=ls())
>   q()