1###{{{ Objective 2 3IV_method.lvm <- NULL 4IV_objective.lvm <- function(x,p,data,...) { 5 IV2(x,data,...) 6} 7IV_variance.lvm <- function(x,p,data,opt,...) { 8 opt$vcov 9} 10 11IV0_method.lvm <- NULL 12IV0_objective.lvm <- function(x,p,data,...) { 13 IV2(x,data,type="non-robust",...) 14} 15IV0_variance.lvm <- function(x,p,data,opt,...) { 16 opt$vcov 17} 18 19IV1_method.lvm <- NULL 20IV1_objective.lvm <- function(x,p,data,...) { 21 IV(x,data) 22} 23IV1_variance.lvm <- function(x,p,data,opt,...) { 24 opt$vcov 25} 26 27###}}} Objective 28 29CondVar <- function(S,idx) { 30 idx2 <- setdiff(seq_len(ncol(S)),idx) 31 S11 <- S[idx2,idx2]; 32 S22 <- S[idx,idx] 33 S12 <- S[idx2,idx] 34 S11-S12%*%solve(S22)%*%t(S12) 35} 36 37varest <- function(x,data) { 38 p <- IV(x,data)$estimate 39 idx <- match(names(p),coef(x,mean=TRUE)) 40 x0 <- parfix(Model(x),idx,p) 41 index(x0) <- reindex(x0,zeroones=TRUE,deriv=TRUE) 42 43 A <- t(index(x)$A) 44 Afix <- A; Afix[t(index(x)$M0)==1] <- 0 45 A[A!=0] <- 1 46 k <- nrow(A) 47 I <- diag(nrow=k) 48 Ap <- modelVar(x)$A ## Estimated parameter matrix 49 50 indicators <- setdiff(vars(x)[rowSums(A)==1],exogenous(x)) 51 responses <- endogenous(x,top=TRUE) 52 y.indicators <- responses[rowSums(A[responses,])==1] 53 Sigma <- var(data[,manifest(x)]) 54 55 var.eta <- c() 56 for (eta in latent(x)) { 57 reachable <- acc(x$M,eta) 58 ys <- intersect(names(reachable),y.indicators) 59 lambdas <- c() 60 for (y in ys) { 61 pp <- path(Model(x), from=eta, to=y) 62 lambda1 <- 0 63 for (i in seq_along(pp)) { 64 lambda <- 1 65 for (j in seq_len(length(pp[[i]])-1)) 66 lambda <- lambda*Ap[pp[[i]][j],pp[[i]][j+1]] 67 lambda1 <- lambda1+lambda 68 } 69 lambdas <- c(lambdas,lambda1) 70 } 71 val <- outer(1/lambdas,1/lambdas)*Sigma[ys,ys] 72 var.eta <- c(var.eta, mean(val[upper.tri(val)])) 73 } 74 75 S <- rep(0,k); S[match(manifest(x),vars(x))] <- diag(Sigma); S[match(latent(x),vars(x))] <- var.eta; names(S) <- vars(x) 76 I <- diag(nrow=k) 77 IA <- (I-t(Ap)) 78 IA%*%cbind(S)%*%t(IA) 79 80} 81 82 83## Instrumental Variable Estimator / 2SLS 84##' @export 85IV <- function(m,data,R2thres=0,type="robust", ...) { 86 if (length(constrain(m))>0) stop("Nonlinear constrains not supported!") 87 if (inherits(m,"lvmfit")) { 88 m <- Model(m) 89 } 90 R2 <- cor(data[,manifest(m)])^2 91 92 A <- t(index(m)$A) 93 Afix <- A; Afix[t(index(m)$M0)==1] <- 0 94 A[A!=0] <- 1 95 P <- index(m)$P 96 k <- nrow(A) 97 I <- diag(nrow=k) 98 B <- rbind(I,solve(I-A)) 99 VV <- B%*%P%*%t(B) 100 u.var <- index(m)$vars 101 lat.idx <- with(index(m), which(vars%in%latent)) 102 if (length(lat.idx)==0) stop("Estimator only defined for models with latent variable") 103 y.var <- endogenous(m) 104 y.idx <- which(index(m)$vars%in%y.var) 105 x.idx <- which(vars(m)%in%exogenous(m)) 106 107 ## Set of Indicator variables: 108 indicators <- c() 109 for (i in seq_len(nrow(A))) { 110 ifix <- (Afix[i,]==1) 111 if ((sum(ifix)==1) & all(A[i,-which(ifix)]==0)) 112 indicators <- c(indicators, i) 113 } 114 y.indicators <- intersect(indicators,y.idx) 115 116 y.scale <- list() 117 for (eta in lat.idx) { 118 pred.eta <- intersect(y.idx, which(Afix[,eta]==1)) ## Candidates for 119 ## eta = y-epsilon 120 if (length(pred.eta)<1) 121 pred.eta <- intersect(lat.idx, which(Afix[,eta]==1)) 122 myidx <- c() 123 for (y in pred.eta) { 124 y.pred <- setdiff(eta,which(A[y,]==1)) ## No other variables predicting y? 125 if (length(y.pred)==0) 126 myidx <- c(myidx,y) 127 } 128 y.scale <- c(y.scale, list(myidx)) 129 } 130 131 if (any(unlist(lapply(y.scale, function(x) length(x)))<1)) stop("At least one scale-measurement pr. latent variable") 132 133 vv <- setdiff(seq_len(k),c(unlist(y.scale),x.idx)) 134 135 Ecor <- list() 136 eta.surrogate <- c() 137 latno <- 0 138 for (e in lat.idx) { 139 latno <- latno+1 140 y0 <- y.scale[[latno]][1] 141 if (!(y0%in%lat.idx)) { 142 eta.surrogate <- c(eta.surrogate,vars(m)[y0]) 143 Ecor <- c(Ecor,list(y0)) 144 } 145 else { 146 v0 <- vars(m)[-c(e,indicators)] 147 ##m.sub <- subset(m,vars(m)[c(e,indicators)]) 148 m.sub <- rmvar(m,v0) 149 i <- 0 150 while (i<length(y.indicators)) { 151 i <- i+1 152 pp <- path(m.sub,from=vars(m)[e],to=vars(m)[y.indicators[i]])[[1]] 153 if (!is.null(pp)) { 154 Ecor <- c(Ecor, 155 list(which(vars(m)%in%pp[-1]))) 156 eta.surrogate <- c(eta.surrogate, tail(pp,1)) 157 } 158 } 159 } 160 }; 161 names(eta.surrogate) <- latent(m) 162 163 coefname <- coef(m,mean=TRUE) 164 P0 <- P 165 V. <- list() 166 Z. <- list() 167 Y. <- list() 168 count <- 0 169 ff <- list() 170 instruments <- c() 171 parname <- c() 172 for (v in vv) { 173 pred <- which(A[v,]==1) 174 if (sum(pred)>0) { 175 Debug(vars(m)[v]) 176 pred.lat <- intersect(pred,lat.idx) # Any latent predictors? 177 lpos <- match(v,lat.idx) 178 lppos <- match(pred.lat,lat.idx) 179 ecor <- c(v,unlist(Ecor[lppos])) 180 if (!is.na(lpos)) { 181 v0 <- match(eta.surrogate[lpos],vars(m)) 182 ecor <- c(ecor,Ecor[[lpos]]) 183 } else { 184 v0 <- v 185 } 186 187 ecor <- unique(c(v0,ecor)) 188 XX <- vars(m)[A[v,]==1] 189 intpred <- exogenous(m) 190 newf <- c() 191 if (length(pred.lat)>0) { 192 intpred <- vars(m) 193 for (i in seq_along(pred.lat)) { 194 uncor <- which(colSums(VV[ecor,k+seq_len(k),drop=FALSE])==0) 195 uncor <- setdiff(uncor,c(lat.idx)) 196 mypred <- vars(m)[uncor] 197 XX[XX==vars(m)[pred.lat[i]]] <- eta.surrogate[lppos[i]] 198 intpred <- intersect(intpred,mypred) 199 f <- toformula(eta.surrogate[lppos[i]],mypred) 200 ff <- c(ff, 201 f) 202 f2 <- list(f) 203 names(f2) <- vars(m)[i] 204 newf <- c(newf,f2) 205 } 206 } 207 208 intpred <- intersect(intpred,manifest(m)) 209 R2max <- apply(R2[XX,intpred,drop=FALSE],2,max) 210 if (any(R2max<R2thres)) intpred <- intpred[R2max>=R2thres] 211 newf <- list(intpred); names(newf) <- vars(m)[v] 212 instruments <- c(instruments, newf) 213 covariates <- unique(c(setdiff(colnames(A)[A[v,]==1],latent(m)),intpred)) 214 if (length(covariates)==0) stop("No instruments") 215 Z <- model.matrix(toformula("",c("1",XX)),data) 216 Y <- as.matrix(data[,vars(m)[v0]]) 217 V <- model.matrix(toformula("",c("1",unique(covariates))),data) 218 count <- count+1 219 V. <- c(V.,list(V)) 220 Z. <- c(Z.,list(Z)) 221 Y. <- c(Y.,list(Y)) 222 XX <- vars(m)[A[v,]==1 & Afix[v,]!=1] 223 parname <- c(parname, c(vars(m)[v0],paste(vars(m)[v],XX,sep=lava.options()$symbol[1]))) 224 } else { 225 if (vars(m)[v]%in%latent(m)) { 226 lpos <- match(v,lat.idx) 227 v0 <- match(eta.surrogate[lpos],vars(m)) 228 Y <- matrix(data[,vars(m)[v0]],ncol=1) 229 Y. <- c(Y.,list(Y)) 230 V. <- c(V.,list(cbind(rep(1,nrow(Y))))) 231 Z. <- c(Z.,list(cbind(rep(1,nrow(Y))))) 232 parname <- c(parname, names(eta.surrogate)[lpos]) 233 } 234 } 235 } 236 237 LS <- function(X) { 238 with(svd(X), v%*%diag(1/d,nrow=length(d))%*%t(u)) 239 } 240 ##projection <- function(X) X%*%LS(X) 241 P0 <- lapply(V.,LS) 242 Zhat <- list(); for (i in seq_along(Z.)) Zhat <- c(Zhat, list(V.[[i]]%*%(P0[[i]]%*%Z.[[i]]))) 243 ZhatLS <- lapply(Zhat,LS) 244 theta <- list(); for (i in seq_along(Y.)) theta <- c(theta, list(ZhatLS[[i]]%*%Y.[[i]])) 245 u <- c() 246 for (i in seq_along(Y.)) 247 u <- cbind(u, Y.[[i]]-Z.[[i]]%*%theta[[i]]) 248 249 covu <- crossprod(u)/nrow(u) 250 251 theta.npar <- unlist(lapply(theta,length)) 252 theta.ncum <- c(0,cumsum(theta.npar)) 253 vartheta <- matrix(0,ncol=sum(theta.npar),nrow=sum(theta.npar)) 254 for (i in seq_along(theta)) { 255 for (j in seq(i,length(theta))) { 256 idx1 <- seq_len(theta.npar[i]) + theta.ncum[i] 257 idx2 <- seq_len(theta.npar[j]) + theta.ncum[j] 258 if (type=="robust") { 259 zi <- ZhatLS[[i]] 260 for (k in seq(nrow(zi))) zi[k,] <- zi[k,]*u[,i] 261 zj <- ZhatLS[[j]] 262 for (k in seq(nrow(zj))) zj[k,] <- zj[k,]*u[,j] 263 uZZ <- zi%*%t(zj) 264 } else { 265 uZZ <- covu[i,j]* (ZhatLS[[i]]%*%t(ZhatLS[[j]])) 266 } 267 vartheta[idx1,idx2] <- uZZ 268 if (i!=j) { 269 vartheta[idx2,idx1] <- t(uZZ) 270 } 271 } 272 } 273 274 parname[which(parname%in%eta.surrogate)] <- names(eta.surrogate)[which(eta.surrogate%in%parname)] 275 276 coef <- cbind(unlist(theta),diag(vartheta)^0.5); rownames(coef) <- parname; colnames(coef) <- c("Estimate","Std.Err") 277 res <- list(estimate=coef[,1], vcov=vartheta) 278 attributes(res)$surrogates <- eta.surrogate 279 attributes(res)$instruments <- instruments 280 return(res) 281} 282 283IV2 <- function(m,data,control=list(),...) { 284 if (is.null(control$R2thres)) control$R2thres <- 0 285 res <- IV(m,data,R2thres=control$R2thres,...) 286 p <- res$estimate 287 idx <- match(names(p),coef(m,mean=TRUE)) 288 x0 <- parfix(m,idx,p) 289 index(x0) <- reindex(x0,zeroones=TRUE,deriv=TRUE) 290 idx0 <- order(idx) 291 p0 <- p[idx0] 292 V0 <- res$vcov[idx0,idx0] 293 if (is.null(control$variance) || control$variance) { 294 suppressWarnings(e0 <- estimate(x0,data,...,messages=0,quick=TRUE)) 295 p0 <- c(p0,e0) 296 V0 <- V0%++%matrix(0,ncol=length(e0),nrow=length(e0)) 297 } 298 R2 <- noquote(formatC(cor(data[,manifest(m)])^2)) 299 colnames(R2) <- rownames(R2) <- manifest(m) 300 l1 <- noquote(rbind(paste(latent(m),collapse=","), 301 paste(attributes(res)$surrogates,collapse=","), 302 "")) 303 rownames(l1) <- c("Latent variables","Surrogate variables:","") 304 colnames(l1) <- "" 305 ii <- attributes(res)$instruments 306 I <- noquote(matrix(NA,ncol=2,nrow=length(ii))) 307 rownames(I) <- rep("",nrow(I)) 308 colnames(I) <- c("Response","Instruments") 309 for (i in seq_along(ii)) { 310 I[i,] <- c(names(ii)[i],paste(ii[[i]],collapse=",")) 311 } 312 mymsg <- list(l1,I); 313 list(estimate=p0,vcov=V0,summary.message=function(...) { 314 mymsg }) 315} 316