Can't work out what is missing from my R function code - stopping it from running well

Question

I'm trying to do Bayesian Occupancy Analysis with site covariates. My first step is making a function. I keep on getting the + in my R console indicating it thinks my code is incomplete. Having ran lines individually I am pretty certain the issue lies in the first line of code. However, I can't work out where exactly I've missed something out and hence where the problem originally lies.

data.fn <- function(R = 39, T = 14, xmin= 0, xmax= 1, alpha.psi = 0.4567,
                beta.psi = 0.0338, alpha.p = 0.4, beta.p = 0.4) {
y <- array(dim = c(R,T)) #This creates an array for counts
  #Ecological Process
  #Covariate values
X <- sort(runif(n=R, min = xmin, max = xmax))
#Expected occurence-covariate relationship
 psi <- plogis(alpha.psi + beta.psi *X) #this applies the inverse logit

 #Add Bernoulli Noise - drawing indicator of occurence (z) from bernoulli psi 
 z<- rbinom(n = R, size = 1, prob = psi)
 occ.fs <- sum(z) #"Finite Sample Occupancy"
 "Make a census"
  p.eff <- z*p
  for (i in 1:T) {
    y[,i] <- rbinom(n=R, size = 1, prob = p.eff)
  }
}

There's more code - i.e. the {} function is complete but the issue started before that was ran and I keep having issues uploading the code into Stack. The error message is simply + all down the left hand side of the R console

EDIT

Could there be something wrong with how R is sensing stuff? For instance with the following code

naive.pred <- plogis(predict(glm(apply(y, 1, max) ~ X + I (X^2),
                             family = binomial)))

I got the error message - unexpected symbol (the bracket) in the family = binomial yet each bracket is paired correctly- there's no extra unnecessary brackets?

Answer 1

While I did not see a + issue when I looked over your code, you are not simulating the observed data correctly and there was a p object inside your function that did not have an argument passed to it. You did create a logit linear predictor for psi using alpha.psi and beta.psi, however you are lacking a logit linear predictor for the probability of detecting a species given that they are there using alpha.p and beta.p. Assuming that the covariate X is used for both the latent occupancy state and the observation model the code should become.

data.fn <- function(R = 39, T = 14, xmin= 0, xmax= 1, alpha.psi = 0.4567,
                                        beta.psi = 0.0338, alpha.p = 0.4, beta.p = 0.4) {
    y <- array(dim = c(R,T)) #This creates an array for counts
    #Ecological Process
    #Covariate values
    X <- sort(runif(n=R, min = xmin, max = xmax))
    #Expected occurence-covariate relationship
    psi <- plogis(alpha.psi + beta.psi *X) #this applies the inverse logit

    #Add Bernoulli Noise - drawing indicator of occurence (z) from bernoulli psi 
    z<- rbinom(n = R, size = 1, prob = psi)
    occ.fs <- sum(z) #"Finite Sample Occupancy"
    # Linear predictor for detection,
    #  assuming the same covariate is used for detection
    p.eff <- plogis(alpha.p + beta.p * X)
    for (i in 1:T) {
        y[,i] <- rbinom(n=R, size = 1, prob = p.eff * z)
    }
    return(list(y = y, z = z, X = X, occ.fs = occ.fs))
}

This code assumes that you are passing logit scale parameters to the data, so if you are trying to simulate data such that X has a very marginal and positive influence on occupancy then you are off to the races, so to speak. If you are looking for a more pronounced effect, than you should increase the effect size. Finally, 39 site is very few for an occupancy analysis given that binary deteciton / non-detection data is quite information poor. Don't be surprised if you posterior estimates you get from analyzing the dataset do not actually return the parameters used to simulate the data.