{"package":"smoothr","topic":"densify","snippet":"### Name: densify\n### Title: Densify spatial lines or polygons\n### Aliases: densify\n\n### ** Examples\n\nlibrary(sf)\nl <- jagged_lines$geometry[[2]]\nl_dense <- densify(l, n = 2)\nplot(l, lwd = 5)\nplot(l_dense, col = \"red\", lwd = 2, lty = 2, add = TRUE)\nplot(l_dense %>% st_cast(\"MULTIPOINT\"), col = \"red\", pch = 19,\n add = TRUE)\n\n\n"} {"package":"smoothr","topic":"drop_crumbs","snippet":"### Name: drop_crumbs\n### Title: Remove small polygons or line segments\n### Aliases: drop_crumbs\n\n### ** Examples\n\n# remove polygons smaller than 200km2\np <- jagged_polygons$geometry[7]\narea_thresh <- units::set_units(200, km^2)\np_dropped <- drop_crumbs(p, threshold = area_thresh)\n# plot\npar(mar = c(0, 0, 1, 0), mfrow = c(1, 2))\nplot(p, col = \"black\", main = \"Original\")\nif (length(p_dropped) > 0) {\n plot(p_dropped, col = \"black\", main = \"After drop_crumbs()\")\n}\n\n\n# remove lines less than 25 miles\nl <- jagged_lines$geometry[8]\n# note that any units can be used\n# conversion to units of projection happens automatically\nlength_thresh <- units::set_units(25, miles)\nl_dropped <- drop_crumbs(l, threshold = length_thresh)\n# plot\npar(mar = c(0, 0, 1, 0), mfrow = c(1, 2))\nplot(l, lwd = 5, main = \"Original\")\nif (length(l_dropped)) {\n plot(l_dropped, lwd = 5, main = \"After drop_crumbs()\")\n}\n\n\n\n"} {"package":"smoothr","topic":"fill_holes","snippet":"### Name: fill_holes\n### Title: Fill small holes in polygons\n### Aliases: fill_holes\n\n### ** Examples\n\n# fill holes smaller than 1000km2\np <- jagged_polygons$geometry[5]\narea_thresh <- units::set_units(1000, km^2)\np_dropped <- fill_holes(p, threshold = area_thresh)\n# plot\npar(mar = c(0, 0, 1, 0), mfrow = c(1, 2))\nplot(p, col = \"black\", main = \"Original\")\nplot(p_dropped, col = \"black\", main = \"After fill_holes()\")\n\n\n"} {"package":"smoothr","topic":"smooth","snippet":"### Name: smooth\n### Title: Smooth a spatial feature\n### Aliases: smooth\n\n### ** Examples\n\nlibrary(sf)\n# compare different smoothing methods\n# polygons\npar(mar = c(0, 0, 0, 0), oma = c(4, 0, 0, 0), mfrow = c(3, 3))\np_smooth_chaikin <- smooth(jagged_polygons, method = \"chaikin\")\np_smooth_ksmooth <- smooth(jagged_polygons, method = \"ksmooth\")\np_smooth_spline <- smooth(jagged_polygons, method = \"spline\")\nfor (i in 1:nrow(jagged_polygons)) {\n plot(st_geometry(p_smooth_spline[i, ]), col = NA, border = NA)\n plot(st_geometry(jagged_polygons[i, ]), col = \"grey40\", border = NA, add = TRUE)\n plot(st_geometry(p_smooth_chaikin[i, ]), col = NA, border = \"#E41A1C\", lwd = 2, add = TRUE)\n plot(st_geometry(p_smooth_ksmooth[i, ]), col = NA, border = \"#4DAF4A\", lwd = 2, add = TRUE)\n plot(st_geometry(p_smooth_spline[i, ]), col = NA, border = \"#377EB8\", lwd = 2, add = TRUE)\n}\npar(fig = c(0, 1, 0, 1), oma = c(0, 0, 0, 0), new = TRUE)\nplot(0, 0, type = \"n\", bty = \"n\", xaxt = \"n\", yaxt = \"n\", axes = FALSE)\nlegend(\"bottom\", legend = c(\"chaikin\", \"ksmooth\", \"spline\"),\n col = c(\"#E41A1C\", \"#4DAF4A\", \"#377EB8\"),\n lwd = 2, cex = 2, box.lwd = 0, inset = 0, horiz = TRUE)\n\n# lines\npar(mar = c(0, 0, 0, 0), oma = c(4, 0, 0, 0), mfrow = c(3, 3))\nl_smooth_chaikin <- smooth(jagged_lines, method = \"chaikin\")\nl_smooth_ksmooth <- smooth(jagged_lines, method = \"ksmooth\")\nl_smooth_spline <- smooth(jagged_lines, method = \"spline\")\nfor (i in 1:nrow(jagged_lines)) {\n plot(st_geometry(l_smooth_spline[i, ]), col = NA)\n plot(st_geometry(jagged_lines[i, ]), col = \"grey20\", lwd = 3, add = TRUE)\n plot(st_geometry(l_smooth_chaikin[i, ]), col = \"#E41A1C\", lwd = 2, lty = 2, add = TRUE)\n plot(st_geometry(l_smooth_ksmooth[i, ]), col = \"#4DAF4A\", lwd = 2, lty = 2, add = TRUE)\n plot(st_geometry(l_smooth_spline[i, ]), col = \"#377EB8\", lwd = 2, lty = 2, add = TRUE)\n}\npar(fig = c(0, 1, 0, 1), oma = c(0, 0, 0, 0), new = TRUE)\nplot(0, 0, type = \"n\", bty = \"n\", xaxt = \"n\", yaxt = \"n\", axes = FALSE)\nlegend(\"bottom\", legend = c(\"chaikin\", \"smooth\", \"spline\"),\n col = c(\"#E41A1C\", \"#4DAF4A\", \"#377EB8\"),\n lwd = 2, cex = 2, box.lwd = 0, inset = 0, horiz = TRUE)\n\n\n"} {"package":"smoothr","topic":"smooth_chaikin","snippet":"### Name: smooth_chaikin\n### Title: Chaikin's corner cutting algorithm\n### Aliases: smooth_chaikin\n\n### ** Examples\n\n# smooth_chaikin works on matrices of coordinates\n# use the matrix of coordinates defining a polygon as an example\nm <- jagged_polygons$geometry[[2]][[1]]\nm_smooth <- smooth_chaikin(m, wrap = TRUE)\nclass(m)\nclass(m_smooth)\nplot(m, type = \"l\", axes = FALSE, xlab = NA, ylab = NA)\nlines(m_smooth, col = \"red\")\n\n# smooth is a wrapper for smooth_chaikin that works on spatial features\nlibrary(sf)\np <- jagged_polygons$geometry[[2]]\np_smooth <- smooth(p, method = \"chaikin\")\nclass(p)\nclass(p_smooth)\nplot(p)\nplot(p_smooth, border = \"red\", add = TRUE)\n\n\n"} {"package":"smoothr","topic":"smooth_densify","snippet":"### Name: smooth_densify\n### Title: Densify lines or polygons\n### Aliases: smooth_densify\n\n### ** Examples\n\n# smooth_densify works on matrices of coordinates\n# use the matrix of coordinates defining a line as an example\nm <- jagged_lines$geometry[[2]][]\nm_dense <- smooth_densify(m, n = 5)\nclass(m)\nclass(m_dense)\nplot(m, type = \"b\", pch = 19, cex = 1.5, axes = FALSE, xlab = NA, ylab = NA)\npoints(m_dense, col = \"red\", pch = 19, cex = 0.5)\n\n# max_distance can be used to ensure vertices are at most a given dist apart\nm_md <- smooth_densify(m, max_distance = 0.05)\nplot(m, type = \"b\", pch = 19, cex = 1.5, axes = FALSE, xlab = NA, ylab = NA)\npoints(m_md, col = \"red\", pch = 19, cex = 0.5)\n\n# smooth is a wrapper for smooth_densify that works on spatial features\nlibrary(sf)\nl <- jagged_lines$geometry[[2]]\nl_dense <- smooth(l, method = \"densify\", n = 2)\nclass(l)\nclass(l_dense)\nplot(l, lwd = 5)\nplot(l_dense, col = \"red\", lwd = 2, lty = 2, add = TRUE)\nplot(l_dense %>% st_cast(\"MULTIPOINT\"), col = \"red\", pch = 19,\n add = TRUE)\n\n\n"} {"package":"smoothr","topic":"smooth_ksmooth","snippet":"### Name: smooth_ksmooth\n### Title: Kernel smooth\n### Aliases: smooth_ksmooth\n\n### ** Examples\n\n# smooth_ksmooth works on matrices of coordinates\n# use the matrix of coordinates defining a polygon as an example\nm <- jagged_polygons$geometry[[2]][[1]]\nm_smooth <- smooth_ksmooth(m, wrap = TRUE)\nclass(m)\nclass(m_smooth)\nplot(m, type = \"l\", col = \"black\", lwd = 3, axes = FALSE, xlab = NA,\n ylab = NA)\nlines(m_smooth, lwd = 3, col = \"red\")\n\n# lines can also be smoothed\nl <- jagged_lines$geometry[[2]][]\nl_smooth <- smooth_ksmooth(l, wrap = FALSE, max_distance = 0.05)\nplot(l, type = \"l\", col = \"black\", lwd = 3, axes = FALSE, xlab = NA,\n ylab = NA)\nlines(l_smooth, lwd = 3, col = \"red\")\n\n# explore different levels of smoothness\np <- jagged_polygons$geometry[[2]][[1]]\nps1 <- smooth_ksmooth(p, wrap = TRUE, max_distance = 0.01, smoothness = 0.5)\nps2 <- smooth_ksmooth(p, wrap = TRUE, max_distance = 0.01, smoothness = 1)\nps3 <- smooth_ksmooth(p, wrap = TRUE, max_distance = 0.01, smoothness = 2)\n# plot\npar(mar = c(0, 0, 0, 0), oma = c(10, 0, 0, 0))\nplot(p, type = \"l\", col = \"black\", lwd = 3, axes = FALSE, xlab = NA,\n ylab = NA)\nlines(ps1, lwd = 3, col = \"#E41A1C\")\nlines(ps2, lwd = 3, col = \"#4DAF4A\")\nlines(ps3, lwd = 3, col = \"#377EB8\")\npar(fig = c(0, 1, 0, 1), oma = c(0, 0, 0, 0), new = TRUE)\nplot(0, 0, type = \"n\", bty = \"n\", xaxt = \"n\", yaxt = \"n\", axes = FALSE)\nlegend(\"bottom\", legend = c(\"0.5\", \"1\", \"2\"),\n col = c(\"#E41A1C\", \"#4DAF4A\", \"#377EB8\"),\n lwd = 3, cex = 2, box.lwd = 0, inset = 0, horiz = TRUE)\n\nlibrary(sf)\np <- jagged_polygons$geometry[[2]]\np_smooth <- smooth(p, method = \"ksmooth\")\nclass(p)\nclass(p_smooth)\nplot(p_smooth, border = \"red\")\nplot(p, add = TRUE)\n\n\n"} {"package":"smoothr","topic":"smooth_spline","snippet":"### Name: smooth_spline\n### Title: Spline interpolation\n### Aliases: smooth_spline\n\n### ** Examples\n\n# smooth_spline works on matrices of coordinates\n# use the matrix of coordinates defining a polygon as an example\nm <- jagged_polygons$geometry[[2]][[1]]\nm_smooth <- smooth_spline(m, wrap = TRUE)\nclass(m)\nclass(m_smooth)\nplot(m_smooth, type = \"l\", col = \"red\", axes = FALSE, xlab = NA, ylab = NA)\nlines(m, col = \"black\")\n\n# smooth is a wrapper for smooth_spline that works on spatial features\nlibrary(sf)\np <- jagged_polygons$geometry[[2]]\np_smooth <- smooth(p, method = \"spline\")\nclass(p)\nclass(p_smooth)\nplot(p_smooth, border = \"red\")\nplot(p, add = TRUE)\n\n\n"} {"package":"agriwater","topic":"albedo_modis","snippet":"### Name: albedo_modis\n### Title: Surface Albedo using MODIS images.\n### Aliases: albedo_modis\n\n### ** Examples\n\nlibrary(agriwater)\n\n# dependencies of package 'agriwater'\nlibrary(terra)\n\n# Using a temporary folder to run example\nwd <- tempdir()\ninitial = getwd()\nsetwd(wd)\n\n# creating raster which simulate Sentinel-2 reflectances - for using\n# real data, please download:\n# https://drive.google.com/open?id=14E1wHNLxG7_Dh4I-GqNYakj8YJDgKLzk\n\nxy <- matrix(rnorm(4, mean = 0.05, sd = 0.015),2, 2)\nrast <- rast(xy, crs=\"+proj=longlat +datum=WGS84\")\next(rast) <- c(-40.5,-40.45,-9.5,-9.45)\nwriteRaster(rast, file.path(wd, \"B2.tif\"),filetype = \"GTiff\", overwrite=TRUE)\nxy <- matrix(rnorm(4, mean = 0.05, sd = 0.01),2, 2)\nrast <- rast(xy, crs=\"+proj=longlat +datum=WGS84\")\next(rast) <- c(-40.5,-40.45,-9.5,-9.45)\nwriteRaster(rast, file.path(wd, \"B1.tif\"),filetype = \"GTiff\", overwrite=TRUE)\n\n# creating mask of study area\nmask <- as.polygons(rast)\nwriteVector(mask, file.path(getwd(),\"mask.shp\"), overwrite=TRUE)\n\n# using \"agriwater\"\nalbedo_modis()\n\n#Exiting temporary folder and returning to previous workspace\nsetwd(initial)\n\n\n"} {"package":"agriwater","topic":"albedo_s2","snippet":"### Name: albedo_s2\n### Title: Surface Albedo using Sentinel-2 images.\n### Aliases: albedo_s2\n\n### ** Examples\n\nlibrary(agriwater)\n\n# dependencies of package 'agriwater'\nlibrary(terra)\n\n# Using a temporary folder to run example\nwd <- tempdir()\ninitial = getwd()\nsetwd(wd)\n\n# creating raster which simulate Sentinel-2 reflectances - for using\n# real data, please download:\n# https://drive.google.com/open?id=14E1wHNLxG7_Dh4I-GqNYakj8YJDgKLzk\n\nxy <- matrix(rnorm(4, mean = 0.07, sd = 0.01), 2, 2)\nrast <- rast(xy, crs=\"+proj=longlat +datum=WGS84\")\next(rast) <- c(-40.5,-40.45,-9.5,-9.45)\nwriteRaster(rast, file.path(wd, \"B2.tif\"),filetype = \"GTiff\", overwrite=TRUE)\nxy <- matrix(rnorm(4, mean = 0.05, sd = 0.015),2, 2)\nrast <- rast(xy, crs=\"+proj=longlat +datum=WGS84\")\next(rast) <- c(-40.5,-40.45,-9.5,-9.45)\nwriteRaster(rast, file.path(wd, \"B3.tif\"),filetype = \"GTiff\", overwrite=TRUE)\nxy <- matrix(rnorm(4, mean = 0.03, sd = 0.018),2, 2)\nrast <- rast(xy, crs=\"+proj=longlat +datum=WGS84\")\next(rast) <- c(-40.5,-40.45,-9.5,-9.45)\nwriteRaster(rast, file.path(wd, \"B4.tif\"),filetype = \"GTiff\", overwrite=TRUE)\nxy <- matrix(rnorm(4, mean = 0.05, sd = 0.015),2, 2)\nrast <- rast(xy, crs=\"+proj=longlat +datum=WGS84\")\next(rast) <- c(-40.5,-40.45,-9.5,-9.45)\nwriteRaster(rast, file.path(wd, \"B8.tif\"),filetype = \"GTiff\", overwrite=TRUE)\nmask <- as.polygons(rast)\nwriteVector(mask, file.path(getwd(),\"mask.shp\"), overwrite=TRUE)\n\n# using \"agriwater\"\nalbedo_s2()\n\n#Exiting temporary folder and returning to previous workspace\nsetwd(initial)\n\n\n"} {"package":"agriwater","topic":"evapo_modis","snippet":"### Name: evapo_modis\n### Title: Actual evapotranspiration (ETa) using MODIS with single\n### agrometeorological data.\n### Aliases: evapo_modis\n\n### ** Examples\n\nlibrary(agriwater)\n\n# dependencies of package 'agriwater'\nlibrary(terra)\n\n# Using a temporary folder to run example\nwd <- tempdir()\ninitial = getwd()\nsetwd(wd)\n\n\n# creating raster which simulate Sentinel-2 reflectances - for using\n# real data, please download:\n# https://drive.google.com/open?id=14E1wHNLxG7_Dh4I-GqNYakj8YJDgKLzk\n\nxy <- matrix(rnorm(4, mean = 0.05, sd = 0.015),2, 2)\nrast <- rast(xy, crs=\"+proj=longlat +datum=WGS84\")\next(rast) <- c(-40.5,-40.45,-9.5,-9.45)\nwriteRaster(rast, file.path(wd, \"B2.tif\"),filetype = \"GTiff\", overwrite=TRUE)\nxy <- matrix(rnorm(4, mean = 0.05, sd = 0.015),2, 2)\nrast <- rast(xy, crs=\"+proj=longlat +datum=WGS84\")\next(rast) <- c(-40.5,-40.45,-9.5,-9.45)\nwriteRaster(rast, file.path(wd, \"B1.tif\"),filetype = \"GTiff\", overwrite=TRUE)\n\nmask <- as.polygons(rast)\nwriteVector(mask, file.path(getwd(),\"mask.shp\"), overwrite=TRUE)\n\n# using \"agriwater\" - it's the same procedure as the used for\n# evapo_l8(), evapo_l8t(), evapo_modis_grid(), evapo_l8_grid(),\n# evapo_l8t_grid(), evapo_s2() and evapo_s2_grid()\nevapo_modis(doy = 134, RG = 17.6, Ta = 27.9, ET0 = 3.8, a = 1.8, b = -0.008)\n\n#Exiting temporary folder and returning to previous workspace\nsetwd(initial)\n\n\n"} {"package":"agriwater","topic":"evapo_s2","snippet":"### Name: evapo_s2\n### Title: Actual evapotranspiration (ETa) using Sentinel-2 images with\n### single agrometeorological data.\n### Aliases: evapo_s2\n\n### ** Examples\n\nlibrary(agriwater)\n\n# dependencies of package 'agriwater'\nlibrary(terra)\n\n\n# Using a temporary folder to run example\nwd <- tempdir()\ninitial = getwd()\nsetwd(wd)\n\n# creating raster which simulate Sentinel-2 reflectances - for using\n# real data, please download:\n# https://drive.google.com/open?id=14E1wHNLxG7_Dh4I-GqNYakj8YJDgKLzk\n\nxy <- matrix(rnorm(4, mean = 0.07, sd = 0.01), 2, 2)\nrast <- rast(xy, crs=\"+proj=longlat +datum=WGS84\")\next(rast) <- c(-40.5,-40.45,-9.5,-9.45)\nwriteRaster(rast, file.path(wd, \"B2.tif\"),filetype = \"GTiff\", overwrite=TRUE)\nxy <- matrix(rnorm(4, mean = 0.05, sd = 0.015),2, 2)\nrast <- rast(xy, crs=\"+proj=longlat +datum=WGS84\")\next(rast) <- c(-40.5,-40.45,-9.5,-9.45)\nwriteRaster(rast, file.path(wd, \"B3.tif\"),filetype = \"GTiff\", overwrite=TRUE)\nxy <- matrix(rnorm(4, mean = 0.03, sd = 0.018),2, 2)\nrast <- rast(xy, crs=\"+proj=longlat +datum=WGS84\")\next(rast) <- c(-40.5,-40.45,-9.5,-9.45)\nwriteRaster(rast, file.path(wd, \"B4.tif\"),filetype = \"GTiff\", overwrite=TRUE)\nxy <- matrix(rnorm(4, mean = 0.05, sd = 0.015),2, 2)\nrast <- rast(xy, crs=\"+proj=longlat +datum=WGS84\")\next(rast) <- c(-40.5,-40.45,-9.5,-9.45)\nwriteRaster(rast, file.path(wd, \"B8.tif\"),filetype = \"GTiff\", overwrite=TRUE)\nmask <- as.polygons(rast)\nwriteVector(mask, file.path(getwd(),\"mask.shp\"), overwrite=TRUE)\n\n# using \"agriwater\"\nevapo_s2(doy = 134, RG = 17.6, Ta = 27.9, ET0 = 3.8, a = 1.8, b = -0.008)\n\n#Exiting temporary folder and returning to previous workspace\nsetwd(initial)\n\n\n"} {"package":"agriwater","topic":"kc_modis","snippet":"### Name: kc_modis\n### Title: Crop coefficient (ETa / ET0) using MODIS with single\n### agrometeorological data.\n### Aliases: kc_modis\n\n### ** Examples\n\nlibrary(agriwater)\n\n# dependencies of package 'agriwater'\nlibrary(terra)\n\n\n# Using a temporary folder to run example\nwd <- tempdir()\ninitial = getwd()\nsetwd(wd)\n\n# creating raster which simulate MODIS reflectances - for using\n# real data, please download:\n# https://drive.google.com/open?id=14E1wHNLxG7_Dh4I-GqNYakj8YJDgKLzk\n\nxy <- matrix(rnorm(4, mean = 0.07, sd = 0.01), 2, 2)\nrast <- rast(xy, crs=\"+proj=longlat +datum=WGS84\")\next(rast) <- c(-40.5,-40.45,-9.5,-9.45)\nwriteRaster(rast, file.path(wd, \"B2.tif\"), filetype = \"GTiff\", overwrite=TRUE)\nxy <- matrix(rnorm(4, mean = 0.05, sd = 0.015),2, 2)\nrast <- rast(xy, crs=\"+proj=longlat +datum=WGS84\")\next(rast) <- c(-40.5,-40.45,-9.5,-9.45)\nwriteRaster(rast, file.path(wd, \"B1.tif\"), filetype = \"GTiff\", overwrite=TRUE)\n\nmask <- as.polygons(rast)\nwriteVector(mask, file.path(getwd(),\"mask.shp\"), overwrite=TRUE)\n\n# using \"agriwater\"\nkc_modis(doy = 134, RG = 17.6, Ta = 27.9, a = 1.8, b = -0.008)\n\n#Exiting temporary folder and returning to previous workspace\nsetwd(initial)\n\n\n"} {"package":"agriwater","topic":"kc_s2","snippet":"### Name: kc_s2\n### Title: Crop coefficient (ETa / ET0) using Sentinel-2 images with single\n### agrometeorological data.\n### Aliases: kc_s2\n\n### ** Examples\n\nlibrary(agriwater)\n\n# dependencies of package 'agriwater'\nlibrary(terra)\n\n# Using a temporary folder to run example\nwd <- tempdir()\ninitial = getwd()\nsetwd(wd)\n\n# creating raster which simulate Sentinel-2 reflectances - for using\n# real data, please download:\n# https://drive.google.com/open?id=14E1wHNLxG7_Dh4I-GqNYakj8YJDgKLzk\n\nxy <- matrix(rnorm(4, mean = 0.07, sd = 0.01), 2, 2)\nrast <- rast(xy, crs=\"+proj=longlat +datum=WGS84\")\next(rast) <- c(-40.5,-40.45,-9.5,-9.45)\nwriteRaster(rast, file.path(wd, \"B2.tif\"),filetype = \"GTiff\", overwrite=TRUE)\nxy <- matrix(rnorm(4, mean = 0.05, sd = 0.015),2, 2)\nrast <- rast(xy, crs=\"+proj=longlat +datum=WGS84\")\next(rast) <- c(-40.5,-40.45,-9.5,-9.45)\nwriteRaster(rast, file.path(wd, \"B3.tif\"),filetype = \"GTiff\", overwrite=TRUE)\nxy <- matrix(rnorm(4, mean = 0.03, sd = 0.018),2, 2)\nrast <- rast(xy, crs=\"+proj=longlat +datum=WGS84\")\next(rast) <- c(-40.5,-40.45,-9.5,-9.45)\nwriteRaster(rast, file.path(wd, \"B4.tif\"),filetype = \"GTiff\", overwrite=TRUE)\nxy <- matrix(rnorm(4, mean = 0.05, sd = 0.015),2, 2)\nrast <- rast(xy, crs=\"+proj=longlat +datum=WGS84\")\next(rast) <- c(-40.5,-40.45,-9.5,-9.45)\nwriteRaster(rast, file.path(wd, \"B8.tif\"),filetype = \"GTiff\", overwrite=TRUE)\nmask <- as.polygons(rast)\nwriteVector(mask, file.path(getwd(),\"mask.shp\"), overwrite=TRUE)\n\n# using \"agriwater\"\nkc_s2(doy = 134, RG = 17.6, Ta = 27.9, a = 1.8, b = -0.008)\n\n#Exiting temporary folder and returning to previous workspace\nsetwd(initial)\n\n\n"} {"package":"agriwater","topic":"radiation_modis","snippet":"### Name: radiation_modis\n### Title: Energy balance using Landsat-8 images with single\n### agrometeorological data.\n### Aliases: radiation_modis\n\n### ** Examples\n\nlibrary(agriwater)\n\n# dependencies of package 'agriwater'\nlibrary(terra)\n\n# Using a temporary folder to run example\nwd <- tempdir()\ninitial = getwd()\nsetwd(wd)\n\n# creating raster which simulate Sentinel-2 reflectances - for using\n# real data, please download:\n# https://drive.google.com/open?id=14E1wHNLxG7_Dh4I-GqNYakj8YJDgKLzk\n\nxy <- matrix(rnorm(4, mean = 0.05, sd = 0.015),2, 2)\nrast <- rast(xy, crs=\"+proj=longlat +datum=WGS84\")\next(rast) <- c(-40.5,-40.45,-9.5,-9.45)\nwriteRaster(rast, file.path(wd, \"B2.tif\"),filetype = \"GTiff\", overwrite=TRUE)\nxy <- matrix(rnorm(4, mean = 0.05, sd = 0.015),2, 2)\nrast <- rast(xy, crs=\"+proj=longlat +datum=WGS84\")\next(rast) <- c(-40.5,-40.45,-9.5,-9.45)\nwriteRaster(rast, file.path(wd, \"B1.tif\"),filetype = \"GTiff\", overwrite=TRUE)\n\n# creating mask of study area\nmask <- as.polygons(rast)\nwriteVector(mask, file.path(getwd(),\"mask.shp\"), overwrite=TRUE)\n\n# using \"agriwater\" - it's the same procedure as the used for\n# radiation_l8(), radiation_l8t(), radiation_s2(),\n# radiation_l8_grid(), radiation_l8t_grid(),\n# radiation_s2_grid(), radiation_s2() and radiation_modis_grid()\nradiation_modis(doy = 134, RG = 17.6, Ta = 27.9, ET0 = 3.8, a = 1.8, b = -0.008)\n\n#Exiting temporary folder and returning to previous workspace\nsetwd(initial)\n\n\n"} {"package":"agriwater","topic":"radiation_s2","snippet":"### Name: radiation_s2\n### Title: Energy balance using Sentinel-2 images with single\n### agrometeorological data.\n### Aliases: radiation_s2\n\n### ** Examples\n\nlibrary(agriwater)\n\n# dependencies of package 'agriwater'\nlibrary(terra)\n\n\n# Using a temporary folder to run example\nwd <- tempdir()\ninitial = getwd()\nsetwd(wd)\n\n# creating raster which simulate Sentinel-2 reflectances - for using\n# real data, please download:\n# https://drive.google.com/open?id=14E1wHNLxG7_Dh4I-GqNYakj8YJDgKLzk\n\nxy <- matrix(rnorm(4, mean = 0.07, sd = 0.01), 2, 2)\nrast <- rast(xy, crs=\"+proj=longlat +datum=WGS84\")\next(rast) <- c(-40.5,-40.45,-9.5,-9.45)\nwriteRaster(rast, file.path(wd, \"B2.tif\"),filetype = \"GTiff\", overwrite=TRUE)\nxy <- matrix(rnorm(4, mean = 0.05, sd = 0.015),2, 2)\nrast <- rast(xy, crs=\"+proj=longlat +datum=WGS84\")\next(rast) <- c(-40.5,-40.45,-9.5,-9.45)\nwriteRaster(rast, file.path(wd, \"B3.tif\"),filetype = \"GTiff\", overwrite=TRUE)\nxy <- matrix(rnorm(4, mean = 0.03, sd = 0.018),2, 2)\nrast <- rast(xy, crs=\"+proj=longlat +datum=WGS84\")\next(rast) <- c(-40.5,-40.45,-9.5,-9.45)\nwriteRaster(rast, file.path(wd, \"B4.tif\"),filetype = \"GTiff\", overwrite=TRUE)\nxy <- matrix(rnorm(4, mean = 0.05, sd = 0.015),2, 2)\nrast <- rast(xy, crs=\"+proj=longlat +datum=WGS84\")\next(rast) <- c(-40.5,-40.45,-9.5,-9.45)\nwriteRaster(rast, file.path(wd, \"B8.tif\"),filetype = \"GTiff\", overwrite=TRUE)\nmask <- as.polygons(rast)\nwriteVector(mask, file.path(getwd(),\"mask.shp\"), overwrite=TRUE)\n\n# using \"agriwater\"\nradiation_s2(doy = 134, RG = 17.6, Ta = 27.9, ET0 = 3.8, a = 1.8, b = -0.008)\n\n#Exiting temporary folder and returning to previous workspace\nsetwd(initial)\n\n\n"} {"package":"ICAOD","topic":"FIM_exp_2par","snippet":"### Name: FIM_exp_2par\n### Title: Fisher Information Matrix for the 2-Parameter Exponential Model\n### Aliases: FIM_exp_2par\n\n### ** Examples\n\nFIM_exp_2par(x = c(1, 2), w = c(.5, .5), param = c(3, 4))\n\n\n"} {"package":"ICAOD","topic":"FIM_logistic","snippet":"### Name: FIM_logistic\n### Title: Fisher Information Matrix for the 2-Parameter Logistic (2PL)\n### Model\n### Aliases: FIM_logistic\n\n### ** Examples\n\nFIM_logistic(x = c(1, 2), w = c(.5, .5), param = c(2, 1))\n\n\n"} {"package":"ICAOD","topic":"FIM_logistic_4par","snippet":"### Name: FIM_logistic_4par\n### Title: Fisher Information Matrix for the 4-Parameter Logistic Model\n### Aliases: FIM_logistic_4par\n\n### ** Examples\n\nFIM_logistic_4par(x = c(-6.9, -4.6, -3.9, 6.7 ),\n w = c(0.489, 0.40, 0.061, 0.050),\n param = c(1.563, 1.790, 8.442, 0.137))\n\n\n"} {"package":"ICAOD","topic":"FIM_mixed_inhibition","snippet":"### Name: FIM_mixed_inhibition\n### Title: Fisher Information Matrix for the Mixed Inhibition Model.\n### Aliases: FIM_mixed_inhibition\n\n### ** Examples\n\nFIM_mixed_inhibition(S = c(30, 3.86, 30, 4.60),\n I = c(0, 0, 5.11, 4.16), w = rep(.25, 4),\n param = c(1.5, 5.2, 3.4, 5.6))\n\n\n"} {"package":"ICAOD","topic":"ICA.control","snippet":"### Name: ICA.control\n### Title: Returns ICA Control Optimization Parameters\n### Aliases: ICA.control\n\n### ** Examples\n\nICA.control(ncount = 100)\n\n\n"} {"package":"ICAOD","topic":"bayes","snippet":"### Name: bayes\n### Title: Bayesian D-Optimal Designs\n### Aliases: bayes\n\n### ** Examples\n\n#############################################\n# Two parameter logistic model: uniform prior\n#############################################\n# set the unfirom prior\nuni <- uniform(lower = c(-3, .1), upper = c(3, 2))\n# set the logistic model with formula\nres1 <- bayes(formula = ~1/(1 + exp(-b *(x - a))),\n predvars = \"x\", parvars = c(\"a\", \"b\"),\n family = binomial(), lx = -3, ux = 3,\n k = 5, iter = 1, prior = uni,\n ICA.control = list(rseed = 1366))\n\n## Not run: \n##D res1 <- update(res1, 500)\n##D plot(res1)\n## End(Not run)\n# You can also use your Fisher information matrix (FIM) if you think it is faster!\n## Not run: \n##D bayes(fimfunc = FIM_logistic, lx = -3, ux = 3, k = 5, iter = 500,\n##D prior = uni, ICA.control = list(rseed = 1366))\n## End(Not run)\n\n# with fixed x\n## Not run: \n##D res1.1 <- bayes(formula = ~1/(1 + exp(-b *(x - a))),\n##D predvars = \"x\", parvars = c(\"a\", \"b\"),\n##D family = binomial(), lx = -3, ux = 3,\n##D k = 5, iter = 100, prior = uni,\n##D x = c( -3, -1.5, 0, 1.5, 3),\n##D ICA.control = list(rseed = 1366))\n##D plot(res1.1)\n##D # not optimal\n## End(Not run)\n\n# with quadrature formula\n## Not run: \n##D res1.2 <- bayes(formula = ~1/(1 + exp(-b *(x - a))),\n##D predvars = \"x\", parvars = c(\"a\", \"b\"),\n##D family = binomial(), lx = -3, ux = 3,\n##D k = 5, iter = 1, prior = uni,\n##D crt.bayes.control = list(method = \"quadrature\"),\n##D ICA.control = list(rseed = 1366))\n##D res1.2 <- update(res1.2, 500)\n##D plot(res1.2) # not optimal\n##D # compare it with res1 that was found by automatic integration\n##D plot(res1)\n##D \n##D # we increase the number of quadrature nodes\n##D res1.3 <- bayes(formula = ~1/(1 + exp(-b *(x - a))),\n##D predvars = \"x\", parvars = c(\"a\", \"b\"),\n##D family = binomial(), lx = -3, ux = 3,\n##D k = 5, iter = 1, prior = uni,\n##D crt.bayes.control = list(method = \"quadrature\",\n##D quadrature = list(level = 9)),\n##D ICA.control = list(rseed = 1366))\n##D res1.3 <- update(res1.3, 500)\n##D plot(res1.3)\n##D # by automatic integration (method = \"cubature\"),\n##D # we did not need to worry about the number of nodes.\n## End(Not run)\n###############################################\n# Two parameter logistic model: normal prior #1\n###############################################\n# defining the normal prior #1\nnorm1 <- normal(mu = c(0, 1),\n sigma = matrix(c(1, -0.17, -0.17, .5), nrow = 2),\n lower = c(-3, .1), upper = c(3, 2))\n## Not run: \n##D # initializing\n##D res2 <- bayes(formula = ~1/(1 + exp(-b *(x - a))), predvars = \"x\", parvars = c(\"a\", \"b\"),\n##D family = binomial(), lx = -3, ux = 3, k = 4, iter = 1, prior = norm1,\n##D ICA.control = list(rseed = 1366))\n##D res2 <- update(res2, 500)\n##D plot(res2)\n## End(Not run)\n\n###############################################\n# Two parameter logistic model: normal prior #2\n###############################################\n# defining the normal prior #1\nnorm2 <- normal(mu = c(0, 1),\n sigma = matrix(c(1, 0, 0, .5), nrow = 2),\n lower = c(-3, .1), upper = c(3, 2))\n## Not run: \n##D # initializing\n##D res3 <- bayes(formula = ~1/(1 + exp(-b *(x - a))), predvars = \"x\", parvars = c(\"a\", \"b\"),\n##D family = binomial(), lx = -3, ux = 3, k = 4, iter = 1, prior = norm2,\n##D ICA.control = list(rseed = 1366))\n##D \n##D res3 <- update(res3, 700)\n##D plot(res3,\n##D sens.bayes.control = list(cubature = list(maxEval = 3000, tol = 1e-4)),\n##D sens.control = list(optslist = list(maxeval = 3000)))\n## End(Not run)\n\n\n######################################################\n# Two parameter logistic model: skewed normal prior #1\n######################################################\nskew1 <- skewnormal(xi = c(0, 1),\n Omega = matrix(c(1, -0.17, -0.17, .5), nrow = 2),\n alpha = c(1, 0), lower = c(-3, .1), upper = c(3, 2))\n## Not run: \n##D res4 <- bayes(formula = ~1/(1 + exp(-b *(x - a))), predvars = \"x\", parvars = c(\"a\", \"b\"),\n##D family = binomial(), lx = -3, ux = 3, k = 4, iter = 700, prior = skew1,\n##D ICA.control = list(rseed = 1366, ncount = 60))\n##D plot(res4,\n##D sens.bayes.control = list(cubature = list(maxEval = 3000, tol = 1e-4)),\n##D sens.control = list(optslist = list(maxeval = 3000)))\n## End(Not run)\n\n\n######################################################\n# Two parameter logistic model: skewed normal prior #2\n######################################################\nskew2 <- skewnormal(xi = c(0, 1),\n Omega = matrix(c(1, -0.17, -0.17, .5), nrow = 2),\n alpha = c(-1, 0), lower = c(-3, .1), upper = c(3, 2))\n## Not run: \n##D res5 <- bayes(formula = ~1/(1 + exp(-b *(x - a))), predvars = \"x\", parvars = c(\"a\", \"b\"),\n##D family = binomial(), lx = -3, ux = 3, k = 4, iter = 700, prior = skew2,\n##D ICA.control = list(rseed = 1366, ncount = 60))\n##D plot(res5,\n##D sens.bayes.control = list(cubature = list(maxEval = 3000, tol = 1e-4)),\n##D sens.control = list(optslist = list(maxeval = 3000)))\n## End(Not run)\n\n###############################################\n# Two parameter logistic model: t student prior\n###############################################\n# set the prior\nstud <- student(mean = c(0, 1), S = matrix(c(1, -0.17, -0.17, .5), nrow = 2),\n df = 3, lower = c(-3, .1), upper = c(3, 2))\n## Not run: \n##D res6 <- bayes(formula = ~1/(1 + exp(-b *(x - a))), predvars = \"x\", parvars = c(\"a\", \"b\"),\n##D family = binomial(), lx = -3, ux = 3, k = 5, iter = 500, prior = stud,\n##D ICA.control = list(ncount = 50, rseed = 1366))\n##D plot(res6)\n## End(Not run)\n# not bad, but to find a very accurate designs we increase\n# the ncount to 200 and repeat the optimization\n## Not run: \n##D res6 <- bayes(formula = ~1/(1 + exp(-b *(x - a))),\n##D predvars = \"x\", parvars = c(\"a\", \"b\"),\n##D family = binomial(), lx = -3, ux = 3, k = 5, iter = 1000, prior = stud,\n##D ICA.control = list(ncount = 200, rseed = 1366))\n##D plot(res6)\n## End(Not run)\n\n\n##############################################\n# 4-parameter sigmoid Emax model: unform prior\n##############################################\nlb <- c(4, 11, 100, 5)\nub <- c(8, 15, 130, 9)\n## Not run: \n##D res7 <- bayes(formula = ~ theta1 + (theta2 - theta1)*(x^theta4)/(x^theta4 + theta3^theta4),\n##D predvars = c(\"x\"), parvars = c(\"theta1\", \"theta2\", \"theta3\", \"theta4\"),\n##D lx = .001, ux = 500, k = 5, iter = 200, prior = uniform(lb, ub),\n##D ICA.control = list(rseed = 1366, ncount = 60))\n##D plot(res7,\n##D sens.bayes.control = list(cubature = list(maxEval = 500, tol = 1e-3)),\n##D sens.control = list(optslist = list(maxeval = 500)))\n## End(Not run)\n\n#######################################################################\n# 2-parameter Cox Proportional-Hazards Model for type one cenosred data\n#######################################################################\n# The Fisher information matrix is available here with name FIM_2par_exp_censor1\n# However, we should reparameterize the function to match the standard of the argument 'fimfunc'\nmyfim <- function(x, w, param)\n FIM_2par_exp_censor1(x = x, w = w, param = param, tcensor = 30)\n## Not run: \n##D res8 <- bayes(fimfunc = myfim, lx = 0, ux = 1, k = 4,\n##D iter = 1, prior = uniform(c(-11, -11), c(11, 11)),\n##D ICA.control = list(rseed = 1366))\n##D \n##D res8 <- update(res8, 200)\n##D plot(res8,\n##D sens.bayes.control = list(cubature = list(maxEval = 500, tol = 1e-3)),\n##D sens.control = list(optslist = list(maxeval = 500)))\n## End(Not run)\n\n\n#######################################################################\n# 2-parameter Cox Proportional-Hazards Model for random cenosred data\n#######################################################################\n# The Fisher information matrix is available here with name FIM_2par_exp_censor2\n# However, we should reparameterize the function to match the standard of the argument 'fimfunc'\nmyfim <- function(x, w, param)\n FIM_2par_exp_censor2(x = x, w = w, param = param, tcensor = 30)\n## Not run: \n##D res9 <- bayes(fimfunc = myfim, lx = 0, ux = 1, k = 2,\n##D iter = 200, prior = uniform(c(-11, -11), c(11, 11)),\n##D ICA.control = list(rseed = 1366))\n##D plot(res9,\n##D sens.bayes.control = list(cubature = list(maxEval = 100, tol = 1e-3)),\n##D sens.control = list(optslist = list(maxeval = 100)))\n## End(Not run)\n\n#################################\n# Weibull model: Uniform prior\n################################\n# see Dette, H., & Pepelyshev, A. (2008).\n# Efficient experimental designs for sigmoidal growth models.\n# Journal of statistical planning and inference, 138(1), 2-17.\n\n## See how we fixed a some parameters in Bayesian designs\n## Not run: \n##D res10 <- bayes(formula = ~a - b * exp(-lambda * t ^h),\n##D predvars = c(\"t\"),\n##D parvars = c(\"a=1\", \"b=1\", \"lambda\", \"h=1\"),\n##D lx = .00001, ux = 20,\n##D prior = uniform(.5, 2.5), k = 5, iter = 400,\n##D ICA.control = list(rseed = 1366))\n##D plot(res10)\n## End(Not run)\n\n#################################\n# Weibull model: Normal prior\n################################\nnorm3 <- normal(mu = 1, sigma = .1, lower = .5, upper = 2.5)\nres11 <- bayes(formula = ~a - b * exp(-lambda * t ^h),\n predvars = c(\"t\"),\n parvars = c(\"a=1\", \"b=1\", \"lambda\", \"h=1\"),\n lx = .00001, ux = 20, prior = norm3, k = 4, iter = 1,\n ICA.control = list(rseed = 1366))\n\n## Not run: \n##D res11 <- update(res11, 400)\n##D plot(res11)\n## End(Not run)\n\n#################################\n# Richards model: Normal prior\n#################################\nnorm4 <- normal(mu = c(1, 1), sigma = matrix(c(.2, 0.1, 0.1, .4), 2, 2),\n lower = c(.4, .4), upper = c(1.6, 1.6))\n## Not run: \n##D res12 <- bayes(formula = ~a/(1 + b * exp(-lambda*t))^h,\n##D predvars = c(\"t\"),\n##D parvars = c(\"a=1\", \"b\", \"lambda\", \"h=1\"),\n##D lx = .00001, ux = 10,\n##D prior = norm4,\n##D k = 5, iter = 400,\n##D ICA.control = list(rseed = 1366))\n##D plot(res12,\n##D sens.bayes.control = list(cubature = list(maxEval = 1000, tol = 1e-3)),\n##D sens.control = list(optslist = list(maxeval = 1000)))\n##D ## or we can use the quadrature formula to plot the derivative function\n##D plot(res12,\n##D sens.bayes.control = list(method = \"quadrature\"),\n##D sens.control = list(optslist = list(maxeval = 1000)))\n##D \n## End(Not run)\n\n#################################\n# Exponential model: Uniform prior\n#################################\n## Not run: \n##D res13 <- bayes(formula = ~a + exp(-b*x), predvars = \"x\",\n##D parvars = c(\"a = 1\", \"b\"),\n##D lx = 0.0001, ux = 1,\n##D prior = uniform(lower = 1, upper = 20),\n##D iter = 300, k = 3,\n##D ICA.control= list(rseed = 100))\n##D plot(res13)\n## End(Not run)\n\n#################################\n# Power logistic model\n#################################\n# See, Duarte, B. P., & Wong, W. K. (2014).\n# A Semidefinite Programming based approach for finding\n# Bayesian optimal designs for nonlinear models\nuni1 <- uniform(lower = c(-.3, 6, .5), upper = c(.3, 8, 1))\n## Not run: \n##D res14 <- bayes(formula = ~1/(1 + exp(-b *(x - a)))^s, predvars = \"x\",\n##D parvars = c(\"a\", \"b\", \"s\"),\n##D lx = -1, ux = 1, prior = uni1, k = 5, iter = 1)\n##D res14 <- update(res14, 300)\n##D plot(res14)\n## End(Not run)\n\n############################################################################\n# A two-variable generalized linear model with a gamma distributed response\n############################################################################\nlb <- c(.5, 0, 0, 0, 0, 0)\nub <- c(2, 1, 1, 1, 1, 1)\nmyformula1 <- ~beta0+beta1*x1+beta2*x2+beta3*x1^2+beta4*x2^2+beta5*x1*x2\n## Not run: \n##D res15 <- bayes(formula = myformula1,\n##D predvars = c(\"x1\", \"x2\"), parvars = paste(\"beta\", 0:5, sep = \"\"),\n##D family = Gamma(),\n##D lx = rep(0, 2), ux = rep(1, 2),\n##D prior = uniform(lower = lb, upper = ub),\n##D k = 7,iter = 1, ICA.control = list(rseed = 1366))\n##D res14 <- update(res14, 500)\n##D plot(res14,\n##D sens.bayes.control = list(cubature = list(maxEval = 5000, tol = 1e-4)),\n##D sens.control = list(optslist = list(maxeval = 3000)))\n## End(Not run)\n\n#################################\n# Three parameter logistic model\n#################################\n## Not run: \n##D sigma1 <- matrix(-0.1, nrow = 3, ncol = 3)\n##D diag(sigma1) <- c(.5, .4, .1)\n##D norm5 <- normal(mu = c(0, 1, .2), sigma = sigma1,\n##D lower = c(-3, .1, 0), upper = c(3, 2, .7))\n##D res16 <- bayes(formula = ~ c + (1-c)/(1 + exp(-b *(x - a))), predvars = \"x\",\n##D parvars = c(\"a\", \"b\", \"c\"),\n##D family = binomial(), lx = -3, ux = 3,\n##D k = 4, iter = 500, prior = norm5,\n##D ICA.control = list(rseed = 1366, ncount = 50),\n##D crt.bayes.control = list(cubature = list(maxEval = 2500, tol = 1e-4)))\n##D plot(res16,\n##D sens.bayes.control = list(cubature = list(maxEval = 3000, tol = 1e-4)),\n##D sens.control = list(optslist = list(maxeval = 3000)))\n##D # took 925 second on my system\n## End(Not run)\n\n\n\n\n"} {"package":"ICAOD","topic":"bayescomp","snippet":"### Name: bayescomp\n### Title: Bayesian Compound DP-Optimal Designs\n### Aliases: bayescomp\n\n### ** Examples\n\n##########################################################################\n# DP-optimal design for a logitic model with two predictors: with formula\n##########################################################################\np <- c(1, -2, 1, -1)\nmyprior <- uniform(p -1.5, p + 1.5)\nmyformula1 <- ~exp(b0+b1*x1+b2*x2+b3*x1*x2)/(1+exp(b0+b1*x1+b2*x2+b3*x1*x2))\nres1 <- bayescomp(formula = myformula1,\n predvars = c(\"x1\", \"x2\"),\n parvars = c(\"b0\", \"b1\", \"b2\", \"b3\"),\n family = binomial(),\n lx = c(-1, -1), ux = c(1, 1),\n prior = myprior, iter = 1, k = 7,\n prob = ~1-1/(1+exp(b0 + b1 * x1 + b2 * x2 + b3 * x1 * x2)),\n alpha = .5, ICA.control = list(rseed = 1366),\n crt.bayes.control = list(cubature = list(tol = 1e-4, maxEval = 1000)))\n\n\n## Not run: \n##D res1 <- update(res1, 1000)\n##D plot(res1, sens.bayes.control = list(cubature = list(tol = 1e-3, maxEval = 1000)))\n##D # or use quadrature method\n##D plot(res1, sens.bayes.control= list(method = \"quadrature\"))\n## End(Not run)\n\n##########################################################################\n# DP-optimal design for a logitic model with two predictors: with fimfunc\n##########################################################################\n# The function of the Fisher information matrix for this model is 'FIM_logistic_2pred'\n# We should reparameterize it to match the standard of the argument 'fimfunc'\n## Not run: \n##D myfim <- function(x, w, param){\n##D npoint <- length(x)/2\n##D x1 <- x[1:npoint]\n##D x2 <- x[(npoint+1):(npoint*2)]\n##D FIM_logistic_2pred(x1 = x1,x2 = x2, w = w, param = param)\n##D }\n##D \n##D ## The following function is equivalent to the function created\n##D # by the formula: ~1-1/(1+exp(b0 + b1 * x1 + b2 * x2 + b3 * x1 * x2))\n##D # It returns probability of success given x and param\n##D # x = c(x1, x2) and param = c()\n##D \n##D myprob <- function(x, param){\n##D npoint <- length(x)/2\n##D x1 <- x[1:npoint]\n##D x2 <- x[(npoint+1):(npoint*2)]\n##D b0 <- param[1]\n##D b1 <- param[2]\n##D b2 <- param[3]\n##D b3 <- param[4]\n##D out <- 1-1/(1+exp(b0 + b1 * x1 + b2 * x2 + b3 * x1 * x2))\n##D return(out)\n##D }\n##D \n##D res2 <- bayescomp(fimfunc = myfim,\n##D lx = c(-1, -1), ux = c(1, 1),\n##D prior = myprior, iter = 1000, k = 7,\n##D prob = myprob, alpha = .5,\n##D ICA.control = list(rseed = 1366))\n##D plot(res2, sens.bayes.control = list(cubature = list(maxEval = 1000, tol = 1e-4)))\n##D # quadrature with 6 nodes (default)\n##D plot(res2, sens.bayes.control= list(method = \"quadrature\"))\n## End(Not run)\n\n\n\n\n"} {"package":"ICAOD","topic":"beff","snippet":"### Name: beff\n### Title: Calculates Relative Efficiency for Bayesian Optimal Designs\n### Aliases: beff\n\n### ** Examples\n\n#############################\n# 2PL model\n############################\nformula4.1 <- ~ 1/(1 + exp(-b *(x - a)))\npredvars4.1 <- \"x\"\nparvars4.1 <- c(\"a\", \"b\")\n\n# des4.1 is a list of Bayesian optimal designs with corresponding priors.\n\n\ndes4.1 <- vector(\"list\", 6)\ndes4.1[[1]]$x <- c(-3, -1.20829, 0, 1.20814, 3)\ndes4.1[[1]]$w <- c(.24701, .18305, .13988, .18309, .24702)\ndes4.1[[1]]$prior <- uniform(lower = c(-3, .1), upper = c(3, 2))\n\ndes4.1[[2]]$x <- c(-2.41692, -1.16676, .04386, 1.18506, 2.40631)\ndes4.1[[2]]$w <- c(.26304, .18231, .14205, .16846, .24414)\ndes4.1[[2]]$prior <- student(mean = c(0, 1), S = matrix(c(1, -0.17, -0.17, .5), nrow = 2),\n df = 3, lower = c(-3, .1), upper = c(3, 2))\n\ndes4.1[[3]]$x <- c(-2.25540, -.76318, .54628, 2.16045)\ndes4.1[[3]]$w <- c(.31762, .18225, .18159, .31853)\ndes4.1[[3]]$prior <- normal(mu = c(0, 1),\n sigma = matrix(c(1, -0.17, -0.17, .5), nrow = 2),\n lower = c(-3, .1), upper = c(3, 2))\n\ndes4.1[[4]]$x <- c(-2.23013, -.66995, .67182, 2.23055)\ndes4.1[[4]]$w <- c(.31420, .18595, .18581, .31404)\ndes4.1[[4]]$prior <- normal(mu = c(0, 1),\n sigma = matrix(c(1, 0, 0, .5), nrow = 2),\n lower = c(-3, .1), upper = c(3, 2))\n\ndes4.1[[5]]$x <- c(-1.51175, .12043, 1.05272, 2.59691)\ndes4.1[[5]]$w <- c(.37679, .14078, .12676, .35567)\ndes4.1[[5]]$prior <- skewnormal(xi = c(0, 1),\n Omega = matrix(c(1, -0.17, -0.17, .5), nrow = 2),\n alpha = c(1, 0), lower = c(-3, .1), upper = c(3, 2))\n\n\ndes4.1[[6]]$x <- c(-2.50914, -1.16780, -.36904, 1.29227)\ndes4.1[[6]]$w <- c(.35767, .11032, .15621, .37580)\ndes4.1[[6]]$prior <- skewnormal(xi = c(0, 1),\n Omega = matrix(c(1, -0.17, -0.17, .5), nrow = 2),\n alpha = c(-1, 0), lower = c(-3, .1), upper = c(3, 2))\n\n## now we want to find the relative efficiency of\n## all Bayesian optimal designs assuming different priors (6 * 6)\neff4.1 <- matrix(NA, 6, 6)\ncolnames(eff4.1) <- c(\"uni\", \"t\", \"norm1\", \"norm2\", \"skew1\", \"skew2\")\nrownames(eff4.1) <- colnames(eff4.1)\n## Not run: \n##D for (i in 1:6)\n##D for(j in 1:6)\n##D eff4.1[i, j] <- beff(formula = formula4.1,\n##D predvars = predvars4.1,\n##D parvars = parvars4.1,\n##D family = binomial(),\n##D prior = des4.1[[i]]$prior,\n##D x2 = des4.1[[i]]$x,\n##D w2 = des4.1[[i]]$w,\n##D x1 = des4.1[[j]]$x,\n##D w1 = des4.1[[j]]$w)\n##D # For example the first row represents Bayesian D-efficiencies of different\n##D # Bayesian optimal design found assuming different priors with respect to\n##D # the Bayesian D-optimal design found under uniform prior distribution.\n##D eff4.1\n## End(Not run)\n\n#############################\n# Relative efficiency for the DP-Compund criterion\n############################\np <- c(1, -2, 1, -1)\nprior4.4 <- uniform(p -1.5, p + 1.5)\nformula4.4 <- ~exp(b0+b1*x1+b2*x2+b3*x1*x2)/(1+exp(b0+b1*x1+b2*x2+b3*x1*x2))\nprob4.4 <- ~1-1/(1+exp(b0 + b1 * x1 + b2 * x2 + b3 * x1 * x2))\npredvars4.4 <- c(\"x1\", \"x2\")\nparvars4.4 <- c(\"b0\", \"b1\", \"b2\", \"b3\")\nlb <- c(-1, -1)\nub <- c(1, 1)\n\n\n\n## des4.4 is a list of DP-optimal designs found using different values for alpha\ndes4.4 <- vector(\"list\", 5)\ndes4.4[[1]]$x <- c(-1, 1)\ndes4.4[[1]]$w <- c(1)\ndes4.4[[1]]$alpha <- 0\n\n\ndes4.4[[2]]$x <- c(1, -.62534, .11405, -1, 1, .28175, -1, -1, 1, -1, -1, 1, 1, .09359)\ndes4.4[[2]]$w <- c(.08503, .43128, .01169, .14546, .05945, .08996, .17713)\ndes4.4[[2]]$alpha <- .25\n\n\ndes4.4[[3]]$x <- c(-1, .30193, 1, 1, .07411, -1, -.31952, -.08251, 1, -1, 1, -1, -1, 1)\ndes4.4[[3]]$w <- c(.09162, .10288, .15615, .13123, .01993, .22366, .27454)\ndes4.4[[3]]$alpha <- .5\n\ndes4.4[[4]]$x <- c(1, -1, .28274, 1, -1, -.19674, .03288, 1, -1, 1, -1, -.16751, 1, -1)\ndes4.4[[4]]$w <- c(.19040, .24015, .10011, .20527, .0388, .20075, .02452)\ndes4.4[[4]]$alpha <- .75\n\ndes4.4[[5]]$x <- c(1, -1, .26606, -.13370, 1, -.00887, -1, 1, -.2052, 1, 1, -1, -1, -1)\ndes4.4[[5]]$w <- c(.23020, .01612, .09546, .16197, .23675, .02701, .2325)\ndes4.4[[5]]$alpha <- 1\n\n# D-efficiency of the DP-optimal designs:\n# des4.4[[5]]$x and des4.4[[5]]$w is the D-optimal design\n\nbeff(formula = formula4.4,\n predvars = predvars4.4,\n parvars = parvars4.4,\n family = binomial(),\n prior = prior4.4,\n x2 = des4.4[[5]]$x,\n w2 = des4.4[[5]]$w,\n x1 = des4.4[[2]]$x,\n w1 = des4.4[[2]]$w)\n\nbeff(formula = formula4.4,\n predvars = predvars4.4,\n parvars = parvars4.4,\n family = binomial(),\n prior = prior4.4,\n x2 = des4.4[[5]]$x,\n w2 = des4.4[[5]]$w,\n x1 = des4.4[[3]]$x,\n w1 = des4.4[[3]]$w)\n\nbeff(formula = formula4.4,\n predvars = predvars4.4,\n parvars = parvars4.4,\n family = binomial(),\n prior = prior4.4,\n x2 = des4.4[[5]]$x,\n w2 = des4.4[[5]]$w,\n x1 = des4.4[[4]]$x,\n w1 = des4.4[[4]]$w)\n\n# must be one!\nbeff(formula = formula4.4,\n predvars = predvars4.4,\n parvars = parvars4.4,\n family = binomial(),\n prior = prior4.4,\n prob = prob4.4,\n type = \"PA\",\n x2 = des4.4[[5]]$x,\n w2 = des4.4[[5]]$w,\n x1 = des4.4[[5]]$x,\n w1 = des4.4[[5]]$w)\n\n## P-efficiency\n# reported in Table 4 as eff_P\n# des4.4[[1]]$x and des4.4[[1]]$w is the P-optimal design\nbeff(formula = formula4.4,\n predvars = predvars4.4,\n parvars = parvars4.4,\n family = binomial(),\n prior = prior4.4,\n prob = prob4.4,\n type = \"PA\",\n x2 = des4.4[[1]]$x,\n w2 = des4.4[[1]]$w,\n x1 = des4.4[[2]]$x,\n w1 = des4.4[[2]]$w)\n\nbeff(formula = formula4.4,\n predvars = predvars4.4,\n parvars = parvars4.4,\n family = binomial(),\n prior = prior4.4,\n prob = prob4.4,\n type = \"PA\",\n x2 = des4.4[[1]]$x,\n w2 = des4.4[[1]]$w,\n x1 = des4.4[[3]]$x,\n w1 = des4.4[[3]]$w)\n\nbeff(formula = formula4.4,\n predvars = predvars4.4,\n parvars = parvars4.4,\n family = binomial(),\n prior = prior4.4,\n prob = prob4.4,\n type = \"PA\",\n x2 = des4.4[[1]]$x,\n w2 = des4.4[[1]]$w,\n x1 = des4.4[[4]]$x,\n w1 = des4.4[[4]]$w)\n\nbeff(formula = formula4.4,\n predvars = predvars4.4,\n parvars = parvars4.4,\n family = binomial(),\n prior = prior4.4,\n prob = prob4.4,\n type = \"PA\",\n x2 = des4.4[[1]]$x,\n w2 = des4.4[[1]]$w,\n x1 = des4.4[[5]]$x,\n w1 = des4.4[[5]]$w)\n\n\n\n\n\n\n\n"} {"package":"ICAOD","topic":"crt.bayes.control","snippet":"### Name: crt.bayes.control\n### Title: Returns Control Parameters for Approximating Bayesian Criteria\n### Aliases: crt.bayes.control\n\n### ** Examples\n\ncrt.bayes.control()\ncrt.bayes.control(cubature = list(tol = 1e-4))\ncrt.bayes.control(quadrature = list(level = 4))\n\n\n"} {"package":"ICAOD","topic":"crt.minimax.control","snippet":"### Name: crt.minimax.control\n### Title: Returns Control Parameters for Optimizing Minimax Criteria Over\n### The Parameter Space\n### Aliases: crt.minimax.control\n\n### ** Examples\n\ncrt.minimax.control(optslist = list(maxeval = 2000))\n\n\n"} {"package":"ICAOD","topic":"leff","snippet":"### Name: leff\n### Title: Calculates Relative Efficiency for Locally Optimal Designs\n### Aliases: leff\n\n### ** Examples\n\np <- c(1, -2, 1, -1)\nprior4.4 <- uniform(p -1.5, p + 1.5)\nformula4.4 <- ~exp(b0+b1*x1+b2*x2+b3*x1*x2)/(1+exp(b0+b1*x1+b2*x2+b3*x1*x2))\nprob4.4 <- ~1-1/(1+exp(b0 + b1 * x1 + b2 * x2 + b3 * x1 * x2))\npredvars4.4 <- c(\"x1\", \"x2\")\nparvars4.4 <- c(\"b0\", \"b1\", \"b2\", \"b3\")\n\n\n# Locally D-optimal design is as follows:\n## weight and point of D-optimal design\n# Point1 Point2 Point3 Point4\n# /1.00000 \\ /-1.00000\\ /0.06801 \\ /1.00000 \\\n# \\-1.00000/ \\-1.00000/ \\1.00000 / \\1.00000 /\n# Weight1 Weight2 Weight3 Weight4\n# 0.250 0.250 0.250 0.250\n\nxopt_D <- c(1, -1, .0680, 1, -1, -1, 1, 1)\nwopt_D <- rep(.25, 4)\n\n# Let see if we use only three of the design points, what is the relative efficiency.\nleff(formula = formula4.4, predvars = predvars4.4, parvars = parvars4.4, family = binomial(),\n x1 = c(1, -1, .0680, -1, -1, 1), w1 = c(.33, .33, .33),\n inipars = p,\n x2 = xopt_D, w2 = wopt_D)\n# Wow, it heavily drops!\n\n\n# Locally P-optimal design has only one support point and is -1 and 1\nxopt_P <- c(-1, 1)\nwopt_P <- 1\n\n# What is the relative P-efficiency of the D-optimal design with respect to P-optimal design?\nleff(formula = formula4.4, predvars = predvars4.4, parvars = parvars4.4, family = binomial(),\n x1 = xopt_D, w1 = wopt_D,\n inipars = p,\n type = \"PA\",\n prob = prob4.4,\n x2 = xopt_P, w2 = wopt_P)\n# .535\n\n\n\n\n"} {"package":"ICAOD","topic":"locally","snippet":"### Name: locally\n### Title: Locally D-Optimal Designs\n### Aliases: locally\n\n### ** Examples\n\n#################################\n# Exponential growth model\n################################\n# See how we set stopping rule by adjusting 'stop_rule', 'checkfreq' and 'stoptol'\n# It calls the 'senslocally' function every checkfreq = 50 iterations to\n# calculate the ELB. if ELB is greater than stoptol = .95, then the algoithm stops.\n\n# initializing by one iteration\nres1 <- locally(formula = ~a + exp(-b*x), predvars = \"x\", parvars = c(\"a\", \"b\"),\n lx = 0, ux = 1, inipars = c(1, 10),\n iter = 1, k = 2,\n ICA.control= ICA.control(rseed = 100,\n stop_rule = \"equivalence\",\n checkfreq = 20, stoptol = .95))\n## Not run: \n##D # update the algorithm\n##D res1 <- update(res1, 150)\n##D #stops at iteration 21 because ELB is greater than .95\n## End(Not run)\n\n### fixed x, lx and ux are only required for equivalence theorem\n## Not run: \n##D res1.1 <- locally(formula = ~a + exp(-b*x), predvars = \"x\", parvars = c(\"a\", \"b\"),\n##D lx = 0, ux = 1, inipars = c(1, 10),\n##D iter = 100,\n##D x = c(.25, .5, .75),\n##D ICA.control= ICA.control(rseed = 100))\n##D plot(res1.1)\n##D # we can not have an optimal design using this x\n## End(Not run)\n\n################################\n## two parameter logistic model\n################################\nres2 <- locally(formula = ~1/(1 + exp(-b *(x - a))),\n predvars = \"x\", parvars = c(\"a\", \"b\"),\n family = binomial(), lx = -3, ux = 3,\n inipars = c(1, 3), iter = 1, k = 2,\n ICA.control= list(rseed = 100, stop_rule = \"equivalence\",\n checkfreq = 50, stoptol = .95))\n## Not run: \n##D res2 <- update(res2, 100)\n##D # stops at iteration 51\n## End(Not run)\n\n\n\n\n################################\n# A model with two predictors\n################################\n# mixed inhibition model\n## Not run: \n##D res3 <- locally(formula = ~ V*S/(Km * (1 + I/Kic)+ S * (1 + I/Kiu)),\n##D predvars = c(\"S\", \"I\"),\n##D parvars = c(\"V\", \"Km\", \"Kic\", \"Kiu\"),\n##D family = gaussian(),\n##D lx = c(0, 0), ux = c(30, 60),\n##D k = 4,\n##D iter = 300,\n##D inipars = c(1.5, 5.2, 3.4, 5.6),\n##D ICA.control= list(rseed = 100, stop_rule = \"equivalence\",\n##D checkfreq = 50, stoptol = .95))\n##D # stops at iteration 100\n## End(Not run)\n\n\n## Not run: \n##D # fixed x\n##D res3.1 <- locally(formula = ~ V*S/(Km * (1 + I/Kic)+ S * (1 + I/Kiu)),\n##D predvars = c(\"S\", \"I\"),\n##D parvars = c(\"V\", \"Km\", \"Kic\", \"Kiu\"),\n##D family = gaussian(),\n##D lx = c(0, 0), ux = c(30, 60),\n##D iter = 100,\n##D x = c(20, 4, 20, 4, 10, 0, 0, 30, 3, 2),\n##D inipars = c(1.5, 5.2, 3.4, 5.6),\n##D ICA.control= list(rseed = 100))\n## End(Not run)\n\n\n###################################\n# user-defined optimality criterion\n##################################\n# When the model is defined by the formula interface\n# A-optimal design for the 2PL model.\n# the criterion function must have argument x, w fimfunc and the parameters defined in 'parvars'.\n# use 'fimfunc' as a function of the design points x, design weights w and\n# the 'parvars' parameters whenever needed.\nAopt <-function(x, w, a, b, fimfunc){\n sum(diag(solve(fimfunc(x = x, w = w, a = a, b = b))))\n}\n## the sensitivtiy function\n# xi_x is a design that put all its mass on x in the definition of the sensitivity function\n# x is a vector of design points\nAopt_sens <- function(xi_x, x, w, a, b, fimfunc){\n fim <- fimfunc(x = x, w = w, a = a, b = b)\n M_inv <- solve(fim)\n M_x <- fimfunc(x = xi_x, w = 1, a = a, b = b)\n sum(diag(M_inv %*% M_x %*% M_inv)) - sum(diag(M_inv))\n}\n\nres4 <- locally(formula = ~1/(1 + exp(-b * (x-a))), predvars = \"x\",\n parvars = c(\"a\", \"b\"), family = \"binomial\",\n lx = -3, ux = 3, inipars = c(1, 1.25),\n iter = 1, k = 2,\n crtfunc = Aopt,\n sensfunc = Aopt_sens,\n ICA.control = list(checkfreq = Inf))\n## Not run: \n##D res4 <- update(res4, 50)\n## End(Not run)\n\n# When the FIM of the model is defined directly via the argument 'fimfunc'\n# the criterion function must have argument x, w fimfunc and param.\n# use 'fimfunc' as a function of the design points x, design weights w\n# and param whenever needed.\nAopt2 <-function(x, w, param, fimfunc){\n sum(diag(solve(fimfunc(x = x, w = w, param = param))))\n}\n## the sensitivtiy function\n# xi_x is a design that put all its mass on x in the definition of the sensitivity function\n# x is a vector of design points\nAopt_sens2 <- function(xi_x, x, w, param, fimfunc){\n fim <- fimfunc(x = x, w = w, param = param)\n M_inv <- solve(fim)\n M_x <- fimfunc(x = xi_x, w = 1, param = param)\n sum(diag(M_inv %*% M_x %*% M_inv)) - sum(diag(M_inv))\n}\n\nres4.1 <- locally(fimfunc = FIM_logistic,\n lx = -3, ux = 3, inipars = c(1, 1.25),\n iter = 1, k = 2,\n crtfunc = Aopt2,\n sensfunc = Aopt_sens2,\n ICA.control = list(checkfreq = Inf))\n## Not run: \n##D res4.1 <- update(res4.1, 50)\n##D plot(res4.1)\n## End(Not run)\n\n\n# locally c-optimal design\n# example from Chaloner and Larntz (1989) Figure 3\nc_opt <-function(x, w, a, b, fimfunc){\n gam <- log(.95/(1-.95))\n M <- fimfunc(x = x, w = w, a = a, b = b)\n c <- matrix(c(1, -gam * b^(-2)), nrow = 1)\n B <- t(c) %*% c\n sum(diag(B %*% solve(M)))\n}\n\nc_sens <- function(xi_x, x, w, a, b, fimfunc){\n gam <- log(.95/(1-.95))\n M <- fimfunc(x = x, w = w, a = a, b = b)\n M_inv <- solve(M)\n M_x <- fimfunc(x = xi_x, w = 1, a = a, b = b)\n c <- matrix(c(1, -gam * b^(-2)), nrow = 1)\n B <- t(c) %*% c\n sum(diag(B %*% M_inv %*% M_x %*% M_inv)) - sum(diag(B %*% M_inv))\n}\n\n\nres4.2 <- locally(formula = ~1/(1 + exp(-b * (x-a))), predvars = \"x\",\n parvars = c(\"a\", \"b\"), family = \"binomial\",\n lx = -1, ux = 1, inipars = c(0, 7),\n iter = 1, k = 2,\n crtfunc = c_opt, sensfunc = c_sens,\n ICA.control = list(rseed = 1, checkfreq = Inf))\n## Not run: \n##D res4.2 <- update(res4.2, 100)\n## End(Not run)\n\n\n\n"} {"package":"ICAOD","topic":"locallycomp","snippet":"### Name: locallycomp\n### Title: Locally DP-Optimal Designs\n### Aliases: locallycomp\n\n### ** Examples\n\n## Here we produce the results of Table 2 in in McGree and Eccleston (2008)\n# For D- and P-efficiency see, ?leff and ?peff\n\np <- c(1, -2, 1, -1)\nprior4.4 <- uniform(p -1.5, p + 1.5)\nformula4.4 <- ~exp(b0+b1*x1+b2*x2+b3*x1*x2)/(1+exp(b0+b1*x1+b2*x2+b3*x1*x2))\nprob4.4 <- ~1-1/(1+exp(b0 + b1 * x1 + b2 * x2 + b3 * x1 * x2))\npredvars4.4 <- c(\"x1\", \"x2\")\nparvars4.4 <- c(\"b0\", \"b1\", \"b2\", \"b3\")\nlb <- c(-1, -1)\nub <- c(1, 1)\n\n\n# set checkfreq = Inf to ask for equivalence theorem at final step.\nres.0 <- locallycomp(formula = formula4.4, predvars = predvars4.4, parvars = parvars4.4,\n family = binomial(), prob = prob4.4, lx = lb, ux = ub,\n alpha = 0, k = 1, inipars = p, iter = 10,\n ICA.control = ICA.control(checkfreq = Inf))\n\n## Not run: \n##D res.25 <- locallycomp(formula = formula4.4, predvars = predvars4.4, parvars = parvars4.4,\n##D family = binomial(), prob = prob4.4, lx = lb, ux = ub,\n##D alpha = .25, k = 4, inipars = p, iter = 350,\n##D ICA.control = ICA.control(checkfreq = Inf))\n##D \n##D res.5 <- locallycomp(formula = formula4.4, predvars = predvars4.4, parvars = parvars4.4,\n##D family = binomial(), prob = prob4.4, lx = lb, ux = ub,\n##D alpha = .5, k = 4, inipars = p, iter = 350,\n##D ICA.control = ICA.control(checkfreq = Inf))\n##D res.75 <- locallycomp(formula = formula4.4, predvars = predvars4.4, parvars = parvars4.4,\n##D family = binomial(), prob = prob4.4, lx = lb, ux = ub,\n##D alpha = .75, k = 4, inipars = p, iter = 350,\n##D ICA.control = ICA.control(checkfreq = Inf))\n##D \n##D res.1 <- locallycomp(formula = formula4.4, predvars = predvars4.4, parvars = parvars4.4,\n##D family = binomial(), prob = prob4.4, lx = lb, ux = ub,\n##D alpha = 1, k = 4, inipars = p, iter = 350,\n##D ICA.control = ICA.control(checkfreq = Inf))\n##D \n##D #### computing the D-efficiency\n##D # locally D-optimal design is locally DP-optimal design when alpha = 1.\n##D \n##D leff(formula = formula4.4, predvars = predvars4.4, parvars = parvars4.4, family = binomial(),\n##D x1 = res.0$evol[[10]]$x, w1 = res.0$evol[[10]]$w,\n##D inipars = p,\n##D x2 = res.1$evol[[350]]$x, w2 = res.1$evol[[350]]$w)\n##D \n##D leff(formula = formula4.4, predvars = predvars4.4, parvars = parvars4.4, family = binomial(),\n##D x1 = res.25$evol[[350]]$x, w1 = res.25$evol[[350]]$w,\n##D inipars = p,\n##D x2 = res.1$evol[[350]]$x, w2 = res.1$evol[[350]]$w)\n##D \n##D leff(formula = formula4.4, predvars = predvars4.4, parvars = parvars4.4, family = binomial(),\n##D x1 = res.5$evol[[350]]$x, w1 = res.5$evol[[350]]$w,\n##D inipars = p,\n##D x2 = res.1$evol[[350]]$x, w2 = res.1$evol[[350]]$w)\n##D \n##D \n##D leff(formula = formula4.4, predvars = predvars4.4, parvars = parvars4.4, family = binomial(),\n##D x1 = res.75$evol[[350]]$x, w1 = res.75$evol[[350]]$w,\n##D inipars = p,\n##D x2 = res.1$evol[[350]]$x, w2 = res.1$evol[[350]]$w)\n##D \n##D \n##D \n##D #### computing the P-efficiency\n##D # locally p-optimal design is locally DP-optimal design when alpha = 0.\n##D \n##D leff(formula = formula4.4, predvars = predvars4.4, parvars = parvars4.4, family = binomial(),\n##D x2 = res.0$evol[[10]]$x, w2 = res.0$evol[[10]]$w,\n##D prob = prob4.4,\n##D type = \"PA\",\n##D inipars = p,\n##D x1 = res.25$evol[[350]]$x, w1 = res.25$evol[[350]]$w)\n##D \n##D leff(formula = formula4.4, predvars = predvars4.4, parvars = parvars4.4, family = binomial(),\n##D x2 = res.0$evol[[10]]$x, w2 = res.0$evol[[10]]$w,\n##D prob = prob4.4,\n##D inipars = p,\n##D type = \"PA\",\n##D x1 = res.5$evol[[350]]$x, w1 = res.5$evol[[350]]$w)\n##D \n##D leff(formula = formula4.4, predvars = predvars4.4, parvars = parvars4.4, family = binomial(),\n##D x2 = res.0$evol[[10]]$x, w2 = res.0$evol[[10]]$w,\n##D prob = prob4.4,\n##D inipars = p,\n##D type = \"PA\",\n##D x1 = res.75$evol[[350]]$x, w1 = res.75$evol[[350]]$w)\n##D \n##D \n##D leff(formula = formula4.4, predvars = predvars4.4, parvars = parvars4.4, family = binomial(),\n##D x2 = res.0$evol[[10]]$x, w2 = res.1$evol[[10]]$w,\n##D prob = prob4.4,\n##D type = \"PA\",\n##D inipars = p,\n##D x1 = res.1$evol[[350]]$x, w1 = res.1$evol[[350]]$w)\n##D \n##D \n## End(Not run)\n\n\n"} {"package":"ICAOD","topic":"meff","snippet":"### Name: meff\n### Title: Calculates Relative Efficiency for Minimax Optimal Designs\n### Aliases: meff\n\n### ** Examples\n\n# Relative D-efficiency with respect to the minimax criterion\nmeff(formula = ~1/(1 + exp(-b * (x-a))), predvars = \"x\",\n parvars = c(\"a\", \"b\"), family = \"binomial\",\n lp = c(-3, .5), up = c(3, 2),\n x2 = c(-3, -1.608782, 0, 1.608782, 3),\n w2 = c(0.22291601, 0.26438449, 0.02539899, 0.26438449, 0.22291601),\n x1 = c(-1, 1), w1 = c(.5, .5))\n\n\n\n# A function to calculate the locally D-optimal design for the 2PL model\nDopt_2pl <- function(a, b){\n x <- c(a + (1/b) * 1.5434046, a - (1/b) * 1.5434046)\n return(list(x = x, w = c(.5, .5)))\n}\n# Relative D-efficiency with respect to the standardized maximin criterion\nmeff (formula = ~1/(1 + exp(-b * (x-a))), predvars = \"x\",\n parvars = c(\"a\", \"b\"), family = \"binomial\",\n lp = c(-3, .5), up = c(3, 2),\n x2 = c(-3, -1.611255, 0, 1.611255, 3),\n w2 = c(0.22167034, 0.26592974, 0.02479984, 0.26592974, 0.22167034),\n x1 = c(0, -1), w1 = c(.5, .5),\n standardized = TRUE,\n localdes = Dopt_2pl)\n\n\n\n\n"} {"package":"ICAOD","topic":"minimax","snippet":"### Name: minimax\n### Title: Minimax and Standardized Maximin D-Optimal Designs\n### Aliases: minimax\n\n### ** Examples\n\n########################################\n# Two-parameter exponential growth model\n########################################\nres1 <- minimax (formula = ~a + exp(-b*x), predvars = \"x\", parvars = c(\"a\", \"b\"),\n lx = 0, ux = 1, lp = c(1, 1), up = c(1, 10),\n iter = 1, k = 4,\n ICA.control= ICA.control(rseed = 100),\n crt.minimax.control = list(optslist = list(maxeval = 100)))\n# The optimal design has 3 points, but we set k = 4 for illustration purpose to\n# show how the algorithm modifies the design by adjusting the weights\n# The value of maxeval is changed to reduce the CPU time\n## Not run: \n##D res1 <- update(res1, 150)\n##D # iterating the algorithm up to 150 more iterations\n## End(Not run)\n\nres1 # print method\nplot(res1) # Veryfying the general equivalence theorem\n\n## Not run: \n##D ## fixed x\n##D res1.1 <- minimax (formula = ~a + exp(-b*x), predvars = \"x\", parvars = c(\"a\", \"b\"),\n##D lx = 0, ux = 1, lp = c(1, 1), up = c(1, 10),\n##D x = c(0, .5, 1),\n##D iter = 150, k = 3, ICA.control= ICA.control(rseed = 100))\n##D # not optimal\n## End(Not run)\n\n########################################\n# Two-parameter logistic model.\n########################################\n# A little playing with the tuning parameters\n# The value of maxeval is reduced to 200 to increase the speed\ncont1 <- crt.minimax.control(optslist = list(maxeval = 200))\ncont2 <- ICA.control(rseed = 100, checkfreq = Inf, ncount = 60)\n\n## Not run: \n##D res2 <- minimax (formula = ~1/(1 + exp(-b *(x - a))), predvars = \"x\",\n##D parvars = c(\"a\", \"b\"),\n##D family = binomial(), lx = -3, ux = 3,\n##D lp = c(0, 1), up = c(1, 2.5), iter = 200, k = 3,\n##D ICA.control= cont2, crt.minimax.control = cont1)\n##D print(res2)\n##D plot(res2)\n## End(Not run)\n\n############################################\n# An example of a model with two predictors\n############################################\n# Mixed inhibition model\nlower <- c(1, 4, 2, 4)\nupper <- c(1, 5, 3, 5)\ncont <- crt.minimax.control(optslist = list(maxeval = 100)) # to be faster\n## Not run: \n##D res3 <- minimax(formula = ~ V*S/(Km * (1 + I/Kic)+ S * (1 + I/Kiu)),\n##D predvars = c(\"S\", \"I\"),\n##D parvars = c(\"V\", \"Km\", \"Kic\", \"Kiu\"),\n##D lx = c(0, 0), ux = c(30, 60), k = 4,\n##D iter = 100, lp = lower, up = upper,\n##D ICA.control= list(rseed = 100),\n##D crt.minimax.control = cont)\n##D \n##D res3 <- update(res3, 100)\n##D print(res3)\n##D plot(res3) # sensitivity plot\n##D res3$arg$time\n## End(Not run)\n\n# Now consider grid points instead of assuming continuous parameter space\n# set n.grid to 5\n## Not run: \n##D res4 <- minimax(formula = ~ V*S/(Km * (1 + I/Kic)+ S * (1 + I/Kiu)),\n##D predvars = c(\"S\", \"I\"),\n##D parvars = c(\"V\", \"Km\", \"Kic\", \"Kiu\"),\n##D lx = c(0, 0), ux = c(30, 60),\n##D k = 4, iter = 130, n.grid = 5, lp = lower, up = upper,\n##D ICA.control= list(rseed = 100, checkfreq = Inf),\n##D crt.minimax.control = cont)\n##D print(res4)\n##D plot(res4) # sensitivity plot\n## End(Not run)\n\n############################################\n# Standardized maximin D-optimal designs\n############################################\n# Assume the purpose is finding STANDARDIZED designs\n# We know from literature that the locally D-optimal design (LDOD)\n# for this model has an analytical solution.\n# The follwoing function takes the parameter as input and returns\n# the design points and weights of LDOD.\n# x and w are exactly similar to the arguments of 'fimfunc'.\n# x is a vector and returns the design points 'dimension-wise'.\n# see explanation of the arguments of 'fimfunc' in 'Details'.\n\nLDOD <- function(V, Km, Kic, Kiu){\n #first dimention is for S and the second one is for I.\n S_min <- 0\n S_max <- 30\n I_min <- 0\n I_max <- 60\n s2 <- max(S_min, S_max*Km*Kiu*(Kic+I_min)/\n (S_max*Kic*I_min+S_max*Kic*Kiu+2*Km*Kiu*I_min+2*Km*Kiu*Kic))\n i3 <- min((2*S_max*Kic*I_min + S_max*Kic*Kiu+2*Km*Kiu*I_min+Km*Kiu*Kic)/\n (Km*Kiu+S_max*Kic), I_max)\n i4 <- min(I_min + (sqrt((Kic+I_min)*(Km*Kic*Kiu+Km*Kiu*I_min+\n S_max*Kic*Kiu+S_max*Kic*I_min)/\n (Km*Kiu+S_max*Kic))), I_max )\n s4 <- max(-Km*Kiu*(Kic+2*I_min-i4)/(Kic*(Kiu+2*I_min-i4)), S_min)\n x <- c(S_max, s2, S_max, s4, I_min, I_min, i3, i4)\n return(list(x = x, w =rep(1/4, 4)))\n\n}\nformalArgs(LDOD)\n## Not run: \n##D minimax(formula = ~ V*S/(Km * (1 + I/Kic)+ S * (1 + I/Kiu)),\n##D predvars = c(\"S\", \"I\"),\n##D parvars = c(\"V\", \"Km\", \"Kic\", \"Kiu\"),\n##D lx = c(0, 0), ux = c(30, 60),\n##D k = 4, iter = 300,\n##D lp = lower, up = upper,\n##D ICA.control= list(rseed = 100, checkfreq = Inf),\n##D crt.minimax.control = cont,\n##D standardized = TRUE,\n##D localdes = LDOD)\n## End(Not run)\n\n\n################################################################\n# Not necessary!\n# The rest of the examples here are only for professional uses.\n################################################################\n# Imagine you have written your own FIM, say in Rcpp that is faster than\n# the FIM created by the formula interface above.\n\n###########################################\n# An example of a model with two predictors\n###########################################\n# For example, th cpp FIM function for the mixed inhibition model is named:\nformalArgs(FIM_mixed_inhibition)\n\n# We should reparamterize the arguments to match the standard of the\n# argument 'fimfunc' (see 'Details').\nmyfim <- function(x, w, param){\n npoint <- length(x)/2\n S <- x[1:npoint]\n I <- x[(npoint+1):(npoint*2)]\n out <- FIM_mixed_inhibition(S = S, I = I, w = w, param = param)\n return(out)\n}\nformalArgs(myfim)\n\n# Finds minimax optimal design, exactly as before, but NOT using the\n# formula interface.\n## Not run: \n##D res5 <- minimax(fimfunc = myfim,\n##D lx = c(0, 0), ux = c(30, 60), k = 4,\n##D iter = 100, lp = lower, up = upper,\n##D ICA.control= list(rseed = 100),\n##D crt.minimax.control = cont)\n##D print(res5)\n##D plot(res5) # sensitivity plot\n## End(Not run)\n#########################################\n# Standardized maximin D-optimal designs\n#########################################\n# To match the argument 'localdes' when no formula inteface is used,\n# we should reparameterize LDOD.\n# The input must be 'param' same as the argument of 'fimfunc'\nLDOD2 <- function(param)\n LDOD(V = param[1], Km = param[2], Kic = param[3], Kiu = param[4])\n\n# compare these two:\nformalArgs(LDOD)\nformalArgs(LDOD2)\n## Not run: \n##D res6 <- minimax(fimfunc = myfim,\n##D lx = c(0, 0), ux = c(30, 60), k = 4,\n##D iter = 300, lp = lower, up = upper,\n##D ICA.control= list(rseed = 100, checkfreq = Inf),\n##D crt.minimax.control = cont,\n##D standardized = TRUE,\n##D localdes = LDOD2)\n##D res6\n##D plot(res6)\n## End(Not run)\n\n###################################\n# user-defined optimality criterion\n##################################\n# When the model is defined by the formula interface\n# A-optimal design for the 2PL model.\n# the criterion function must have argument x, w fimfunc and the parameters defined in 'parvars'.\n# use 'fimfunc' as a function of the design points x, design weights w and\n# the 'parvars' parameters whenever needed.\nAopt <-function(x, w, a, b, fimfunc){\n sum(diag(solve(fimfunc(x = x, w = w, a = a, b = b))))\n}\n## the sensitivtiy function\n# xi_x is a design that put all its mass on x in the definition of the sensitivity function\n# x is a vector of design points\nAopt_sens <- function(xi_x, x, w, a, b, fimfunc){\n fim <- fimfunc(x = x, w = w, a = a, b = b)\n M_inv <- solve(fim)\n M_x <- fimfunc(x = xi_x, w = 1, a = a, b = b)\n sum(diag(M_inv %*% M_x %*% M_inv)) - sum(diag(M_inv))\n}\n## Not run: \n##D res7 <- minimax(formula = ~1/(1 + exp(-b * (x-a))), predvars = \"x\",\n##D parvars = c(\"a\", \"b\"), family = \"binomial\",\n##D lx = -2, ux = 2,\n##D lp = c(-2, 1), up = c(2, 1.5),\n##D iter = 400, k = 3,\n##D crtfunc = Aopt,\n##D sensfunc = Aopt_sens,\n##D crt.minimax.control = list(optslist = list(maxeval = 200)),\n##D ICA.control = list(rseed = 1))\n##D plot(res7)\n## End(Not run)\n# with grid points\nres7.1 <- minimax(formula = ~1/(1 + exp(-b * (x-a))), predvars = \"x\",\n parvars = c(\"a\", \"b\"), family = \"binomial\",\n lx = -2, ux = 2,\n lp = c(-2, 1), up = c(2, 1.5),\n iter = 1, k = 3,\n crtfunc = Aopt,\n sensfunc = Aopt_sens,\n n.grid = 9,\n ICA.control = list(rseed = 1))\n## Not run: \n##D res7.1 <- update(res7.1, 400)\n##D plot(res7.1)\n## End(Not run)\n\n# When the FIM of the model is defined directly via the argument 'fimfunc'\n# the criterion function must have argument x, w fimfunc and param.\n# use 'fimfunc' as a function of the design points x, design weights w and\n# the 'parvars' parameters whenever needed.\nAopt2 <-function(x, w, param, fimfunc){\n sum(diag(solve(fimfunc(x = x, w = w, param = param))))\n}\n## the sensitivtiy function\n# xi_x is a design that put all its mass on x in the definition of the sensitivity function\n# x is a vector of design points\nAopt_sens2 <- function(xi_x, x, w, param, fimfunc){\n fim <- fimfunc(x = x, w = w, param = param)\n M_inv <- solve(fim)\n M_x <- fimfunc(x = xi_x, w = 1, param = param)\n sum(diag(M_inv %*% M_x %*% M_inv)) - sum(diag(M_inv))\n}\n## Not run: \n##D res7.2 <- minimax(fimfunc = FIM_logistic,\n##D lx = -2, ux = 2,\n##D lp = c(-2, 1), up = c(2, 1.5),\n##D iter = 1, k = 3,\n##D crtfunc = Aopt2,\n##D sensfunc = Aopt_sens2,\n##D crt.minimax.control = list(optslist = list(maxeval = 200)),\n##D ICA.control = list(rseed = 1))\n##D res7.2 <- update(res7.2, 200)\n##D plot(res7.2)\n## End(Not run)\n# with grid points\nres7.3 <- minimax(fimfunc = FIM_logistic,\n lx = -2, ux = 2,\n lp = c(-2, 1), up = c(2, 1.5),\n iter = 1, k = 3,\n crtfunc = Aopt2,\n sensfunc = Aopt_sens2,\n n.grid = 9,\n ICA.control = list(rseed = 1))\n## Not run: \n##D res7.3 <- update(res7.2, 200)\n##D plot(res7.3)\n## End(Not run)\n\n\n# robust c-optimal design\n# example from Chaloner and Larntz (1989), Figure 3, but robust design\nc_opt <-function(x, w, a, b, fimfunc){\n gam <- log(.95/(1-.95))\n M <- fimfunc(x = x, w = w, a = a, b = b)\n c <- matrix(c(1, -gam * b^(-2)), nrow = 1)\n B <- t(c) %*% c\n sum(diag(B %*% solve(M)))\n}\n\nc_sens <- function(xi_x, x, w, a, b, fimfunc){\n gam <- log(.95/(1-.95))\n M <- fimfunc(x = x, w = w, a = a, b = b)\n M_inv <- solve(M)\n M_x <- fimfunc(x = xi_x, w = 1, a = a, b = b)\n c <- matrix(c(1, -gam * b^(-2)), nrow = 1)\n B <- t(c) %*% c\n sum(diag(B %*% M_inv %*% M_x %*% M_inv)) - sum(diag(B %*% M_inv))\n}\n\n## Not run: \n##D res8 <- minimax(formula = ~1/(1 + exp(-b * (x-a))), predvars = \"x\",\n##D parvars = c(\"a\", \"b\"), family = \"binomial\",\n##D lx = -1, ux = 1,\n##D lp = c(-.3, 6), up = c(.3, 8),\n##D iter = 500, k = 3,\n##D crtfunc = c_opt, sensfunc = c_sens,\n##D ICA.control = list(rseed = 1, ncount = 100),\n##D n.grid = 12)\n##D plot(res8)\n## End(Not run)\n\n\n\n\n\n"} {"package":"ICAOD","topic":"multiple","snippet":"### Name: multiple\n### Title: Locally Multiple Objective Optimal Designs for the 4-Parameter\n### Hill Model\n### Aliases: multiple\n\n### ** Examples\n\n# All the examples are available in Hyun and Wong (2015)\n\n#################################\n# 4-parameter logistic model\n# Example 1, Table 3\n#################################\nlam <- c(0.05, 0.05, .90)\n# Initial estimates are derived from Table 1\n# See how the stopping rules are set via 'stop_rul', checkfreq' and 'stoptol'\nTheta1 <- c(1.563, 1.790, 8.442, 0.137)\nres1 <- multiple(minDose = log(.001), maxDose = log(1000),\n inipars = Theta1, k = 4, lambda = lam, delta = -1,\n Hill_par = FALSE,\n iter = 1,\n ICA.control = list(rseed = 1366, ncount = 100,\n stop_rule = \"equivalence\",\n checkfreq = 100, stoptol = .95))\n## Not run: \n##D res1 <- update(res1, 1000)\n##D # stops at iteration 101\n## End(Not run)\n\n#################################\n# 4-parameter Hill model\n#################################\n## initial estimates for the parameters of Hill model:\na <- 0.008949 # ED50\nb <- -1.79 # Hill constant\nc <- 0.137 # lower limit\nd <- 1.7 # upper limit\n# D belongs to c(.001, 1000) ## dose in mg\n## the vector of Hill parameters are now c(a, b, c, d)\n## Not run: \n##D res2 <- multiple(minDose = .001, maxDose = 1000,\n##D inipars = c(a, b, c, d),\n##D Hill_par = TRUE, k = 4, lambda = lam,\n##D delta = -1, iter = 1000,\n##D ICA.control = list(rseed = 1366, ncount = 100,\n##D stop_rule = \"equivalence\",\n##D checkfreq = 100, stoptol = .95))\n##D # stops at iteration 100\n## End(Not run)\n\n\n\n# use x argument to provide fix number of dose levels.\n# In this case, the optimization is only over weights\n## Not run: \n##D res3 <- multiple(minDose = log(.001), maxDose = log(1000),\n##D inipars = Theta1, k = 4, lambda = lam, delta = -1,\n##D iter = 300,\n##D Hill_par = FALSE,\n##D x = c(-6.90, -4.66, -3.93, 3.61),\n##D ICA.control = list(rseed = 1366))\n##D res3$evol[[300]]$w\n##D # if the user provide the desugn points via x, there is no guarantee\n##D # that the resulted design is optimal. It only provides the optimal weights given\n##D # the x points of the design.\n##D plot(res3)\n## End(Not run)\n\n\n\n"} {"package":"ICAOD","topic":"normal","snippet":"### Name: normal\n### Title: Assumes A Multivariate Normal Prior Distribution for The Model\n### Parameters\n### Aliases: normal\n\n### ** Examples\n\nnormal(mu = c(0, 1), sigma = matrix(c(1, -0.17, -0.17, .5), nrow = 2),\n lower = c(-3, .1), upper = c(3, 2))\n\n\n"} {"package":"ICAOD","topic":"robust","snippet":"### Name: robust\n### Title: Robust D-Optimal Designs\n### Aliases: robust\n\n### ** Examples\n\n# Finding a robust design for the two-parameter logistic model\n# See how we set a stopping rule.\n# The ELB is computed every checkfreq = 30 iterations\n# The optimization stops when the ELB is larger than stoptol = .95\nres1 <- robust(formula = ~1/(1 + exp(-b *(x - a))),\n predvars = c(\"x\"), parvars = c(\"a\", \"b\"),\n family = binomial(),\n lx = -5, ux = 5, prob = rep(1/4, 4),\n parset = matrix(c(0.5, 1.5, 0.5, 1.5, 4.0, 4.0, 5.0, 5.0), 4, 2),\n iter = 1, k =3,\n ICA.control = list(stop_rule = \"equivalence\",\n stoptol = .95, checkfreq = 30))\n\n## Not run: \n##D res1 <- update(res1, 100)\n##D # stops at iteration 51\n## End(Not run)\n\n\n## Not run: \n##D res1.1 <- robust(formula = ~1/(1 + exp(-b *(x - a))),\n##D predvars = c(\"x\"), parvars = c(\"a\", \"b\"),\n##D family = binomial(),\n##D lx = -5, ux = 5, prob = rep(1/4, 4),\n##D parset = matrix(c(0.5, 1.5, 0.5, 1.5, 4.0, 4.0, 5.0, 5.0), 4, 2),\n##D x = c(-3, 0, 3),\n##D iter = 150, k =3)\n##D plot(res1.1)\n##D # not optimal\n## End(Not run)\n\n\n###################################\n# user-defined optimality criterion\n##################################\n# When the model is defined by the formula interface\n# A-optimal design for the 2PL model.\n# the criterion function must have argument x, w fimfunc and the parameters defined in 'parvars'.\n# use 'fimfunc' as a function of the design points x, design weights w and\n# the 'parvars' parameters whenever needed.\nAopt <-function(x, w, a, b, fimfunc){\n sum(diag(solve(fimfunc(x = x, w = w, a = a, b = b))))\n}\n## the sensitivtiy function\n# xi_x is a design that put all its mass on x in the definition of the sensitivity function\n# x is a vector of design points\nAopt_sens <- function(xi_x, x, w, a, b, fimfunc){\n fim <- fimfunc(x = x, w = w, a = a, b = b)\n M_inv <- solve(fim)\n M_x <- fimfunc(x = xi_x, w = 1, a = a, b = b)\n sum(diag(M_inv %*% M_x %*% M_inv)) - sum(diag(M_inv))\n}\n\nres2 <- robust(formula = ~1/(1 + exp(-b * (x-a))), predvars = \"x\",\n parvars = c(\"a\", \"b\"), family = \"binomial\",\n lx = -3, ux = 3,\n iter = 1, k = 4,\n crtfunc = Aopt,\n sensfunc = Aopt_sens,\n prob = c(.25, .5, .25),\n parset = matrix(c(-2, 0, 2, 1.25, 1.25, 1.25), 3, 2),\n ICA.control = list(checkfreq = 50, stoptol = .999,\n stop_rule = \"equivalence\",\n rseed = 1))\n## Not run: \n##D res2 <- update(res2, 500)\n## End(Not run)\n\n\n\n\n\n# robust c-optimal design\n# example from Chaloner and Larntz (1989), Figure 3, but robust design\nc_opt <-function(x, w, a, b, fimfunc){\n gam <- log(.95/(1-.95))\n M <- fimfunc(x = x, w = w, a = a, b = b)\n c <- matrix(c(1, -gam * b^(-2)), nrow = 1)\n B <- t(c) %*% c\n sum(diag(B %*% solve(M)))\n}\n\nc_sens <- function(xi_x, x, w, a, b, fimfunc){\n gam <- log(.95/(1-.95))\n M <- fimfunc(x = x, w = w, a = a, b = b)\n M_inv <- solve(M)\n M_x <- fimfunc(x = xi_x, w = 1, a = a, b = b)\n c <- matrix(c(1, -gam * b^(-2)), nrow = 1)\n B <- t(c) %*% c\n sum(diag(B %*% M_inv %*% M_x %*% M_inv)) - sum(diag(B %*% M_inv))\n}\n\n\nres3 <- robust(formula = ~1/(1 + exp(-b * (x-a))), predvars = \"x\",\n parvars = c(\"a\", \"b\"), family = \"binomial\",\n lx = -1, ux = 1,\n parset = matrix(c(0, 7, .2, 6.5), 2, 2, byrow = TRUE),\n prob = c(.5, .5),\n iter = 1, k = 3,\n crtfunc = c_opt, sensfunc = c_sens,\n ICA.control = list(rseed = 1, checkfreq = Inf))\n\n## Not run: \n##D res3 <- update(res3, 300)\n## End(Not run)\n\n\n\n"} {"package":"ICAOD","topic":"sens.bayes.control","snippet":"### Name: sens.bayes.control\n### Title: Returns Control Parameters for Approximating The Integrals In\n### The Bayesian Sensitivity Functions\n### Aliases: sens.bayes.control\n\n### ** Examples\n\nsens.bayes.control()\nsens.bayes.control(cubature = list(maxEval = 50000))\nsens.bayes.control(quadrature = list(level = 4))\n\n\n"} {"package":"ICAOD","topic":"sens.control","snippet":"### Name: sens.control\n### Title: Returns Control Parameters To Find Maximum of The Sensitivity\n### (Derivative) Function Over The Design Space\n### Aliases: sens.control\n\n### ** Examples\n\nsens.control()\nsens.control(optslist = list(maxeval = 1000))\n\n\n"} {"package":"ICAOD","topic":"sens.minimax.control","snippet":"### Name: sens.minimax.control\n### Title: Returns Control Parameters for Verifying General Equivalence\n### Theorem For Minimax Optimal Designs\n### Aliases: sens.minimax.control\n\n### ** Examples\n\nsens.minimax.control()\nsens.minimax.control(n_seg = 4)\n\n\n"} {"package":"ICAOD","topic":"sensbayes","snippet":"### Name: sensbayes\n### Title: Verifying Optimality of Bayesian D-optimal Designs\n### Aliases: sensbayes\n\n### ** Examples\n\n##################################################################\n# Checking the Bayesian D-optimality of a design for the 2Pl model\n##################################################################\nskew2 <- skewnormal(xi = c(0, 1), Omega = matrix(c(1, -0.17, -0.17, .5), nrow = 2),\n alpha = c(-1, 0), lower = c(-3, .1), upper = c(3, 2))\n## Not run: \n##D sensbayes(formula = ~1/(1 + exp(-b *(x - a))),\n##D predvars = \"x\", parvars = c(\"a\", \"b\"),\n##D family = binomial(),\n##D x= c(-2.50914, -1.16780, -0.36904, 1.29227),\n##D w =c(0.35767, 0.11032, 0.15621, 0.37580),\n##D lx = -3, ux = 3,\n##D prior = skew2)\n##D # took 29 seconds on my system!\n## End(Not run)\n\n# It took very long.\n# We re-adjust the tuning parameters in sens.bayes.control to be faster\n# See how we drastically reduce the maxEval and increase the tolerance\n## Not run: \n##D sensbayes(formula = ~1/(1 + exp(-b *(x - a))),\n##D predvars = \"x\", parvars = c(\"a\", \"b\"),\n##D family = binomial(),\n##D x= c(-2.50914, -1.16780, -0.36904, 1.29227),\n##D w =c(0.35767, 0.11032, 0.15621, 0.37580),\n##D lx = -3, ux = 3,prior = skew2,\n##D sens.bayes.control = list(cubature = list(tol = 1e-4, maxEval = 300)))\n##D # took 5 Seconds on my system!\n## End(Not run)\n\n\n\n# Compare it with the following:\nsensbayes(formula = ~1/(1 + exp(-b *(x - a))),\n predvars = \"x\", parvars = c(\"a\", \"b\"),\n family = binomial(),\n x= c(-2.50914, -1.16780, -0.36904, 1.29227),\n w =c(0.35767, 0.11032, 0.15621, 0.37580),\n lx = -3, ux = 3,prior = skew2,\n sens.bayes.control = list(cubature = list(tol = 1e-4, maxEval = 200)))\n# Look at the plot!\n# took 3 seconds on my system\n\n\n########################################################################################\n# Checking the Bayesian D-optimality of a design for the 4-parameter sigmoid emax model\n########################################################################################\nlb <- c(4, 11, 100, 5)\nub <- c(9, 17, 140, 10)\n## Not run: \n##D sensbayes(formula = ~ theta1 + (theta2 - theta1)*(x^theta4)/(x^theta4 + theta3^theta4),\n##D predvars = c(\"x\"), parvars = c(\"theta1\", \"theta2\", \"theta3\", \"theta4\"),\n##D x = c(0.78990, 95.66297, 118.42964,147.55809, 500),\n##D w = c(0.23426, 0.17071, 0.17684, 0.1827, 0.23549),\n##D lx = .001, ux = 500, prior = uniform(lb, ub))\n##D # took 200 seconds on my system\n## End(Not run)\n\n# Re-adjust the tuning parameters to have it faster\n## Not run: \n##D sensbayes(formula = ~ theta1 + (theta2 - theta1)*(x^theta4)/(x^theta4 + theta3^theta4),\n##D predvars = c(\"x\"), parvars = c(\"theta1\", \"theta2\", \"theta3\", \"theta4\"),\n##D x = c(0.78990, 95.66297, 118.42964,147.55809, 500),\n##D w = c(0.23426, 0.17071, 0.17684, 0.1827, 0.23549),\n##D lx = .001, ux = 500, prior = uniform(lb, ub),\n##D sens.bayes.control = list(cubature = list(tol = 1e-3, maxEval = 300)))\n##D # took 4 seconds on my system. See how much it makes difference\n## End(Not run)\n\n## Not run: \n##D # Now we try it with quadrature. Default is 6 nodes\n##D sensbayes(formula = ~ theta1 + (theta2 - theta1)*(x^theta4)/(x^theta4 + theta3^theta4),\n##D predvars = c(\"x\"), parvars = c(\"theta1\", \"theta2\", \"theta3\", \"theta4\"),\n##D x = c(0.78990, 95.66297, 118.42964,147.55809, 500),\n##D w = c(0.23426, 0.17071, 0.17684, 0.1827, 0.23549),\n##D sens.bayes.control = list(method = \"quadrature\"),\n##D lx = .001, ux = 500, prior = uniform(lb, ub))\n##D # 166.519 s\n##D \n##D # use less number of nodes to see if we can reduce the CPU time\n##D sensbayes(formula = ~ theta1 + (theta2 - theta1)*(x^theta4)/(x^theta4 + theta3^theta4),\n##D predvars = c(\"x\"), parvars = c(\"theta1\", \"theta2\", \"theta3\", \"theta4\"),\n##D x = c(0.78990, 95.66297, 118.42964,147.55809, 500),\n##D w = c(0.23426, 0.17071, 0.17684, 0.1827, 0.23549),\n##D sens.bayes.control = list(method = \"quadrature\",\n##D quadrature = list(level = 3)),\n##D lx = .001, ux = 500, prior = uniform(lb, ub))\n##D # we don't have an accurate plot\n##D \n##D # use less number of levels: use 4 nodes\n##D sensbayes(formula = ~ theta1 + (theta2 - theta1)*(x^theta4)/(x^theta4 + theta3^theta4),\n##D predvars = c(\"x\"), parvars = c(\"theta1\", \"theta2\", \"theta3\", \"theta4\"),\n##D x = c(0.78990, 95.66297, 118.42964,147.55809, 500),\n##D w = c(0.23426, 0.17071, 0.17684, 0.1827, 0.23549),\n##D sens.bayes.control = list(method = \"quadrature\",\n##D quadrature = list(level = 4)),\n##D lx = .001, ux = 500, prior = uniform(lb, ub))\n##D \n## End(Not run)\n\n\n"} {"package":"ICAOD","topic":"sensbayescomp","snippet":"### Name: sensbayescomp\n### Title: Verifying Optimality of Bayesian Compound DP-optimal Designs\n### Aliases: sensbayescomp\n\n### ** Examples\n\n##########################################\n# Verifing the DP-optimality of a design\n# The logistic model with two predictors\n##########################################\n\n# The design points and corresponding weights are as follows:\n# Point1 Point2 Point3 Point4 Point5 Point6 Point7\n# 0.07410 -0.31953 -1.00000 1.00000 -1.00000 1.00000 0.30193\n# -1.00000 1.00000 -1.00000 1.00000 -0.08251 -1.00000 1.00000\n# Weight1 Weight2 Weight3 Weight4 Weight5 Weight6 Weight7\n# 0.020 0.275 0.224 0.131 0.092 0.156 0.103\n\n# It should be given to the function as two seperate vectors:\nx1 <- c(0.07409639, -0.3195265, -1, 1, -1, 1, 0.3019317, -1, 1, -1, 1, -0.08251169, -1, 1)\nw1 <- c(0.01992863, 0.2745394, 0.2236575, 0.1312331, 0.09161503, 0.1561454, 0.1028811)\n\n\np <- c(1, -2, 1, -1)\n\n## Not run: \n##D sensbayescomp(formula = ~exp(b0+b1*x1+b2*x2+b3*x1*x2)/(1+exp(b0+b1*x1+b2*x2+b3*x1*x2)),\n##D predvars = c(\"x1\", \"x2\"),\n##D parvars = c(\"b0\", \"b1\", \"b2\", \"b3\"),\n##D family = binomial(),\n##D x = x1, w = w1,\n##D lx = c(-1, -1), ux = c(1, 1),\n##D prior = uniform(p -1.5, p + 1.5),\n##D prob = ~1-1/(1+exp(b0 + b1 * x1 + b2 * x2 + b3 * x1 * x2)),\n##D alpha = .5, plot_3d = \"rgl\",\n##D sens.bayes.control = list(cubature = list(tol = 1e-3, maxEval = 1000)))\n## End(Not run)\n\n\n\n\n\n\n\n"} {"package":"ICAOD","topic":"senslocally","snippet":"### Name: senslocally\n### Title: Verifying Optimality of The Locally D-optimal Designs\n### Aliases: senslocally\n\n### ** Examples\n\n############################\n# Exponential growth model\n############################\n# Verifying optimailty of a locally D-optimal design\nsenslocally(formula = ~a + exp(-b*x),\n predvars = \"x\", parvars = c(\"a\", \"b\"),\n x = c(.1, 1), w = c(.5, .5),\n lx = 0, ux = 1, inipars = c(1, 10))\n\n\n##############################\n# A model with two predictors\n##############################\nx0 <- c(30, 3.861406, 30, 4.600633, 0, 0, 5.111376, 4.168798)\nw0 <- rep(.25, 4)\nsenslocally(formula = ~ V*S/(Km * (1 + I/Kic)+ S * (1 + I/Kiu)),\n predvars = c(\"S\", \"I\"),\n parvars = c(\"V\", \"Km\", \"Kic\", \"Kiu\"),\n x = x0, w = w0,\n lx = c(0, 0), ux = c(30, 60),\n inipars = c(1.5, 5.2, 3.4, 5.6))\n## Not run: \n##D # using package rgl for 3d plot:\n##D res<- senslocally(formula = ~ V*S/(Km * (1 + I/Kic)+ S * (1 + I/Kiu)),\n##D predvars = c(\"S\", \"I\"),\n##D parvars = c(\"V\", \"Km\", \"Kic\", \"Kiu\"),\n##D x = x0, w = w0,\n##D lx = c(0, 0), ux = c(30, 60),\n##D inipars = c(1.5, 5.2, 3.4, 5.6),\n##D plot_3d = \"rgl\")\n##D \n## End(Not run)\n\n###################################\n# user-defined optimality criterion\n##################################\n# When the model is defined by the formula interface\n# Checking the A-optimality for the 2PL model.\n# the criterion function must have argument x, w fimfunc and the parameters defined in 'parvars'.\n# use 'fimfunc' as a function of the design points x, design weights w and\n# the 'parvars' parameters whenever needed.\nAopt <-function(x, w, a, b, fimfunc){\n sum(diag(solve(fimfunc(x = x, w = w, a = a, b = b))))\n}\n## the sensitivtiy function\n# xi_x is a design that put all its mass on x in the definition of the sensitivity function\n# x is a vector of design points\nAopt_sens <- function(xi_x, x, w, a, b, fimfunc){\n fim <- fimfunc(x = x, w = w, a = a, b = b)\n M_inv <- solve(fim)\n M_x <- fimfunc(x = xi_x, w = 1, a = a, b = b)\n sum(diag(M_inv %*% M_x %*% M_inv)) - sum(diag(M_inv))\n}\n\nsenslocally(formula = ~1/(1 + exp(-b * (x-a))), predvars = \"x\",\n parvars = c(\"a\", \"b\"), family = \"binomial\",\n inipars = c(0, 1.5),\n crtfunc = Aopt,\n lx = -2, ux = 2,\n sensfunc = Aopt_sens,\n x = c(-1, 1), w = c(.5, .5))\n# not optimal\n\n\n"} {"package":"ICAOD","topic":"senslocallycomp","snippet":"### Name: senslocallycomp\n### Title: Verifying Optimality of The Locally DP-optimal Designs\n### Aliases: senslocallycomp\n\n### ** Examples\n\n\np <- c(1, -2, 1, -1)\nprior4.4 <- uniform(p -1.5, p + 1.5)\nformula4.4 <- ~exp(b0+b1*x1+b2*x2+b3*x1*x2)/(1+exp(b0+b1*x1+b2*x2+b3*x1*x2))\nprob4.4 <- ~1-1/(1+exp(b0 + b1 * x1 + b2 * x2 + b3 * x1 * x2))\npredvars4.4 <- c(\"x1\", \"x2\")\nparvars4.4 <- c(\"b0\", \"b1\", \"b2\", \"b3\")\nlb <- c(-1, -1)\nub <- c(1, 1)\n\n## That is the optimal design when alpha = .25, see ?locallycomp on how to find it\nxopt <- c(-1, -0.389, 1, 0.802, -1, 1, -1, 1)\nwopt <- c(0.198, 0.618, 0.084, 0.1)\n\n# We want to verfiy the optimality of the optimal design by the general equivalence theorem.\n\nsenslocallycomp(formula = formula4.4, predvars = predvars4.4, parvars = parvars4.4,\n family = binomial(), prob = prob4.4, lx = lb, ux = ub,\n alpha = .25, inipars = p, x = xopt, w = wopt)\n\n## Not run: \n##D # is this design also optimal when alpha = .3\n##D \n##D senslocallycomp(formula = formula4.4, predvars = predvars4.4, parvars = parvars4.4,\n##D family = binomial(), prob = prob4.4, lx = lb, ux = ub,\n##D alpha = .3, inipars = p, x = xopt, w = wopt)\n##D \n##D # when alpha = .3\n##D senslocallycomp(formula = formula4.4, predvars = predvars4.4, parvars = parvars4.4,\n##D family = binomial(), prob = prob4.4, lx = lb, ux = ub,\n##D alpha = .5, inipars = p, x = xopt, w = wopt)\n##D # when alpha = .8\n##D senslocallycomp(formula = formula4.4, predvars = predvars4.4, parvars = parvars4.4,\n##D family = binomial(), prob = prob4.4, lx = lb, ux = ub,\n##D alpha = .8, inipars = p, x = xopt, w = wopt)\n##D \n##D \n##D # when alpha = .9\n##D senslocallycomp(formula = formula4.4, predvars = predvars4.4, parvars = parvars4.4,\n##D family = binomial(), prob = prob4.4, lx = lb, ux = ub,\n##D alpha = .9, inipars = p, x = xopt, w = wopt)\n##D \n##D ## As can be seen, the design looses efficiency as alpha increases.\n## End(Not run)\n\n\n"} {"package":"ICAOD","topic":"sensminimax","snippet":"### Name: sensminimax\n### Title: Verifying Optimality of The Minimax and Standardized maximin\n### D-optimal Designs\n### Aliases: sensminimax\n\n### ** Examples\n\n##########################\n# Power logistic model\n##########################\n# verifying the minimax D-optimality of a design with points x0 and weights w0\nx0 <- c(-4.5515, 0.2130, 2.8075)\nw0 <- c(0.4100, 0.3723, 0.2177)\n# Power logistic model when s = .2\nsensminimax(formula = ~ (1/(1 + exp(-b * (x-a))))^.2,\n predvars = \"x\",\n parvars = c(\"a\", \"b\"),\n family = binomial(),\n x = x0, w = w0,\n lx = -5, ux = 5,\n lp = c(0, 1), up = c(3, 1.5))\n\n##############################\n# A model with two predictors\n##############################\n# Verifying the minimax D-optimality of a design for a model with two predictors\n# The model is the mixed inhibition model.\n# X0 is the vector of four design points that are:\n# (3.4614, 0) (4.2801, 3.1426) (30, 0) (30, 4.0373)\nx0 <- c(3.4614, 4.2801, 30, 30, 0, 3.1426, 0, 4.0373)\nw0 <- rep(1/4, 4)\nsensminimax(formula = ~ V*S/(Km * (1 + I/Kic)+ S * (1 + I/Kiu)),\n predvars = c(\"S\", \"I\"),\n parvars = c(\"V\", \"Km\", \"Kic\", \"Kiu\"),\n family = \"gaussian\",\n x = x0, w = w0,\n lx = c(0, 0), ux = c(30, 60),\n lp = c(1, 4, 2, 4), up = c(1, 5, 3, 5))\n\n##########################################\n# Standardized maximin D-optimal designs\n##########################################\n# Verifying the standardized maximin D-optimality of a design for\n# the loglinear model\n# First we should define the function for 'localdes' argument\n# The function LDOD takes the parameters and returns the points and\n# weights of the locally D-optimal design\nLDOD <- function(theta0, theta1, theta2){\n ## param is the vector of theta = (theta0, theta1, theta2)\n lx <- 0 # lower bound of the design space\n ux <- 150 # upper bound of the design space\n param <- c()\n param[1] <- theta0\n param[2] <- theta1\n param[3] <- theta2\n xstar <- (ux+param[3]) * (lx + param[3]) *\n (log(ux + param[3]) - log(lx + param[3]))/(ux - lx) - param[3]\n return(list(x = c(lx, xstar, ux) , w = rep(1/3, 3)))\n}\nx0 <- c(0, 4.2494, 17.0324, 149.9090)\nw0 <- c(0.3204, 0.1207, 0.2293, 0.3296)\n## Not run: \n##D sensminimax(formula = ~theta0 + theta1* log(x + theta2),\n##D predvars = c(\"x\"),\n##D parvars = c(\"theta0\", \"theta1\", \"theta2\"),\n##D x = x0, w = w0,\n##D lx = 0, ux = 150,\n##D lp = c(2, 2, 1), up = c(2, 2, 15),\n##D localdes = LDOD,\n##D standardized = TRUE,\n##D sens.minimax.control = list(n_seg = 10))\n## End(Not run)\n################################################################\n# Not necessary!\n# The rest of the examples here are only for professional uses.\n################################################################\n# Imagine you have written your own FIM, say in Rcpp that is faster than\n# the FIM created by the formula interface here.\n\n##########################\n# Power logistic model\n##########################\n# For example, th cpp FIM function for the power logistic model is named:\nFIM_power_logistic\nargs(FIM_power_logistic)\n# The arguments do not match the standard of the argument 'fimfunc'\n# in 'sensminimax'\n# So we reparameterize it:\nmyfim1 <- function(x, w, param)\n FIM_power_logistic(x = x, w = w, param =param, s = .2)\n\nargs(myfim1)\n## Not run: \n##D # Verify minimax D-optimality of a design\n##D sensminimax(fimfunc = myfim1,\n##D x = c(-4.5515, 0.2130, 2.8075),\n##D w = c(0.4100, 0.3723, 0.2177),\n##D lx = -5, ux = 5,\n##D lp = c(0, 1), up = c(3, 1.5))\n## End(Not run)\n##############################\n# A model with two predictors\n##############################\n# An example of a model with two-predictors: mixed inhibition model\n# Fisher information matrix:\nFIM_mixed_inhibition\nargs(FIM_mixed_inhibition)\n\n# We should first reparameterize the FIM to match the standard of the\n# argument 'fimfunc'\nmyfim2 <- function(x, w, param){\n npoint <- length(x)/2\n S <- x[1:npoint]\n I <- x[(npoint+1):(npoint*2)]\n out <- FIM_mixed_inhibition(S = S, I = I, w = w, param = param)\n return(out)\n}\nargs(myfim2)\n## Not run: \n##D # Verifyng minimax D-optimality of a design\n##D sensminimax(fimfunc = myfim2,\n##D x = c(3.4614, 4.2801, 30, 30, 0, 3.1426, 0, 4.0373),\n##D w = rep(1/4, 4),\n##D lx = c(0, 0), ux = c(30, 60),\n##D lp = c(1, 4, 2, 4), up = c(1, 5, 3, 5))\n## End(Not run)\n\n#########################################\n# Standardized maximin D-optimal designs\n#########################################\n# An example of a user-written FIM function:\nhelp(FIM_loglin)\n# An example of verfying standardaized maximin D-optimality for a design\n# Look how we re-define the function LDOD above\nLDOD2 <- function(param){\n ## param is the vector of theta = (theta0, theta1, theta2)\n lx <- 0 # lower bound of the design space\n ux <- 150 # upper bound of the design space\n xstar <- (ux + param[3]) * (lx + param[3]) *\n (log(ux + param[3]) - log(lx + param[3]))/(ux - lx) - param[3]\n return(list(x = c(lx, xstar, ux) , w = rep(1/3, 3)))\n}\n\nargs(LDOD2)\n\nsensminimax(fimfunc = FIM_loglin,\n x = x0,\n w = w0,\n lx = 0, ux = 150,\n lp = c(2, 2, 1), up = c(2, 2, 15),\n localdes = LDOD2,\n standardized = TRUE)\n\n\n\n###################################\n# user-defined optimality criterion\n##################################\n# When the model is defined by the formula interface\n# Checking the A-optimality for the 2PL model.\n# the criterion function must have argument x, w fimfunc and the parameters defined in 'parvars'.\n# use 'fimfunc' as a function of the design points x, design weights w and\n# the 'parvars' parameters whenever needed.\nAopt <-function(x, w, a, b, fimfunc){\n sum(diag(solve(fimfunc(x = x, w = w, a = a, b = b))))\n}\n## the sensitivtiy function\n# xi_x is a design that put all its mass on x in the definition of the sensitivity function\n# x is a vector of design points\nAopt_sens <- function(xi_x, x, w, a, b, fimfunc){\n fim <- fimfunc(x = x, w = w, a = a, b = b)\n M_inv <- solve(fim)\n M_x <- fimfunc(x = xi_x, w = 1, a = a, b = b)\n sum(diag(M_inv %*% M_x %*% M_inv)) - sum(diag(M_inv))\n}\n\nsensminimax(formula = ~1/(1 + exp(-b * (x-a))), predvars = \"x\",\n parvars = c(\"a\", \"b\"), family = \"binomial\",\n lp = c(-2, 1), up = c(2, 1.5),\n crtfunc = Aopt,\n lx = -2, ux = 2,\n sensfunc = Aopt_sens,\n x = c(-2, .0033, 2), w = c(.274, .452, .274))\n\n\n\n"} {"package":"ICAOD","topic":"sensmultiple","snippet":"### Name: sensmultiple\n### Title: Verifying Optimality of The Multiple Objective Designs for The\n### 4-Parameter Hill Model\n### Aliases: sensmultiple\n\n### ** Examples\n\n#################################################################\n# Verifying optimality of a design for the 4-parameter Hill model\n#################################################################\n\n## initial estiamtes for the parameters of the Hill model\na <- 0.008949 # ED50\nb <- -1.79 # Hill constant\nc <- 0.137 # lower limit\nd <- 1.7 # upper limit\n# D belongs to c(.001, 1000) ## dose in mg\n## Hill parameters are c(a, b, c, d)\n# dose, minDose and maxDose vector in mg scale\n\nsensmultiple (dose = c(0.001, 0.009426562, 0.01973041, 999.9974),\n w = c(0.4806477, 0.40815, 0.06114173, 0.05006055),\n minDose = .001, maxDose = 1000,\n Hill_par = TRUE,\n inipars = c(a, b, c, d),\n lambda = c(0.05, 0.05, .90),\n delta = -1)\n\n\n\n\n\n## Don't show: \n\n## examples fof using this function for c-optimal designs\n# first row second column: c-optimal design for estimating ED50 of the 4-parameter logistic model\nsensmultiple (dose = c(log(.001), -4.80, log(1000)),\n w = c(.276, .500, .224),\n Hill_par = FALSE,\n minDose = log(.001), maxDose = log(1000),\n inipars = c(d - c, -b, b * log(a), c),\n lambda = c(0, 1, 0),\n delta = -1)\n## criterion value is 1e+24 which will be returned when variance for estimating ED50 is comutationaly negative!\n## if we change the tolerance for finding Moore-Penrose Matrix Inverse to .Machine$double.eps\n# when get 2.201179 for the criterion value\n\n\nsensmultiple (dose = c(-6.907755, -4.664224, -3.925594, 6.907753 ),\n w = c(0.4806477, 0.40815, 0.06114173, 0.05006055 ),\n minDose = log(.001), maxDose = log(1000),\n Hill_par = FALSE,\n inipars = c(d - c, -b, b * log(a), c),\n lambda = c(0.05, 0.05, .90),\n delta = -1)\n## End(Don't show)\n\n\n"} {"package":"ICAOD","topic":"sensrobust","snippet":"### Name: sensrobust\n### Title: Verifying Optimality of The Robust Designs\n### Aliases: sensrobust\n\n### ** Examples\n\n# Verifying a robust design for the two-parameter logistic model\nsensrobust(formula = ~1/(1 + exp(-b *(x - a))),\n predvars = c(\"x\"),\n parvars = c(\"a\", \"b\"),\n family = binomial(),\n prob = rep(1/4, 4),\n parset = matrix(c(0.5, 1.5, 0.5, 1.5, 4.0, 4.0, 5.0, 5.0), 4, 2),\n x = c(0.260, 1, 1.739), w = c(0.275, 0.449, 0.275),\n lx = -5, ux = 5)\n\n\n###################################\n# user-defined optimality criterion\n##################################\n# When the model is defined by the formula interface\n# Checking the A-optimality for the 2PL model.\n# the criterion function must have argument x, w fimfunc and the parameters defined in 'parvars'.\n# use 'fimfunc' as a function of the design points x, design weights w and\n# the 'parvars' parameters whenever needed.\nAopt <-function(x, w, a, b, fimfunc){\n sum(diag(solve(fimfunc(x = x, w = w, a = a, b = b))))\n}\n## the sensitivtiy function\n# xi_x is a design that put all its mass on x in the definition of the sensitivity function\n# x is a vector of design points\nAopt_sens <- function(xi_x, x, w, a, b, fimfunc){\n fim <- fimfunc(x = x, w = w, a = a, b = b)\n M_inv <- solve(fim)\n M_x <- fimfunc(x = xi_x, w = 1, a = a, b = b)\n sum(diag(M_inv %*% M_x %*% M_inv)) - sum(diag(M_inv))\n}\n\nsensrobust(formula = ~1/(1 + exp(-b * (x-a))), predvars = \"x\",\n parvars = c(\"a\", \"b\"), family = \"binomial\",\n crtfunc = Aopt,\n sensfunc = Aopt_sens,\n lx = -3, ux = 3,\n prob = c(.25, .5, .25),\n parset = matrix(c(-2, 0, 2, 1.25, 1.25, 1.25), 3, 2),\n x = c(-2.469, 0, 2.469), w = c(.317, .365, .317))\n# not optimal. the optimal design has four points. see the last example in ?robust\n\n\n"} {"package":"ICAOD","topic":"skewnormal","snippet":"### Name: skewnormal\n### Title: Assumes A Multivariate Skewed Normal Prior Distribution for The\n### Model Parameters\n### Aliases: skewnormal\n\n### ** Examples\n\nskewnormal(xi = c(0, 1),\n Omega = matrix(c(1, -0.17, -0.17, .5), nrow = 2),\n alpha = c(1, 0), lower = c(-3, .1), upper = c(3, 2))\n\n\n"} {"package":"ICAOD","topic":"student","snippet":"### Name: student\n### Title: Multivariate Student's t Prior Distribution for Model Parameters\n### Aliases: student\n\n### ** Examples\n\nskewnormal(xi = c(0, 1),\n Omega = matrix(c(1, -0.17, -0.17, .5), nrow = 2),\n alpha = c(1, 0), lower = c(-3, .1), upper = c(3, 2))\n\n\n"} {"package":"ICAOD","topic":"uniform","snippet":"### Name: uniform\n### Title: Assume A Multivariate Uniform Prior Distribution for The Model\n### Parameters\n### Aliases: uniform\n\n### ** Examples\n\nuniform(lower = c(-3, .1), upper = c(3, 2))\n\n\n"} {"package":"higrad","topic":"higrad","snippet":"### Name: higrad\n### Title: Fitting HiGrad\n### Aliases: higrad\n\n### ** Examples\n\n# fitting linear regression on a simulated dataset\nn <- 1e3\nd <- 10\nsigma <- 0.1\ntheta <- rep(1, d)\nx <- matrix(rnorm(n * d), n, d)\ny <- as.numeric(x %*% theta + rnorm(n, 0, sigma))\nfit <- higrad(x, y, model = \"lm\")\nprint(fit)\n# predict for 10 new samples\nnewx <- matrix(rnorm(10 * d), 10, d)\npred <- predict(fit, newx)\npred\n\n\n\n"} {"package":"ktsolve","topic":"ktsolve","snippet":"### Name: ktsolve\n### Title: Configurable Function for Solving Families of Nonlinear\n### Equations Version: 1.3\n### Aliases: ktsolve\n\n### ** Examples\n\nzfunc<-function(x) {\n\tz<-vector()\nz[1]<- 4*var1 -3*var2 +5*var3\nz[2]<-8*var1 +5*var2 -2*var3\nz\n}\n\n known=list(var2=5)\n guess=list(var1=2,var3=0)\n solv1 <- ktsolve(zfunc,known=known,guess=guess)\n# Successful convergence.\n# solution is:\n# var1 var3 \n# -1.979167 4.583333 \n# \"known\" inputs were:\n# var2\n# known 5 \n eval(solv1$yfunc)(solv1$results$par)\n\n \n known=list(var1=5)\n guess=list(var2=2,var3=0)\n solv2<- ktsolve(zfunc,known=known,guess=guess)\n# Successful convergence.\n# solution is:\n# var2 var3 \n# -12.63158 -11.57895 \n# \"known\" inputs were:\n# var1\n# known 5 \neval(solv2$yfunc)(solv2$results$par)\n \n\n"} {"package":"rabhit","topic":"createFullHaplotype","snippet":"### Name: createFullHaplotype\n### Title: Anchor gene haplotype inference\n### Aliases: createFullHaplotype\n\n### ** Examples\n\n# Load example data and germlines\ndata(samples_db, HVGERM, HDGERM)\n\n# Selecting a single individual\nclip_db = samples_db[samples_db$subject=='I5', ]\n\n# Infering haplotype\nhaplo_db = createFullHaplotype(clip_db,toHap_col=c('v_call','d_call'),\nhapBy_col='j_call',hapBy='IGHJ6',toHap_GERM=c(HVGERM,HDGERM))\n\n\n\n\n"} {"package":"rabhit","topic":"deletionHeatmap","snippet":"### Name: deletionHeatmap\n### Title: Graphical output of single chromosome deletions\n### Aliases: deletionHeatmap\n\n### ** Examples\n\n# Plotting single choromosme deletion from haplotype inference\ndeletionHeatmap(samplesHaplotype)\n\n\n"} {"package":"rabhit","topic":"deletionsByBinom","snippet":"### Name: deletionsByBinom\n### Title: Double chromosome deletion by relative gene usage\n### Aliases: deletionsByBinom\n\n### ** Examples\n\n# Load example data and germlines\ndata(samples_db)\n\n# Selecting a single individual\nclip_db = samples_db[samples_db$subject=='I5', ]\n# Infering haplotype\ndel_binom_df = deletionsByBinom(clip_db)\nhead(del_binom_df)\n\n\n\n"} {"package":"rabhit","topic":"deletionsByVpooled","snippet":"### Name: deletionsByVpooled\n### Title: Single chromosomal D or J gene deletions inferred by the V\n### pooled method\n### Aliases: deletionsByVpooled\n\n### ** Examples\n\n## No test: \ndata(samples_db)\n\n# Infering V pooled deletions\ndel_db <- deletionsByVpooled(samples_db)\nhead(del_db)\n## End(No test)\n\n\n"} {"package":"rabhit","topic":"hapDendo","snippet":"### Name: hapDendo\n### Title: Hierarchical clustering of haplotypes graphical output\n### Aliases: hapDendo\n\n### ** Examples\n\n# Plotting haplotype hierarchical clustering based on the Jaccard distance\n## No test: \nhapDendo(samplesHaplotype)\n## End(No test)\n\n\n\n"} {"package":"rabhit","topic":"hapHeatmap","snippet":"### Name: hapHeatmap\n### Title: Graphical output of alleles division by chromosome\n### Aliases: hapHeatmap\n\n### ** Examples\n\n# Plotting haplotpe heatmap\np <- hapHeatmap(samplesHaplotype)\np$p\n\n\n"} {"package":"rabhit","topic":"nonReliableVGenes","snippet":"### Name: nonReliableVGenes\n### Title: Detect non reliable gene assignment\n### Aliases: nonReliableVGenes\n\n### ** Examples\n\n# Example IGHV call data frame\nclip_db <- data.frame(subject=rep('S1',6),\nv_call=c('IGHV1-69*01','IGHV1-69*01','IGHV1-69*01,IGHV1-69*02',\n'IGHV4-59*01,IGHV4-61*01','IGHV4-59*01,IGHV4-31*02','IGHV4-59*01'))\n# Detect non reliable genes\nnonReliableVGenes(clip_db)\n\n\n"} {"package":"rabhit","topic":"plotDeletionsByBinom","snippet":"### Name: plotDeletionsByBinom\n### Title: Graphical output of double chromosome deletions\n### Aliases: plotDeletionsByBinom\n\n### ** Examples\n\n\n# Load example data and germlines\ndata(samples_db)\n\n# Infering haplotype\ndeletions_db = deletionsByBinom(samples_db);\nplotDeletionsByBinom(deletions_db)\n\n\n\n"} {"package":"rabhit","topic":"plotDeletionsByVpooled","snippet":"### Name: plotDeletionsByVpooled\n### Title: Graphical output for single chromosome D or J gene deletions\n### according to V pooled method\n### Aliases: plotDeletionsByVpooled\n\n### ** Examples\n\n## No test: \n# Load example data and germlines\ndata(samples_db)\ndel_db <- deletionsByVpooled(samples_db)\nplotDeletionsByVpooled(del_db)\n## End(No test)\n\n\n"} {"package":"rabhit","topic":"plotHaplotype","snippet":"### Name: plotHaplotype\n### Title: Graphical output of an inferred haplotype\n### Aliases: plotHaplotype\n\n### ** Examples\n\n\n# Selecting a single individual from the haplotype samples data\nhaplo_db = samplesHaplotype[samplesHaplotype$subject=='I5', ]\n\n# plot haplotype\nplotHaplotype(haplo_db)\n\n\n\n"} {"package":"gJLS2","topic":"gJLS2","snippet":"### Name: gJLS2\n### Title: A Generalized Joint-Location-Scale (gJLS) Test\n### Aliases: gJLS2\n\n### ** Examples\n\nN <- 1000\ngenDAT <- rbinom(N, 2, 0.3)\nsex <- rbinom(N, 1, 0.5)+1\ny <- rnorm(N)\ncovar <- matrix(rnorm(N*10), ncol=10)\n\ngJLS2(GENO=data.frame(\"SNP1\" = genDAT, \"aSNP1\" = genDAT), SEX=sex, Y=y, COVAR=covar)\n\ngJLS2(GENO=genDAT, SEX=sex, Y=y, COVAR=covar , Xchr=TRUE)\n\n\n\n\n"} {"package":"gJLS2","topic":"gJLS2s","snippet":"### Name: gJLS2s\n### Title: generalized Joint-Location-Scale (gJLS) test with summary\n### statistics\n### Aliases: gJLS2s\n\n### ** Examples\n\ngL <- data.frame(\"SNP\" = paste(\"rs\", 1:100, sep=\"\"), \"gL\"=runif(100))\ngS <- runif(100)\n\ngJLS2s(gL = gL, gS=gS)\n\n\n\n\n"} {"package":"gJLS2","topic":"leveneRegA_per_SNP","snippet":"### Name: leveneRegA_per_SNP\n### Title: The generalized Levene's test via a two-stage regression for\n### variance homogeneity by SNP genotype (autosomes)\n### Aliases: leveneRegA_per_SNP\n\n### ** Examples\n\nN <- 100\ngenDAT <- rbinom(N, 2, 0.3)\nY <- rnorm(N)\ncovar <- matrix(rnorm(N*10), ncol=10)\n\n# vanilla example:\nleveneRegA_per_SNP(geno_one=genDAT, Y=Y, COVAR=covar)\n\n# relatedness samples:\nleveneRegA_per_SNP(geno_one=genDAT, Y=Y, COVAR=covar,\nrelated=TRUE)\nleveneRegA_per_SNP(geno_one=genDAT, Y=Y, COVAR=covar,\nrelated=TRUE, clust = factor(rbinom(N, 2, 0.6)))\n\n\n# dosage genotypes example:\nlibrary(\"MCMCpack\")\na <- 0.3\ngeno <- rbinom(N, 2, 0.3)\na <- 0.3 ## uncertainty\ngenPP <- rbind(rdirichlet(sum(geno==0),c(a,(1-a)/2,(1-a)/2)),\n rdirichlet(sum(geno==1),c((1-a)/2,a,(1-a)/2)),\n rdirichlet(sum(geno==2),c((1-a)/2,(1-a)/2,a)))\n\nleveneRegA_per_SNP(geno_one=genPP, Y=Y, COVAR=covar)\nleveneRegA_per_SNP(geno_one=genPP, Y=Y, COVAR=covar,\ngenotypic=TRUE)\n\n# dosage and related samples:\nleveneRegA_per_SNP(geno_one=genPP, Y=Y, COVAR=covar,\nrelated=TRUE, clust = factor(rbinom(N, 1, 0.3)))\nleveneRegA_per_SNP(geno_one=genPP, Y=Y, COVAR=covar,\nrelated=TRUE, clust = factor(rbinom(N, 1, 0.3)), genotypic=TRUE)\n\n\n\n\n\n"} {"package":"gJLS2","topic":"leveneRegX_per_SNP","snippet":"### Name: leveneRegX_per_SNP\n### Title: Levene's regression tests for variance homogeneity by SNP\n### genotype (X-chromosome specific)\n### Aliases: leveneRegX_per_SNP\n\n### ** Examples\n\nN <- 1000\nsex <- rbinom(N, 1, 0.5)+1\nY <- rnorm(N)\ngenDAT <- NA\ngenDAT[sex==2] <- rbinom(sum(sex==2), 2, 0.3)\ntable(genDAT, sex)\ngenDAT[sex==1] <- rbinom(sum(sex==1), 1, 0.3)\ntable(genDAT, sex)\n\nleveneRegX_per_SNP(geno_one=genDAT, SEX=sex, Y=Y)\nleveneRegX_per_SNP(geno_one=genDAT, SEX=sex, Y=Y, genotypic=TRUE)\nleveneRegX_per_SNP(geno_one=genDAT, SEX=sex, Y=Y, loc_alg=\"OLS\")\n\n\n\n"} {"package":"gJLS2","topic":"leveneTests_per_SNP","snippet":"### Name: leveneTests_per_SNP\n### Title: Levene's test for variance homogeneity by SNP genotypes\n### (sex-specific p-values)\n### Aliases: leveneTests_per_SNP\n\n### ** Examples\n\nN <- 5000\nsex <- rbinom(N, 1, 0.5)+1\ngenDAT <- rbinom(N, 2, 0.3)\ny <- rnorm(N);\n\ngenDAT[sex==2] <- rbinom(sum(sex==2), 1, 0.3)\ntable(genDAT, sex)\nleveneTests_per_SNP(geno_one=genDAT, SEX=sex, Y=y^2, transform=TRUE)\n\ngenDAT[sex==2] <- rbinom(sum(sex==2), 1, 0.01)\ntable(genDAT, sex)\nleveneTests_per_SNP(geno_one=genDAT, SEX=sex, Y=y^2, transform=FALSE)\n\nleveneTests_per_SNP(geno_one=rep(0, N), SEX=sex, Y=y^2, transform=TRUE)\nleveneTests_per_SNP(geno_one=rep(0, N), Y=y^2, transform=TRUE)\n\n\n\n\n"} {"package":"gJLS2","topic":"locReg","snippet":"### Name: locReg\n### Title: Location (mean-based association) test\n### Aliases: locReg\n\n### ** Examples\n\nN <- 100\ngenDAT <- rbinom(N, 2, 0.3)\nsex <- rbinom(N, 1, 0.5)+1\ny <- rnorm(N)\nCOVAR <- matrix(rnorm(N*10), ncol=10)\n\nlocReg(GENO=genDAT, SEX=sex, Y=y, COVAR=COVAR)\n\n# correlated example:\nlibrary(\"MASS\")\nyy <- mvrnorm(1, mu= rep(0, N), Sigma = matrix(0.3, N, N) + diag(0.7, N))\nlocReg(GENO=list(\"SNP1\"= genDAT, \"SNP2\" = genDAT[sample(1:100)]),\nSEX=sex, Y=as.numeric(yy), COVAR=COVAR, related = TRUE,\nclust = rep(1, 100))\n\n# sibpair example:\npairedY <- mvrnorm(N/2,rep(0,2),matrix(c(1,0.2,0.2,1), 2))\nyy <- c(pairedY[,1], pairedY[,2])\nlocReg(GENO=list(\"SNP1\"= genDAT, \"SNP2\" = genDAT[sample(1:100)]),\nSEX=sex, Y=as.numeric(yy), COVAR=COVAR, related = TRUE,\nclust = rep(c(1:50), 2))\n\n# Xchr data example:\ngenDAT1 <- rep(NA, N)\ngenDAT1[sex==1] <- rbinom(sum(sex==1), 1, 0.5)\ngenDAT1[sex==2] <-rbinom(sum(sex==2), 2, 0.5)\nlocReg(GENO=genDAT1, SEX=sex, Y=y, COVAR=COVAR, Xchr=TRUE)\n\n\n\n\n"} {"package":"gJLS2","topic":"scaleReg","snippet":"### Name: scaleReg\n### Title: Scale (variance-based association) test\n### Aliases: scaleReg\n\n### ** Examples\n\nN <- 1000\ngenoDAT <- rbinom(N, 2, 0.3)\nsex <- rbinom(N, 1, 0.5)+1\nY <- rnorm(N)\ncovar <- matrix(rnorm(N*10), ncol=10)\n\n# vanilla example:\nscaleReg(GENO=list(genoDAT, genoDAT), Y=Y, COVAR=covar)\nscaleReg(GENO=list(genoDAT, genoDAT), Y=Y, COVAR=covar, genotypic=TRUE)\nscaleReg(GENO=list(genoDAT, genoDAT), Y=Y, COVAR=covar, origLev = TRUE)\nscaleReg(GENO=list(genoDAT, genoDAT), Y=Y, COVAR=covar, origLev = TRUE, SEX=sex)\n\n\n\n"} {"package":"swipeR","topic":"swipeR","snippet":"### Name: swipeR\n### Title: HTML widget displaying a carousel\n### Aliases: swipeR\n\n### ** Examples\n\nlibrary(swipeR)\nlibrary(htmltools)\n\nwrapper <- swipeRwrapper(\n tags$img(src = \"https://swiperjs.com/demos/images/nature-1.jpg\"),\n tags$img(src = \"https://swiperjs.com/demos/images/nature-2.jpg\"),\n tags$img(src = \"https://swiperjs.com/demos/images/nature-3.jpg\"),\n tags$img(src = \"https://swiperjs.com/demos/images/nature-4.jpg\"),\n tags$img(src = \"https://swiperjs.com/demos/images/nature-5.jpg\"),\n tags$img(src = \"https://swiperjs.com/demos/images/nature-6.jpg\"),\n tags$img(src = \"https://swiperjs.com/demos/images/nature-7.jpg\"),\n tags$img(src = \"https://swiperjs.com/demos/images/nature-8.jpg\")\n)\n\nswipeR(\n wrapper, height = \"400px\", width = \"70%\", thumbs = TRUE, keyboard = TRUE,\n on = list(reachEnd = htmlwidgets::JS(\"function() {alert('the end');}\"))\n)\n\n# Shiny example ####\nlibrary(swipeR)\nlibrary(shiny)\nlibrary(ggplot2)\n\nwrapper <- swipeRwrapper(\n div(\n plotOutput(\"ggplot1\", width = \"500px\", height = \"400px\"),\n align = \"center\"\n ),\n div(\n plotOutput(\"ggplot2\", width = \"500px\", height = \"400px\"),\n align = \"center\"\n ),\n div(\n plotOutput(\"ggplot3\", width = \"500px\", height = \"400px\"),\n align = \"center\"\n ),\n div(\n plotOutput(\"ggplot4\", width = \"500px\", height = \"400px\"),\n align = \"center\"\n )\n)\n\nui <- fluidPage(\n tags$head(\n tags$style(HTML(\n \".shiny-plot-output {border: 2px solid royalblue;}\"\n ))\n ),\n br(),\n fluidRow(\n column(\n 12,\n swipeR(\n wrapper, height = \"450px\", width = \"80%\", effect = \"cube\", speed = 2000,\n navigationColor = \"black\", rewind = TRUE, id = \"CAROUSEL\"\n )\n ),\n column(\n 12,\n br(), br(), br(),\n ),\n column(\n 3, align = \"center\",\n actionButton(\n \"btn1\", \"Scatter plot\", class = \"btn-primary\",\n onclick = \"document.getElementById('CAROUSEL').swiper.slideTo(0);\"\n )\n ),\n column(\n 3, align = \"center\",\n actionButton(\n \"btn2\", \"Line chart\", class = \"btn-primary\",\n onclick = \"document.getElementById('CAROUSEL').swiper.slideTo(1);\"\n )\n ),\n column(\n 3, align = \"center\",\n actionButton(\n \"btn3\", \"Bar chart\", class = \"btn-primary\",\n onclick = \"document.getElementById('CAROUSEL').swiper.slideTo(2);\"\n )\n ),\n column(\n 3, align = \"center\",\n actionButton(\n \"btn4\", \"Boxplots\", class = \"btn-primary\",\n onclick = \"document.getElementById('CAROUSEL').swiper.slideTo(3);\"\n )\n )\n )\n)\n\nserver <- function(input, output, session) {\n output[[\"ggplot1\"]] <- renderPlot({\n ggplot(mtcars, aes(wt, mpg)) + geom_point() +\n theme(panel.border = element_rect(fill = NA, color = \"firebrick\"))\n }, width = 500, height = 400)\n output[[\"ggplot2\"]] <- renderPlot({\n ggplot(economics, aes(date, unemploy)) + geom_line()\n }, width = 500, height = 400)\n output[[\"ggplot3\"]] <- renderPlot({\n ggplot(mpg, aes(class)) + geom_bar()\n }, width = 500, height = 400)\n output[[\"ggplot4\"]] <- renderPlot({\n ggplot(mpg, aes(class, hwy)) + geom_boxplot()\n }, width = 500, height = 400)\n}\n\nif(interactive()) shinyApp(ui, server)\n\n\n# other Shiny example ####\nlibrary(swipeR)\nlibrary(shiny)\nlibrary(shinyWidgets)\nlibrary(ggplot2)\nlibrary(ggthemes)\n\nwrapper <- swipeRwrapper(\n div(\n fluidRow(\n column(\n 6,\n awesomeRadio(\n \"theme\", \"Choose a theme\",\n c(\n \"Calc\",\n \"Clean\",\n \"Economist\",\n \"Excel\",\n \"FiveThirtyEight\",\n \"Foundation\",\n \"Google Docs\",\n \"Highcharts\",\n \"Pander\",\n \"Solarized\",\n \"Stata\",\n \"Wall Street\"\n )\n )\n ),\n column(\n 6,\n tags$p(\"The Shiny slider does not work here...\"),\n tags$label(\"Base font size\"),\n tags$input(\n type = \"range\", min = \"10\", max = \"20\", value = \"12\",\n oninput =\n \"this.nextElementSibling.value = this.value;\n Shiny.setInputValue('slider', this.value);\"\n ),\n tags$output(\"12\", style = \"font-weight: bold; color: blue\"),\n br(), hr(), br(),\n materialSwitch(\"facets\", \"Facets?\", status = \"info\"),\n conditionalPanel(\n condition = \"input.facets\",\n awesomeRadio(\n \"direction\", label = NULL, status = \"info\",\n choices = c(\"by row\" = \"row\", \"by column\" = \"column\"),\n )\n ),\n br(), hr(), br(),\n actionButton(\n \"btn\", \"Add slide\", class = \"btn-primary btn-block\",\n onclick = \"document.getElementById('SWIPER').swiper.appendSlide(\n '