Sim_Fn: Function to simulate data sets for testing SDFA model

Usage Arguments Examples

View source: R/Sim_Fn.R

Usage

1
Sim_Fn(n_species, n_years, n_stations = 20, phi = NULL, n_factors = 2, SpatialScale = 0.1, SD_O = 0.5, SD_E = 0.2, SD_extra = 0.1, rho = 0.8, logMeanDens = 1, Lmat = NULL, Loc = NULL, RandomSeed = NA)

Arguments

n_species
n_years
n_stations
phi
n_factors
SpatialScale
SD_O
SD_E
SD_extra
rho
logMeanDens
Lmat
Loc
RandomSeed

Examples

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
##---- Should be DIRECTLY executable !! ----
##-- ==>  Define data, use random,
##--	or do  help(data=index)  for the standard data sets.

## The function is currently defined as
function (n_species, n_years, n_stations = 20, phi = NULL, n_factors = 2, 
    SpatialScale = 0.1, SD_O = 0.5, SD_E = 0.2, SD_extra = 0.1, 
    rho = 0.8, logMeanDens = 1, Lmat = NULL, Loc = NULL, RandomSeed = NA) 
{
    if (!is.na(RandomSeed)) 
        set.seed(RandomSeed)
    if (is.null(Lmat)) {
        Lmat = matrix(rnorm(n_factors * n_species), nrow = n_species, 
            ncol = n_factors)
        for (i in 1:ncol(Lmat)) {
            Lmat[seq(from = 1, to = i - 1, length = i - 1), i] = 0
            if (Lmat[, i][which.max(abs(Lmat[, i]))] < 0) {
                Lmat[, i] = -1 * Lmat[, i]
            }
        }
    }
    if (is.null(phi)) 
        phi = rnorm(n_factors, mean = 0, sd = 1)
    Beta = rep(logMeanDens, n_species)
    if (is.null(Loc)) 
        Loc = cbind(x = runif(n_stations, min = 0, max = 1), 
            y = runif(n_stations, min = 0, max = 1))
    model_O <- RMgauss(var = SD_O^2, scale = SpatialScale)
    model_E <- RMgauss(var = SD_E^2, scale = SpatialScale)
    Omega = matrix(NA, nrow = n_stations, ncol = n_factors)
    for (i in 1:n_factors) {
        Omega[, i] = RFsimulate(model = model_O, x = Loc[, "x"], 
            y = Loc[, "y"])@data[, 1]
    }
    Epsilon = array(NA, dim = c(n_stations, n_factors, n_years))
    for (i in 1:n_factors) {
        Epsilon[, i, 1] = RFsimulate(model = model_E, x = Loc[, 
            "x"], y = Loc[, "y"])@data[, 1]
        for (t in 2:n_years) {
            Epsilon[, i, t] = rho * Epsilon[, i, t - 1] + RFsimulate(model = model_E, 
                x = Loc[, "x"], y = Loc[, "y"])@data[, 1]
        }
    }
    Psi = array(NA, dim = c(n_stations, n_factors, n_years))
    for (i in 1:n_factors) {
        for (t in 1:n_years) {
            Psi[, i, t] = phi[i] * rho^t + Epsilon[, i, t] + 
                Omega[, i]/(1 - rho)
        }
    }
    Theta = array(NA, dim = c(n_stations, n_species, n_years))
    for (s in 1:n_stations) {
        for (t in 1:n_years) {
            Theta[s, , t] = Lmat %*% Psi[s, , t]
        }
    }
    DF = NULL
    for (s in 1:n_stations) {
        for (p in 1:n_species) {
            for (t in 1:n_years) {
                Tmp = c(sitenum = s, spp = p, year = t, catch = rpois(1, 
                  lambda = exp(Theta[s, p, t] + logMeanDens + 
                    SD_extra * rnorm(1))), waterTmpC = 0)
                DF = rbind(DF, Tmp)
            }
        }
    }
    DF = data.frame(DF, row.names = NULL)
    DF[, "spp"] = factor(letters[DF[, "spp"]])
    if (n_species > 26) 
        stop("problem with using letters")
    Sim_List = list(DF = DF, Psi = Psi, Lmat = Lmat, phi = phi, 
        Loc = Loc, Omega = Omega, Epsilon = Epsilon, Theta = Theta, 
        Psi = Psi)
    return(Sim_List)
  }

James-Thorson/spatial_DFA documentation built on July 9, 2020, 7:56 a.m.