STA313 F 2004 Handout 9: Path Model 2a with R

Simple regression with measurement error: A non-identified model


# path2.R
# source("path2.R")
sigsqf <- 1 ; sigsqe <- 2 ; sigsqe3 <- 3 ; b <- -0.10
n <- 200 ; set.seed(9999)
F <- rnorm(n,0,sqrt(sigsqf)) ; e1 <- rnorm(n,0,sqrt(sigsqe)) 
e2 <- rnorm(n,0,sqrt(sigsqe)) ; e3 <- rnorm(n,0,sqrt(sigsqe3))

x1 <- F + e1
x2 <- F + e2
y <- b*F + e3

path2dat <- cbind(x1,x2,y)
write(t(path2dat),"path2.dat",ncolumns=3)
xydat <- cbind(x1,y)

source("mvnLL.R")
path2a <- function(theta,datta) # Full (unresticted) model
      #             1        2       3      4
    { # theta = (sigsqf, sigsqe, sigsqe3,   b) 4 parameters
      # Elements of the covariance matrix SIGij are functions of the
      # parameter vector.
      SIG11 <- theta[1]+theta[2];  SIG12 <- theta[4]*theta[1]
      SIG22 <- theta[4]^2*theta[1] + theta[3]
      covmat <- rbind( c(SIG11,SIG12),
                       c(SIG12,SIG22)   )
    path2a <-  mvnLL(covmat,datta)
    path2a
    } # End of path2a

start1 <- c(sigsqf,sigsqe,sigsqe3,b) # Start at the truth!
path2.1 <- nlm(path2a,start1,datta=xydat) ; print(path2.1)
cat ("\n")
cat ("---------------------------------------------------\n")
cat ("\n")
start2 <- start1*.9   # Still pretty close to the truth
path2.2 <- nlm(path2a,start2,datta=xydat) ; print(path2.2)


> source("path2.R")
$minimum
[1] 1579.282

$estimate
[1] 0.97578652 1.98931682 3.10719304 0.01092407

$gradient
[1] -3.637979e-06  8.115113e-06 -3.439298e-06 -2.273737e-07

$code
[1] 1

$iterations
[1] 5


---------------------------------------------------

$minimum
[1] 1579.282

$estimate
[1] 1.02417245 1.94092628 3.10720375 0.01040338

$gradient
[1] -1.098936e-04 -9.688056e-05  1.030322e-04 -2.064553e-04

$code
[1] 1

$iterations
[1] 6

Two different answers, same minimum value of -2LL, and no hint of any trouble. This is scary.