【R】SvyPVpm()関数を用いて複雑な調査デザインによって作成された大規模データの処理をしたい

Question

前提・実現したいこと

・RでOECDによる国際成人力調査(PIAAC)の分析をしています
・PIAACは家庭環境や数的思考力(ソースコード中NUMERIC)・読解力(LITERACY)・問題解決力(PROBLEMSOLVING)を含む個票データです

・svyPVpackパッケージに含まれるsvyPVpm()関数を用いて、平均・標準誤差・母集団の合計・ウエイトの合計を計算したいです

・このページを参照して作成しています。参照ページ記載コードからの変更点は以下。残りはコピペです

res1 <- svyPVpm(by = ~ GENDER_R, svydat=sd2, pvs=LITERACY)

↓

res1 <- svyPVpm(by = ~ B_Q01a, svydat=sd2, pvs=PROBLEMSOLVING)

・B_Q01aは教育水準の指標を格納した列です

発生している問題・エラーメッセージ

・参照ページ

res1 <- svyPVpm(by = ~ GENDER_R, svydat=sd2, pvs=LITERACY)

↓

res1 <- svyPVpm(by = ~ B_Q01a, svydat=sd2, pvs=PROBLEMSOLVING)

・以上のようにコードを書き換えると、以下のエラーメッセージが表示されます(元コードは正常に動作します)

 svrVar(repmeans, scale, rscales, mse = design$mse, coef = rval) でエラー: 
  All replicates contained NAs

該当のソースコード

library(memisc)
library(survey)
library(svyPVpack)
library(ggplot2)

dir.create("firstEX")
download.file("https://webfs.oecd.org/piaac/puf-data/SPSS/prgautp1.sav","./firstEX/prgautp1.sav", method="auto", mode="wb")
setwd("./firstEX")

#1
d <- as.data.set(spss.system.file("prgautp1.sav"))
d1 <- data.frame(d)

#2
hmNA <- sapply(d1, function(x) sum(is.na(x)));
onlyNAcols <- which(hmNA == nrow(d1));
length(onlyNAcols);

d2 <- d1[,-onlyNAcols];

# get the columns which contain the weights
nr.SPFWT0    <- grep("^SPFWT0$",colnames(d2),value=F)
nr.SPFWT_rep <- grep("^SPFWT[1-9][0-9]*$",colnames(d2),value=F)

# get the variable names of the plausible values
LITERACY        <- grep("^PVLIT[0-9]+", colnames(d2),value=TRUE)
NUMERACY        <- grep("^PVNUM[0-9]+", colnames(d2),value=TRUE)
PROBLEMSOLVING  <- grep("^PVPSL[0-9]+", colnames(d2),value=TRUE)

# create the survey design
sd2 <- svrepdesign(variables=d2, repweights=d2[,nr.SPFWT_rep],weights=d2[,nr.SPFWT0], type="JK1")

res1 <- svyPVpm(by = ~ B_Q01a, svydat=sd2, pvs=PROBLEMSOLVING)

試したこと

・データ型の確認
GENDER_Rはnumeric
B_Q01aはcharactor

補足情報（FW/ツールのバージョンなど）

R version 3.6.0 (2019-04-26)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 17134)

Matrix products: default

attached base packages:
[1] grid stats graphics grDevices utils
[6] datasets methods base

other attached packages:
[1] ggplot2_3.2.1 svyPVpack_0.1-1 survey_3.36
[4] survival_2.44-1.1 Matrix_1.2-17 memisc_0.99.20.1
[7] MASS_7.3-51.4 lattice_0.20-38

loaded via a namespace (and not attached):
[1] Rcpp_1.0.1 rstudioapi_0.10 splines_3.6.0
[4] munsell_0.5.0 colorspace_1.4-1 rlang_0.4.0
[7] tools_3.6.0 gtable_0.3.0 DBI_1.0.0
[10] withr_2.1.2 htmltools_0.4.0 mitools_2.4
[13] lazyeval_0.2.2 digest_0.6.23 tibble_2.1.1
[16] crayon_1.3.4 repr_1.0.1 base64enc_0.1-3
[19] compiler_3.6.0 pillar_1.4.2 scales_1.0.0
[22] jsonlite_1.6 pkgconfig_2.0.3

Answer 1

"All replicates contained NAs"とエラーが出てますので、NAが含まれているのが原因かと思います。
一度NAが含まれているかどうか確認してみてください

Comments

[Koki328]

```R
d1$B_Q01a <- as.numeric(as.character(d1$B_Q01a))
```

と実行し、B_Q01a列をnumericとして扱うことで、一部の国のデータに対しては正しく出力することができました。

エラーが出ている関数(svrvar)の内容は以下の通りです

```R
function (thetas, scale, rscales, na.action = getOption("na.action"),
mse = getOption("survey.replicates.mse"), coef)
{
thetas <- get(na.action)(thetas)
naa <- attr(thetas, "na.action")
if (!is.null(naa)) {
rscales <- rscales[-naa]
if (length(rscales))
warning(length(naa), " replicates gave NA results and were discarded.")
else stop("All replicates contained NAs")
}
if (is.null(mse))
mse <- FALSE
if (length(dim(thetas)) == 2) {
if (mse) {
meantheta <- coef
}
else {
meantheta <- colMeans(thetas[rscales > 0, , drop = FALSE])
}
v <- crossprod(sweep(thetas, 2, meantheta, "-") *
sqrt(rscales)) * scale
}
else {
if (mse) {
meantheta <- coef
}
else {
meantheta <- mean(thetas[rscales > 0])
}
v <- sum((thetas - meantheta)^2 * rscales) * scale
}
attr(v, "na.replicates") <- naa
attr(v, "means") <- meantheta
return(v)
}
```

現在は、エラーメッセージ部分の理解を進めている状況です。