STA313 F 2004 Handout 4

Simple Regression SEM with R

# simplereg.R               R --vanilla < simplereg.R
#                        Simple regression SEM
#
# Simulate some data
set.seed(444)
x <- rnorm(250,0,2) ; y <- rnorm(250,0,3)
XY <- cbind(x,y)
#
# Define functions
#
source("mvnLL.R")
simplereg1 <- function(theta,datta)  # Full (unrestricted) moel
    { # theta = (sigmaXX,sigmaUU,b)
    # Elements of the covariance matrix are functions of the parameter vector
    sig11 <- theta[1] ; sig12 <- theta[3]*theta[1]
    sig22 <- theta[3]^2*theta[1] + theta[2]
    # Bind the covariance elements into a matrix
    covmat <- rbind( c(sig11,sig12),
                     c(sig12,sig22)   )
    simplereg1 <- mvnLL(covmat,datta)
    simplereg1
    } # End of simplereg1

simplereg0 <- function(theta,datta)  # Reduced (restricted) moel
    { # theta = (sigmaXX,sigmaUU,b)
    # Elements of the covariance matrix are functions of the parameter vector
    sig11 <- theta[1] ; sig22 <- theta[2]
    # Bind the covariance elements into a matrix
    covmat <- rbind( c(sig11,    0),
                     c(    0,sig22)   )
    simplereg0 <- mvnLL(covmat,datta)
    simplereg0
    } # End of simplereg0

# Get starting values for the full model
Shat <- var(XY) ; Shat
# start1 contins MOM estimates for theta = (sigmaXX,sigmaUU,b)
start1 <- numeric(3)
start1[1] <- Shat[1,1]
start1[2] <- Shat[2,2] - Shat[1,2]^2/Shat[1,1]
start1[3] <- Shat[1,2]/Shat[1,1]
start1

Full <- nlm(simplereg1,start1,datta=XY) ; Full
Red  <- nlm(simplereg0,c(var(x),var(y)),datta=XY) ; Red

G <- Red$minimum - Full$minimum
pval <- 1-pchisq(G,1)
G ; pval

# Compare common regression through origin
summary(lm(y~-1+x))

# Use the regression to get nice starting values

var(x)
start2 <- c(4.021413,3.047^2,-0.04720)
rbind(start1,start2)
Full2 <- nlm(simplereg1,start2,datta=XY)
rbind(Full$estimate,Full2$estimate)



> # simplereg.R               R --vanilla < simplereg.R
> #                        Simple regression SEM
> #
> # Simulate some data
> set.seed(444)
> x <- rnorm(250,0,2) ; y <- rnorm(250,0,3)
> XY <- cbind(x,y)
> #
> # Define functions
> #
> source("mvnLL.R")
> simplereg1 <- function(theta,datta)  # Full (unrestricted) moel
+     { # theta = (sigmaXX,sigmaUU,b)
+     # Elements of the covariance matrix are functions of the parameter vector
+     sig11 <- theta[1] ; sig12 <- theta[3]*theta[1]
+     sig22 <- theta[3]^2*theta[1] + theta[2]
+     # Bind the covariance elements into a matrix
+     covmat <- rbind( c(sig11,sig12),
+                      c(sig12,sig22)   )
+     simplereg1 <- mvnLL(covmat,datta)
+     simplereg1
+     } # End of simplereg1
>
> simplereg0 <- function(theta,datta)  # Reduced (restricted) moel
+     { # theta = (sigmaXX,sigmaUU,b)
+     # Elements of the covariance matrix are functions of the parameter vector
+     sig11 <- theta[1] ; sig22 <- theta[2]
+     # Bind the covariance elements into a matrix
+     covmat <- rbind( c(sig11,    0),
+                      c(    0,sig22)   )
+     simplereg0 <- mvnLL(covmat,datta)
+     simplereg0
+     } # End of simplereg0
>
> # Get starting values for the full model
> Shat <- var(XY) ; Shat
           x          y
x  4.0214134 -0.1101015
y -0.1101015  9.1585168
> # start1 contins MOM estimates for theta = (sigmaXX,sigmaUU,b)
> start1 <- numeric(3)
> start1[1] <- Shat[1,1]
> start1[2] <- Shat[2,2] - Shat[1,2]^2/Shat[1,1]
> start1[3] <- Shat[1,2]/Shat[1,1]
> start1
[1]  4.02141339  9.15550239 -0.02737881
>
> Full <- nlm(simplereg1,start1,datta=XY) ; Full
$minimum
[1] 2318.432

$estimate
[1]  4.00532405  9.11887424 -0.02737932

$gradient
[1] -2.611322e-05 -4.687668e-06 -2.728484e-06

$code
[1] 1

$iterations
[1] 6

> Red  <- nlm(simplereg0,c(var(x),var(y)),datta=XY) ; Red
$minimum
[1] 2318.514

$estimate
[1] 4.005338 9.121846

$gradient
[1]  1.938048e-04 -9.726234e-05

$code
[1] 1

$iterations
[1] 6

>
> G <- Red$minimum - Full$minimum
> pval <- 1-pchisq(G,1)
> G ; pval
[1] 0.08229895
[1] 0.7742058
>
> # Compare common regression through origin
> summary(lm(y~-1+x))

Call:
lm(formula = y ~ -1 + x)

Residuals:
    Min      1Q  Median      3Q     Max
-6.6348 -1.8367  0.3149  2.2161  9.3319

Coefficients:
  Estimate Std. Error t value Pr(>|t|)
x -0.04720    0.09568  -0.493    0.622

Residual standard error: 3.047 on 249 degrees of freedom
Multiple R-Squared: 0.0009766,  Adjusted R-squared: -0.003036
F-statistic: 0.2434 on 1 and 249 DF,  p-value: 0.6222

>
> # Use the regression to get nice starting values
>
> var(x)
[1] 4.021413
> start2 <- c(4.021413,3.047^2,-0.04720)
> rbind(start1,start2)
           [,1]     [,2]        [,3]
start1 4.021413 9.155502 -0.02737881
start2 4.021413 9.284209 -0.04720000
> Full2 <- nlm(simplereg1,start2,datta=XY)
> rbind(Full$estimate,Full2$estimate)
         [,1]     [,2]        [,3]
[1,] 4.005324 9.118874 -0.02737932
[2,] 4.005329 9.118906 -0.02737937
>