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.