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.