################################################### set.seed(1231) library("basicspace") N <- 1000 M <- 20 s <- 2 fraction.missing <- 0.3 E <- matrix(runif(N * M, min = -0.5, max = 0.5), nrow = N, ncol = M) ################################################### U <- matrix(runif(N * s), nrow = N, ncol = s) D <- diag(seq(from = 2.1, by = -0.2, length.out = s)) V.prime <- matrix(runif(s * M), nrow = s, ncol = M) c <- rnorm(M) Jn <- rep(1, N) ################################################### X.true <- U %*% D %*% V.prime + Jn %o% c X.0 <- X.true + E Psi.true <- U %*% sqrt(D) W.true <- t(V.prime) %*% sqrt(D) ################################################### missing <- sample(1:(N * M), round(fraction.missing * N * M)) X.0[missing] <- 999 ################################################### rownames(X.0) <- paste("Legis", 1:N, sep = "") colnames(X.0) <- paste("V", 1:M, sep = "") ################################################### result <- blackbox(X.0, missing = c(999), verbose = TRUE, dims = 3, minscale = 8) names(result)[1:4] ################################################### Psi.hat <- cbind(result$individuals[[2]]$c1, result$individuals[[2]]$c2) c.hat <- result$stimuli[[2]]$c xrow <- sapply(1:N, function(i) length(rep(1, s)[!is.na(Psi.hat[i, ])])) Psi.hat <- Psi.hat[!(xrow < 2), ] Psi.true <- Psi.true[!(xrow < 2), ] Psi.hat[, 1] <- Psi.hat[, 1] - mean(Psi.hat[, 1]) Psi.hat[, 2] <- Psi.hat[, 2] - mean(Psi.hat[, 2]) Psi.true[, 1] <- Psi.true[, 1] - mean(Psi.true[, 1]) Psi.true[, 2] <- Psi.true[, 2] - mean(Psi.true[, 2]) C <- t(Psi.true) %*% Psi.hat svddecomp <- svd(C) U.rotate <- svddecomp$u V.rotate <- svddecomp$v T <- V.rotate %*% t(U.rotate) Psi.hatrotate <- Psi.hat %*% T par(mfrow = c(1, 2)) plot(Psi.true[, 1], Psi.hatrotate[, 1], xlim = c(-0.7, 0.7), ylim = c(-0.5, 0.5), pch = 20, cex = 0.4, cex.lab = 1.6, bty = "n", xlab = "True Psi, first dimension", ylab = "Recovered Psi, first dimension") plot(Psi.true[, 2], Psi.hatrotate[, 2], xlim = c(-0.7, 0.7), ylim = c(-0.5, 0.5), pch = 20, cex = 0.4, cex.lab = 1.6, bty = "n", xlab = "True Psi, second dimension", ylab = "Recovered Psi, second dimension") ################################################### W.hat <- cbind(result$stimuli[[2]]$w1, result$stimuli[[2]]$w2) W.hatrotate <- W.hat %*% T par(mfrow = c(1, 2)) plot(W.true[, 1], W.hatrotate[, 1], xlim = c(0, 1.5), ylim = c(-0.75, 2.75), pch = 20, cex = 1.5, cex.lab = 1.6, bty = "n", xlab = "True W, first dimension", ylab = "Recovered W, first dimension") plot(W.true[, 2], W.hatrotate[, 2], xlim = c(0, 1.5), ylim = c(-0.75, 2.75), pch = 20, cex = 1.5, cex.lab = 1.6, bty = "n", xlab = "True W, second dimension", ylab = "Recovered W, second dimension") ################################################### par(mfrow = c(1, 1)) plot(c, c.hat, pch = 20, cex = 1.2, cex.lab = 1.1, bty = "n", xlab = "True C", ylab = "Recovered C") ################################################### W.hat <- cbind(result$stimuli[[2]]$w1, result$stimuli[[2]]$w2) Psi.hat <- cbind(result$individuals[[2]]$c1, result$individuals[[2]]$c2) X.hat <- Psi.hat %*% t(W.hat) + Jn %o% result$stimuli[[2]]$c par(mfrow = c(1, 2)) plot(X.true[missing], X.hat[missing], pch = 20, cex = 0.4, cex.lab = 1.2, bty = "n", xlab = "True X, missing values", ylab = "Recovered X, missing values") plot(X.true[!(1:(N * M) %in% missing)], X.hat[!(1:(N * M) %in% missing)], pch = 20, cex = 0.4, cex.lab = 1.2, bty = "n", xlab = "True X, nonmissing values", ylab = "Recovered X, nonmissing values") ################################################### data("Issues1980") Issues1980[1:10, 1:4] ################################################### Issues1980[Issues1980[, "abortion1"] == 7, "abortion1"] <- 8 Issues1980[Issues1980[, "abortion2"] == 7, "abortion2"] <- 8 ################################################### Issues1980_bb <- blackbox(Issues1980, missing = c(0, 8, 9), verbose = FALSE, dims = 3, minscale = 8) ################################################### summary(Issues1980_bb) ################################################### cor(Issues1980_bb$individuals[[1]]$c1, Issues1980[, "libcon1"], use = "pairwise") ################################################### data("LC1980") LCdat <- LC1980[, -1] LCdat[1:10, ] ## Commented to shorten runtimes LC1980_bbt <- blackbox_transpose(LCdat, ## missing=c(0,8,9), dims=3, minscale=5,verbose=TRUE) data(LC1980_bbt) ################################################### par(mfrow = c(1, 2)) plot(LC1980_bbt) plotcdf.blackbt(LC1980_bbt) ################################################### summary(LC1980_bbt) ################################################### data("LC1980") result <- aldmck(data = LC1980, polarity = 2, respondent = 1, missing = c(0, 8, 9), verbose = TRUE) summary(result) ################################################### plot.aldmck(result) ################################################### result <- boot_aldmck(data = LC1980, polarity = 2, respondent = 1, missing = c(0, 8, 9), iter = 100) apply(result, 2, sd) ################################################### Nstimuli <- 6 Nresp <- 500 Z_j <- rnorm(6) Z_j <- (Z_j - mean(Z_j))/sd(Z_j) respondent.sd <- runif(Nresp, min = 0.3, max = 0.9) error_heteroskedastic <- matrix(NA, Nresp, Nstimuli) for (i in 1:Nresp) error_heteroskedastic <- rnorm(Nstimuli, sd = respondent.sd) w_i <- runif(Nresp, min = 0, max = 1) c_i <- rnorm(Nresp) Y_ij <- rep(1, 500) %o% Z_j Y_ij <- Y_ij + error_heteroskedastic R_ij <- 1/w_i %o% rep(1, Nstimuli) * (Y_ij - c_i %o% rep(1, Nstimuli)) result <- aldmck(R_ij, polarity = 6, missing = c(999)) cor(Z_j, result$stimuli)