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