I am trying to replicate a semi-parametric analysis described in Harezlak,
Jaroslaw, David Ruppert, and Matt P. Wand. Semiparametric regression with R. New
York, NY: Springer, 2018.
(https://link.springer.com/book/10.1007%2F978-1-4939-8853-2).
I can successfully run the analysis, but now I'm trying to move it into my
workflow, which requires that the analysis be conducted within a function (using
targets package), and the `groupedData` function now fails with an error that it
cannot find the `X` matrix (see reprex below). I've tried the reprex on
both my personal Mac (where it works??) and on windows machines (where it does
not) - so the problem is likely specific to Windows computers (yes, this seems
weird to me too).
All packages have been updated, and I'm running the latest version of R on
all machines.
Reprex:
library(HRW) # contains example data and ZOSull function
library(nlme)
data(growthIndiana)
analyze_this <- function(df) {
mean.x <- mean(df$age)
mean.y <- mean(df$height)
sd.x <- sd(df$age)
sd.y <- sd(df$height)
df$x <- (df$age - mean.x) / sd.x
df$y <- (df$height - mean.y) / sd.y
X <- model.matrix(~ x * male * black, data = df)
dummyID <- rep(1, length(nrow(X)))
grouped_data <- groupedData(y ~ X[,-1]|rep(1, length = nrow(X)), data =
data.frame(y = df$y, X, dummyID))
}
# doesn't work on Windows machine, does work on the Mac
analyze_this(growthIndiana)
#> Error in eval(aux[[2]], object): object 'X' not found
# does work
df <- growthIndiana
mean.x <- mean(df$age)
mean.y <- mean(df$height)
sd.x <- sd(df$age)
sd.y <- sd(df$height)
df$x <- (df$age - mean.x) / sd.x
df$y <- (df$height - mean.y) / sd.y
X <- model.matrix(~ x * male * black, data = df)
dummyID <- rep(1, length(nrow(X)))
grouped_data <- groupedData(y ~ X[,-1]|rep(1, length = nrow(X)), data =
data.frame(y = df$y, X, dummyID))
# attempted work-around.
analyze_this2 <- function(df) {
num.global.knots = 20
num.subject.knots = 10
mean.x <- mean(df$age)
mean.y <- mean(df$height)
sd.x <- sd(df$age)
sd.y <- sd(df$height)
df$x <- (df$age - mean.x) / sd.x
df$y <- (df$height - mean.y) / sd.y
X <- model.matrix(~ x * male * black, data = df)
dummyID <- rep(1, length(nrow(X)))
# grouped_data <- groupedData(y ~ X[,-1]|rep(1, length = nrow(X)), data =
data.frame(y = df$y, X, dummyID))
global.knots = quantile(unique(df$x), seq(0, 1, length = num.global.knots +
2)[-c(1, num.global.knots + 2)])
subject.knots = quantile(unique(df$x), seq(0, 1, length = num.subject.knots +
2)[-c(1, num.subject.knots + 2)])
Z.global <- ZOSull(df$x, range.x = range(df$x), global.knots)
Z.group <- df$black * Z.global
Z.subject <- ZOSull(df$x, range.x = range(df$x), subject.knots)
Zblock <- list(
dummyID = pdIdent(~ 0 + Z.global),
dummyID = pdIdent(~ 0 + Z.group),
idnum = pdSymm(~ x),
idnum = pdIdent(~ 0 + Z.subject)
)
df$dummyID <- dummyID
tmp_data <- c(
df,
X = list(X),
Z.global = list(Z.global),
Z.group = list(Z.global),
Z.subject = list(Z.subject)
)
fit <- lme(y ~ 0 + X,
data = tmp_data,
random = Zblock
)
}
# this works (warning - lme takes awhile to fit)
analyze_this2(growthIndiana)
sessionInfo()
#> R version 4.1.2 (2021-11-01)
#> Platform: x86_64-w64-mingw32/x64 (64-bit)
#> Running under: Windows 10 x64 (build 22000)
#>
#> Matrix products: default
#>
#> locale:
#> [1] LC_COLLATE=English_United States.1252
#> [2] LC_CTYPE=English_United States.1252
#> [3] LC_MONETARY=English_United States.1252
#> [4] LC_NUMERIC=C
#> [5] LC_TIME=English_United States.1252
#>
#> attached base packages:
#> [1] stats graphics grDevices utils datasets methods base
#>
#> other attached packages:
#> [1] nlme_3.1-153 HRW_1.0-5
#>
#> loaded via a namespace (and not attached):
#> [1] lattice_0.20-45 digest_0.6.29 withr_2.4.3 grid_4.1.2
#> [5] magrittr_2.0.1 reprex_2.0.1 evaluate_0.14
KernSmooth_2.23-20
#> [9] highr_0.9 stringi_1.7.6 rlang_0.4.12 cli_3.1.0
#> [13] rstudioapi_0.13 fs_1.5.2 rmarkdown_2.11
splines_4.1.2
#> [17] tools_4.1.2 stringr_1.4.0 glue_1.6.0 xfun_0.29
#> [21] yaml_2.2.1 fastmap_1.1.0 compiler_4.1.2
htmltools_0.5.2
#> [25] knitr_1.37
Created on 2022-01-07 by the [reprex package](https://reprex.tidyverse.org)
(v2.0.1)
Here's the session Info for the Mac:
sessionInfo()
R version 4.1.2 (2021-11-01)
Platform: aarch64-apple-darwin20 (64-bit)
Running under: macOS Monterey 12.1
Matrix products: default
LAPACK:
/Library/Frameworks/R.framework/Versions/4.1-arm64/Resources/lib/libRlapack.dylib
locale:[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] nlme_3.1-153 targets_0.8.1 HRW_1.0-5
loaded via a namespace (and not attached):
[1] compiler_4.1.2 pillar_1.6.4 tools_4.1.2 digest_0.6.28
lattice_0.20-45 jsonlite_1.7.2 evaluate_0.14
[8] lifecycle_1.0.1 tibble_3.1.6 pkgconfig_2.0.3 rlang_0.4.12
igraph_1.2.9 cli_3.1.0 yaml_2.2.1
[15] xfun_0.28 fastmap_1.1.0 withr_2.4.2 knitr_1.36
htmlwidgets_1.5.4 vctrs_0.3.8 grid_4.1.2
[22] tidyselect_1.1.1 glue_1.5.0 data.table_1.14.2 R6_2.5.1
processx_3.5.2 fansi_0.5.0 rmarkdown_2.11
[29] callr_3.7.0 purrr_0.3.4 magrittr_2.0.1 scales_1.1.1
codetools_0.2-18 ps_1.6.0 htmltools_0.5.2
[36] ellipsis_0.3.2 splines_4.1.2 colorspace_2.0-2 KernSmooth_2.23-20
utf8_1.2.2 visNetwork_2.1.0 munsell_0.5.0
[43] crayon_1.4.2
Melissa Key, Ph.D.
Statistician
Non-Invasive Brain Stimulation Team
Infoscitex
Cell: 937-550-4981
Email: mkey at infoscitex.com<mailto:mkey at infoscitex.com>
Base email: melissa.key.ctr at us.af.mil<mailto:melissa.key.ctr at
us.af.mil>
[[alternative HTML version deleted]]