R : Copyright 2004, The R Foundation for Statistical Computing Version 1.9.1 (2004-06-21), ISBN 3-900051-00-3 R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R in publications. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for a HTML browser interface to help. Type 'q()' to quit R. > # > # exer5-Q1.q > # > #------------------------------------------------------------ > # > # PURPOSE: compare naive and correct standard errors > # > # AUTHOR: P. Heagerty > # > # DATE: 08 Feb 2007 > # > #------------------------------------------------------------ > # > X.a.0 <- cbind( 1, c(1:12), 0 ) > X.a.1 <- cbind( 1, c(1:12), 1 ) > X.a <- rbind( X.a.0, X.a.1 ) > # > X.b.0 <- cbind( 1, c(1:12), c(rep(0,6),rep(1,6)) ) > X.b.1 <- cbind( 1, c(1:12), c(rep(1,6),rep(0,6)) ) > X.b <- rbind( X.b.0, X.b.1 ) > # > ### > ### (a) random intercepts > ### > # > n <- 12 > sss <- 20 > ttt <- 5 > scenario1 <- "sigma = 20, tau = 5" > # > ss1 <- diag( rep(sss^2,n) ) + ttt^2 > S1 <- matrix( 0, 2*n, 2*n ) > S1[1:n,1:n] <- ss1 > S1[(n+1):(2*n),(n+1):(2*n)] <- ss1 > # > sigma <- sqrt( mean( diag(ss1) ) ) > Ainv <- solve( t(X.a)%*%X.a ) > Ia1.naive <- sigma^2 * Ainv > Ia1.correct <- Ainv %*% ( t(X.a)%*%S1%*%X.a ) %*% Ainv > Ainv <- solve( t(X.b)%*%X.b ) > Ib1.naive <- sigma^2 * Ainv > Ib1.correct <- Ainv %*% ( t(X.b)%*%S1%*%X.b ) %*% Ainv > # > out1 <- cbind( + sqrt(diag(Ia1.naive)), sqrt(diag(Ia1.correct)), + sqrt(diag(Ib1.naive)), sqrt(diag(Ib1.correct)) ) > # > n <- 12 > sss <- 5 > ttt <- 20 > scenario2 <- "sigma = 5, tau = 20" > # > ss2 <- diag( rep(sss^2,n) ) + ttt^2 > S2 <- matrix( 0, 2*n, 2*n ) > S2[1:n,1:n] <- ss2 > S2[(n+1):(2*n),(n+1):(2*n)] <- ss2 > # > sigma <- sqrt( mean( diag(ss2) ) ) > S2 <- matrix( 0, 2*n, 2*n ) > S2[1:n,1:n] <- ss2 > S2[(n+1):(2*n),(n+1):(2*n)] <- ss2 > Ainv <- solve( t(X.a)%*%X.a ) > Ia2.naive <- sigma^2 * Ainv > Ia2.correct <- Ainv %*% ( t(X.a)%*%S2%*%X.a ) %*% Ainv > Ainv <- solve( t(X.b)%*%X.b ) > Ib2.naive <- sigma^2 * Ainv > Ib2.correct <- Ainv %*% ( t(X.b)%*%S2%*%X.b ) %*% Ainv > # > out2 <- cbind( + sqrt(diag(Ia2.naive)), sqrt(diag(Ia2.correct)), + sqrt(diag(Ib2.naive)), sqrt(diag(Ib2.correct)) ) > # > cat("\n random intercepts \n") random intercepts > cat( scenario1 ) sigma = 20, tau = 5> cat("\n") > dimnames(out1) <- list( c("beta0","beta1","beta2"), + c( paste( "A:", c("naive","sandwich") ), + paste( "B:", c("naive","sandwich") ) ) ) > print( round(out1/10,3) ) A: naive A: sandwich B: naive B: sandwich beta0 0.991 1.084 0.991 1.024 beta1 0.122 0.118 0.122 0.118 beta2 0.842 1.080 0.842 0.816 > cat( scenario2 ) sigma = 5, tau = 20> cat("\n") > dimnames(out2) <- list( c("beta0","beta1","beta2"), + c( paste( "A:", c("naive","sandwich") ), + paste( "B:", c("naive","sandwich") ) ) ) > print( round(out2/10,3) ) A: naive A: sandwich B: naive B: sandwich beta0 0.991 2.014 0.991 1.434 beta1 0.122 0.030 0.122 0.030 beta2 0.842 2.836 0.842 0.204 > # > # > ### > ### (b) random intercepts and slopes > ### > # > n <- 12 > D <- diag( c(100,2) ) > scenario1 <- "D = diag( 100, 2 ), sigma = 1" > ss1 <- X.a.0[,1:2]%*%D%*%t(X.a.0[,1:2]) + diag(rep(1,n)) > S1 <- matrix( 0, 2*n, 2*n ) > S1[1:n,1:n] <- ss1 > S1[(n+1):(2*n),(n+1):(2*n)] <- ss1 > # > sigma <- sqrt( mean( diag(ss1) ) ) > Ainv <- solve( t(X.a)%*%X.a ) > Ia1.naive <- sigma^2 * Ainv > Ia1.correct <- Ainv %*% ( t(X.a)%*%S1%*%X.a ) %*% Ainv > Ainv <- solve( t(X.b)%*%X.b ) > Ib1.naive <- sigma^2 * Ainv > Ib1.correct <- Ainv %*% ( t(X.b)%*%S1%*%X.b ) %*% Ainv > # > out1 <- cbind( + sqrt(diag(Ia1.naive)), sqrt(diag(Ia1.correct)), + sqrt(diag(Ib1.naive)), sqrt(diag(Ib1.correct)) ) > # > D <- diag( c(50,4) ) > scenario2 <- "D = diag( 50, 4 ), sigma = 1" > ss2 <- X.a.0[,1:2]%*%D%*%t(X.a.0[,1:2]) + diag(rep(1,n)) > S2 <- matrix( 0, 2*n, 2*n ) > S2[1:n,1:n] <- ss2 > S2[(n+1):(2*n),(n+1):(2*n)] <- ss2 > sigma <- sqrt( mean( diag(ss2) ) ) > Ainv <- solve( t(X.a)%*%X.a ) > Ia2.naive <- sigma^2 * Ainv > Ia2.correct <- Ainv %*% ( t(X.a)%*%S2%*%X.a ) %*% Ainv > Ainv <- solve( t(X.b)%*%X.b ) > Ib2.naive <- sigma^2 * Ainv > Ib2.correct <- Ainv %*% ( t(X.b)%*%S2%*%X.b ) %*% Ainv > # > out2 <- cbind( + sqrt(diag(Ia2.naive)), sqrt(diag(Ia2.correct)), + sqrt(diag(Ib2.naive)), sqrt(diag(Ib2.correct)) ) > # > cat("\n random intercepts and slopes\n") random intercepts and slopes > cat( scenario1 ) D = diag( 100, 2 ), sigma = 1> cat("\n") > dimnames(out1) <- list( c("beta0","beta1","beta2"), + c( paste( "A:", c("naive","sandwich") ), + paste( "B:", c("naive","sandwich") ) ) ) > print( round(out1/10,3) ) A: naive A: sandwich B: naive B: sandwich beta0 0.695 1.194 0.695 0.770 beta1 0.086 0.100 0.086 0.100 beta2 0.591 1.921 0.591 0.601 > cat( scenario2 ) D = diag( 50, 4 ), sigma = 1> cat("\n") > dimnames(out2) <- list( c("beta0","beta1","beta2"), + c( paste( "A:", c("naive","sandwich") ), + paste( "B:", c("naive","sandwich") ) ) ) > print( round(out2/10,3) ) A: naive A: sandwich B: naive B: sandwich beta0 0.786 1.161 0.786 0.658 beta1 0.097 0.142 0.097 0.142 beta2 0.668 2.093 0.668 0.850 > # > # > ### > ### (c) autocorrelated errors > ### > # > create.AR <- function( rho, n ){ + d1 <- rep( c(1:n), n ) + d2 <- rep( c(1:n), rep(n,n) ) + ccc <- rho^( abs(d1-d2) ) + matrix( ccc, n, n ) + } > # > sss <- 10 > ttt <- 20 > rho <- 0.9 > scenario1 <- "sigma = 10, tau = 20, rho = 0.9" > ss1 <- ttt^2 * create.AR( rho, 12 ) + diag( rep(sss^2,n) ) > S1 <- matrix( 0, 2*n, 2*n ) > S1[1:n,1:n] <- ss1 > S1[(n+1):(2*n),(n+1):(2*n)] <- ss1 > # > sigma <- sqrt( mean( diag( ss1 ) ) ) > Ainv <- solve( t(X.a)%*%X.a ) > Ia1.naive <- sigma^2 * Ainv > Ia1.correct <- Ainv %*% ( t(X.a)%*%S1%*%X.a ) %*% Ainv > Ainv <- solve( t(X.b)%*%X.b ) > Ib1.naive <- sigma^2 * Ainv > Ib1.correct <- Ainv %*% ( t(X.b)%*%S1%*%X.b ) %*% Ainv > # > out1 <- cbind( + sqrt(diag(Ia1.naive)), sqrt(diag(Ia1.correct)), + sqrt(diag(Ib1.naive)), sqrt(diag(Ib1.correct)) ) > # > sss <- 10 > ttt <- 20 > rho <- 0.5 > scenario2 <- "sigma = 10, tau = 20, rho = 0.5" > ss2 <- ttt^2 * create.AR( rho, 12 ) + diag( rep(sss^2,n) ) > S2 <- matrix( 0, 2*n, 2*n ) > S2[1:n,1:n] <- ss1 > S2[(n+1):(2*n),(n+1):(2*n)] <- ss2 > sigma <- sqrt( mean( diag( ss2 ) ) ) > Ainv <- solve( t(X.a)%*%X.a ) > Ia2.naive <- sigma^2 * Ainv > Ia2.correct <- Ainv %*% ( t(X.a)%*%S2%*%X.a ) %*% Ainv > Ainv <- solve( t(X.b)%*%X.b ) > Ib2.naive <- sigma^2 * Ainv > Ib2.correct <- Ainv %*% ( t(X.b)%*%S2%*%X.b ) %*% Ainv > # > out2 <- cbind( + sqrt(diag(Ia2.naive)), sqrt(diag(Ia2.correct)), + sqrt(diag(Ib2.naive)), sqrt(diag(Ib2.correct)) ) > # > cat("\n autocorrelation \n") autocorrelation > cat( scenario1 ) sigma = 10, tau = 20, rho = 0.9> cat("\n") > dimnames(out1) <- list( c("beta0","beta1","beta2"), + c( paste( "A:", c("naive","sandwich") ), + paste( "B:", c("naive","sandwich") ) ) ) > print( round(out1/10,3) ) A: naive A: sandwich B: naive B: sandwich beta0 1.075 2.021 1.075 1.728 beta1 0.132 0.172 0.132 0.172 beta2 0.913 2.379 0.913 1.125 > cat( scenario2 ) sigma = 10, tau = 20, rho = 0.5> cat("\n") > dimnames(out2) <- list( c("beta0","beta1","beta2"), + c( paste( "A:", c("naive","sandwich") ), + paste( "B:", c("naive","sandwich") ) ) ) > print( round(out2/10,3) ) A: naive A: sandwich B: naive B: sandwich beta0 1.075 2.034 1.075 1.597 beta1 0.132 0.176 0.132 0.176 beta2 0.913 1.950 0.913 1.179 > # > # > # end-of-file... > > > >