timcdlucas / INLAutils
Showing 3 of 17 files from the diff.

@@ -17,7 +17,7 @@
Loading
17 17
#' @param modform Formula or list of Formula for the model or models, respectively in INLA format
18 18
#' @param family Character string or list of Character string giving the family or families, respectively of the response in INLA format 
19 19
#' @param mesh inla.mesh object, consisting in triangle mesh to be defined using INLA.
20 -
#' @param ntrials Numeric value (1,2,...) setting the number of trials for a Binomial family
20 +
#' @param ntrials vector of positive integers c(1,2,...) setting the number of trials for a Binomial family
21 21
#' @param int.strategy Character string giving INLA integration strategy
22 22
#' @param alpha Numeric value (0,...,1) giving the threshold for computing confidence intervals (1-alpha) of rmse and mae estimation
23 23
#' @param mae If TRUE, compute the mean absolute error (mae) and the root mean square error (rmse). If FALSE, compute the rmse only.
@@ -25,6 +25,7 @@
Loading
25 25
#' @param sqroot If TRUE, compute the square root of the observed and predicted values to generate the rmse and/or mae. 
26 26
#' If FALSE, the rmse and/or mae are computed without transformation of the data.
27 27
#' @param print Logical to determine whether to print a summary of the results.
28 +
#' @param plot Logical to determine whether to print the plot of the results. 
28 29
#' @param ... other arguments passed to inla
29 30
#'
30 31
#' @export
@@ -65,7 +66,7 @@
Loading
65 66
#' dataframe$y <- round(runif(length(dataframe$y), 1, 12))#create positive discrete response
66 67
#' modform <- stats::as.formula(paste('y ~ -1+ y.intercept + x1 + f(spatial.field, model=spde)'))
67 68
#' family <- list('gamma')#one model
68 -
#' ntrials <- round(max(dataframe$y,na.rm=TRUE))# for Binomial family
69 +
#' ntrials <- rep(round(max(dataframe$y,na.rm=TRUE)*2),length(dataframe$y))# create ntrials for Binomial family
69 70
#' alpha <- 0.05#rmse and mae confidence intervals (1-alpha)
70 71
#'
71 72
#' # run the function
@@ -74,228 +75,233 @@
Loading
74 75
#'
75 76
#' # SLOO function with one model formula (Binomial) and linear terms
76 77
#' sloo1<-inlasloo(dataframe=dataframe, long='long', lat='lat', y='y', ss=1, rad=0.1,
77 -
#' ntrials=5,modform='y ~ -1+ x1 + x2 + y.intercept + f(spatial.field, model=spde)',mesh=mesh,
78 +
#' ntrials=ntrials,modform='y ~ -1+ x1 + x2 + y.intercept + f(spatial.field, model=spde)',mesh=mesh,
78 79
#' family='binomial',alpha=0.05,mae=FALSE,ds=FALSE,sqroot=FALSE)
79 80
#' 
80 81
#' # SLOO function with two families (Binomial and Poisson) and linear terms
81 82
#' sloo2<-inlasloo(dataframe=dataframe, long='long', lat='lat', y='y', ss=1, rad=0.1,
82 -
#' ntrials=5,modform='y ~ -1+ x1 + x2 + y.intercept + f(spatial.field, model=spde)',mesh=mesh,
83 +
#' ntrials=ntrials,modform='y ~ -1+ x1 + x2 + y.intercept + f(spatial.field, model=spde)',mesh=mesh,
83 84
#' family=list('binomial','poisson'),alpha=0.05,mae=TRUE,ds=FALSE,sqroot=FALSE)
84 85
#' 
85 86
#' # SLOO function with one family (Binomial) and two model formulas
86 87
#' sloo3<-inlasloo(dataframe=dataframe, long='long', lat='lat', y='y', ss=1,
87 -
#' rad=0.1, ntrials=5,modform=list('y ~ -1+ y.intercept + f(spatial.field, model=spde)',
88 +
#' rad=0.1, ntrials=ntrials,modform=list('y ~ -1+ y.intercept + f(spatial.field, model=spde)',
88 89
#' 'y ~ -1+ x1 + x2 + y.intercept + f(spatial.field, model=spde)'),
89 90
#' mesh=mesh,family='binomial',alpha=0.05,mae=TRUE,ds=FALSE,sqroot=FALSE)
90 91
#'  
91 92
#' # SLOO function with two families (Binomial and Poisson), two model formulas, 
92 93
#' # and MAE, DS, and sqroot activated
93 94
#' sloo4<-inlasloo(dataframe=dataframe, long='long', lat='lat', y='y', ss=1,
94 -
#' rad=0.1, ntrials=5,modform=list('y ~ -1+ y.intercept + f(spatial.field, model=spde)',
95 +
#' rad=0.1, ntrials=ntrials,modform=list('y ~ -1+ y.intercept + f(spatial.field, model=spde)',
95 96
#' 'y ~ -1+ x1 + x2 + y.intercept + f(spatial.field, model=spde)'),
96 97
#' mesh=mesh,list('binomial','poisson'),alpha=0.05,mae=TRUE,ds=TRUE,sqroot=TRUE)
97 98
#' }
98 99
99 -
 
100 +
100 101
inlasloo <- function(dataframe, long, lat, y, ss, rad, modform, 
101 -
                         family, mesh, ntrials = NULL, int.strategy = "eb", alpha = 0.05,
102 -
                         mae = FALSE, ds = FALSE, sqroot = FALSE, print = FALSE, ...) {
103 -
    message("Identification of input parameters values")
104 -
    message("#########################################", "\n")
105 -
    if (class(modform) != "list") {
106 -
        modform <- list(modform)  #if there is only one model make it as list
107 -
    } else {
108 -
        modform <- modform  # keep it as list for multiple models
109 -
    }
110 -
    if (class(family) != "list") {
111 -
        family <- list(family)  #if there is only one model make it as list
112 -
    } else {
113 -
        family <- family  # keep it as list for multiple models
114 -
    }
115 -
    famnames <- (paste(family, collapse = ","))
116 -
    message(paste("number of models =", length(modform) * length(family)))
117 -
    message(ifelse(missing(dataframe), "dataframe not specified", paste("number of rows in dataframe =", nrow(dataframe))))
118 -
    message(ifelse(missing(long), "long not specified", paste("longitude =", long)))
119 -
    message(ifelse(missing(lat), "lat not specified", paste("latitude =", lat)))
120 -
    message(ifelse(missing(y), "y not specified", paste("response =", y)))
121 -
    message(ifelse(missing(ss), "ss not specified", paste("sampling size =", ss)))
122 -
    message(ifelse(missing(rad), "rad not specified", paste("radius of disc of removed observations =", round(rad, 2))))
123 -
    message(ifelse(mae == TRUE, "RMSE and MAE computed", "RMSE computed"))
124 -
    message(ifelse(ds == TRUE, "DS computed", "DS not computed"))
125 -
    message(ifelse(sqroot == TRUE, "square root for RMSE and MAE computed", "square root for RMSE and MAE not computed"))
126 -
    message(ifelse(missing(family), "family not specified", paste("INLA family distribution of response =", famnames)))
127 -
    message(ifelse(is.null(ntrials) == "TRUE", "family is not Binomial so ntrials is not specified", ifelse(ntrials > 0, paste("number of Bernoulli trials =", 
128 -
        ntrials), "ntrials should be positive integer or not specified")))
129 -
    message(ifelse(missing(mesh), "mesh not specified", paste("number of mesh vertices =", mesh$n)))
130 -
    message(ifelse(int.strategy == "eb", "INLA integration strategy =  empirical bayes ", paste(int.strategy, ": user-defined INLA integration strategy")))
131 -
    message(ifelse(alpha != 0.05, paste("(1-alpha)", "", "% credible intervals of scores"), "default 95% credible intervals of scores"))
132 -
    message("End identification of input parameters values", "\n")
133 -
    message("#############################################", "\n")
134 -
    
135 -
    if (class(dataframe) != "data.frame") 
136 -
        stop("unexpected class instead of data.frame class for the argument \"dataframe\"")
137 -
    if (class(long) != "character") 
138 -
        stop("unexpected class instead of character class for the argument \"long\"")
139 -
    if (class(lat) != "character") 
140 -
        stop("unexpected class instead of character class for the argument \"lat\"")
141 -
    if (class(y) != "character") 
142 -
        stop("unexpected class instead of character class for the argument \"y\"")
143 -
    if (class(mesh) != "inla.mesh") 
144 -
        stop("unexpected class instead of inla.mesh class for the argument \"mesh\"")
145 -
    if (ss > nrow(dataframe)) 
146 -
        stop("sample size larger than number of observations in the dataframe")
147 -
    if (ss < 1) 
148 -
        stop("sample size smaller than 1")
149 -
    if (any(is.na(dataframe[, long]) == TRUE)) 
150 -
        stop("NA present in long")
151 -
    if (any(is.na(dataframe[, lat]) == TRUE)) 
152 -
        stop("NA present in lat")
153 -
    if (rad >= min(max(dataframe[, long]) - min(dataframe[, long]), max(dataframe[, lat]) - min(dataframe[, lat]))) 
154 -
        stop("radius (rad) equal or larger than geographic extent of the observed data")
155 -
    if (rad < 0) 
156 -
        stop("only non-negative values for the radius (rad) can be used")
157 -
    # end tests
102 +
                     family, mesh, ntrials = NULL, int.strategy = "eb", alpha = 0.05,
103 +
                     mae = FALSE, ds = FALSE, sqroot = FALSE, print = FALSE, plot= TRUE, ...) {
104 +
  message("Identification of input parameters values")
105 +
  message("#########################################", "\n")
106 +
  if (class(modform) != "list") {
107 +
    modform <- list(modform)  #if there is only one model make it as list
108 +
  } else {
109 +
    modform <- modform  # keep it as list for multiple models
110 +
  }
111 +
  if (class(family) != "list") {
112 +
    family <- list(family)  #if there is only one model make it as list
113 +
  } else {
114 +
    family <- family  # keep it as list for multiple models
115 +
  }
116 +
  famnames <- (paste(family, collapse = ","))
117 +
  message(paste("number of models =", length(modform) * length(family)))
118 +
  message(ifelse(missing(dataframe), "dataframe not specified", paste("number of rows in dataframe =", nrow(dataframe))))
119 +
  message(ifelse(missing(long), "long not specified", paste("longitude =", long)))
120 +
  message(ifelse(missing(lat), "lat not specified", paste("latitude =", lat)))
121 +
  message(ifelse(missing(y), "y not specified", paste("response =", y)))
122 +
  message(ifelse(missing(ss), "ss not specified", paste("sampling size =", ss)))
123 +
  message(ifelse(missing(rad), "rad not specified", paste("radius of disc of removed observations =", round(rad, 2))))
124 +
  message(ifelse(mae == TRUE, "RMSE and MAE computed", "RMSE computed"))
125 +
  message(ifelse(ds == TRUE, "DS computed", "DS not computed"))
126 +
  message(ifelse(sqroot == TRUE, "square root for RMSE and MAE computed", "square root for RMSE and MAE not computed"))
127 +
  message(ifelse(missing(family), "family not specified", paste("INLA family distribution of response =", famnames)))
128 +
  message(ifelse((family=="binomial"), 
129 +
                 ifelse(isTRUE(all(ntrials>0) & all(abs(ntrials - round(ntrials)) < .Machine$double.eps^0.5))==TRUE,
130 +
                        "number of trials provided for binomial family", "ntrials should be a vector of positive integers"),
131 +
                 ""))
132 +
  message(ifelse(missing(mesh), "mesh not specified", paste("number of mesh vertices =", mesh$n)))
133 +
  message(ifelse(int.strategy == "eb", "INLA integration strategy =  empirical bayes ", paste(int.strategy, ": user-defined INLA integration strategy")))
134 +
  message(ifelse(alpha != 0.05, paste("(1-alpha)", "", "% credible intervals of scores"), "default 95% credible intervals of scores"))
135 +
  message("End identification of input parameters values", "\n")
136 +
  message("#############################################", "\n")
137 +
  
138 +
  if (class(dataframe) != "data.frame") 
139 +
    stop("unexpected class instead of data.frame class for the argument \"dataframe\"")
140 +
  if (class(long) != "character") 
141 +
    stop("unexpected class instead of character class for the argument \"long\"")
142 +
  if (class(lat) != "character") 
143 +
    stop("unexpected class instead of character class for the argument \"lat\"")
144 +
  if (class(y) != "character") 
145 +
    stop("unexpected class instead of character class for the argument \"y\"")
146 +
  if (class(mesh) != "inla.mesh") 
147 +
    stop("unexpected class instead of inla.mesh class for the argument \"mesh\"")
148 +
  if (ss > nrow(dataframe)) 
149 +
    stop("sample size larger than number of observations in the dataframe")
150 +
  if (ss < 1) 
151 +
    stop("sample size smaller than 1")
152 +
  if (any(is.na(dataframe[, long]) == TRUE)) 
153 +
    stop("NA present in long")
154 +
  if (any(is.na(dataframe[, lat]) == TRUE)) 
155 +
    stop("NA present in lat")
156 +
  if (rad >= min(max(dataframe[, long]) - min(dataframe[, long]), max(dataframe[, lat]) - min(dataframe[, lat]))) 
157 +
    stop("radius (rad) equal or larger than geographic extent of the observed data")
158 +
  if (rad < 0) 
159 +
    stop("only non-negative values for the radius (rad) can be used")
160 +
  # end tests
161 +
  
162 +
  # start of the function
163 +
  # set values based on user input or default
164 +
  dataframe <- dataframe
165 +
  rad <- rad
166 +
  ss <- ss
167 +
  int.strategy <- int.strategy
168 +
  
169 +
  # data tranformation for further computing
170 +
  
171 +
  # Stupid way of getting correct column names without collisions (almost certainly).
172 +
  colnames(dataframe)[colnames(dataframe) == y] <- "tmp_y_column_name_eeeeee"
173 +
  colnames(dataframe)[colnames(dataframe) == long] <- "tmp_long_column_name_eeeeee"
174 +
  colnames(dataframe)[colnames(dataframe) == lat] <- "tmp_lat_column_name_eeeeee"
175 +
  
176 +
  colnames(dataframe)[colnames(dataframe) == "tmp_y_column_name_eeeeee"] <- "y"
177 +
  colnames(dataframe)[colnames(dataframe) == "tmp_long_column_name_eeeeee"] <- "long"
178 +
  colnames(dataframe)[colnames(dataframe) == "tmp_lat_column_name_eeeeee"] <- "lat"
179 +
  
180 +
  measurevar <- "y"  #replace response by 'y' in model formula
181 +
  for (i in 1:length(modform)) {
182 +
    modform[i] <- gsub(".*~", "", modform[i])
183 +
    modform[i] <- paste(measurevar, paste(modform[i]), sep = " ~ ")
184 +
  }
185 +
  
186 +
  # start of SLOO process
187 +
  test_rows <- sample(nrow(dataframe), ss)
188 +
  test.df <- dataframe[test_rows, ]  #sample subset of size ss from the dataframe
189 +
  test.ntrials <- ntrials[test_rows]
190 +
  slooresults <- list()
191 +
  slooresultsprint<- list()
192 +
  familys <- family
193 +
  modforms <- modform
194 +
  result.list <- expand.grid(familys, modforms)  #create all possible combinations of family and modforms
195 +
  result.list <- lapply(apply(result.list, 1, identity), unlist)  #put into a list
196 +
  for (j in 1:(length(result.list))) {
197 +
    family <- as.character(result.list[[j]][1])  #select family for j element
198 +
    modform <- stats::as.formula(as.character(result.list[[j]][2]))  #select modform for j element
199 +
    spde <- INLA::inla.spde2.matern(mesh, alpha = 2)  #SPDE computation
158 200
    
159 -
   # start of the function
160 -
   # set values based on user input or default
161 -
    dataframe <- dataframe
162 -
    rad <- rad
163 -
    ss <- ss
164 -
    int.strategy <- int.strategy
165 -
    
166 -
    # data tranformation for further computing
201 +
    for (i in 1:nrow(test.df)) {
202 +
      if (i == 1) {
203 +
        p <- c()
204 +
      }
205 +
      # Training data (test point & buffer removed)
206 +
      train <- dataframe[sqrt((dataframe$long - test.df[i, ]$long)^2 + (dataframe$lat - test.df[i, ]$lat)^2) > rad, ]
207 +
      # predictions in test locations
208 +
      predptdf <- test.df[i, ]  #predict on test point (one point only)
209 +
      A.pred <- INLA::inla.spde.make.A(mesh = mesh, loc = cbind(test.df[i, ]$long, test.df[i, ]$lat))
210 +
      stk.pred <- INLA::inla.stack(data = list(y = NA), A = list(A.pred, 1, 1, 1), effects = list(spatial.field = 1:spde$n.spde, 
211 +
                                                                                                  iid = 1:length(predptdf$y), y.intercept = rep(1, length(predptdf$y)), covariate = subset(predptdf, select = -c(y, 
212 +
                                                                                                                                                                                                                 long, lat))), tag = "pred.y")  #covariate: remove resp and lat long to keep just covariates
213 +
      
214 +
      # model at response location (not prediction)
215 +
      A <- INLA::inla.spde.make.A(mesh = mesh, loc = cbind(train$long, train$lat))
216 +
      covariate_all = subset(train, select = -c(y, long, lat))  #remove resp and lat long to keep just covariates
217 +
      stk.y <- INLA::inla.stack(data = list(y = train$y), A = list(A, 1, 1), effects = list(list(spatial.field = 1:spde$n.spde), 
218 +
                                                                                            list(iid = 1:length(train$y)), list(y.intercept = rep(1, length(train$y)), covariate = covariate_all)), tag = "est.y")
219 +
      # put all together prediction + at response level
220 +
      stk.full <- INLA::inla.stack(stk.y, stk.pred)
221 +
      datastk  <-  INLA::inla.stack.data(stk.full, spde = spde)
222 +
      
223 +
      if (family == "bernoulli") {
224 +
        out <- INLA::inla(modform, family = "binomial", Ntrials = 1, data = datastk, control.predictor = list(A = INLA::inla.stack.A(stk.full), 
225 +
                                                                                                              link = 1), control.inla = list(int.strategy = int.strategy), ...)
226 +
        indpred <- INLA::inla.stack.index(stk.full, "pred.y")$data
227 +
        pall <- round(out$summary.fitted.values[indpred, ], 3)  #get mean, sd, etc of predicted reponse at test.df location
228 +
        p <- c(p, as.numeric(pall[1]))  #first element is the mean of predictions
229 +
        
230 +
        
231 +
      } else if (family == "binomial") {
232 +
        ntrials_sub <- c(ntrials[sqrt((dataframe$long - test.df[i, ]$long)^2 + (dataframe$lat - test.df[i, ]$lat)^2) > rad], test.ntrials[i])
233 +
        out <- INLA::inla(modform, family = "binomial", Ntrials = ntrials_sub, data = datastk, control.predictor = list(A = INLA::inla.stack.A(stk.full), 
234 +
                                                                                                                        link = 1), control.inla = list(int.strategy = int.strategy), ...)
235 +
        indpred <- INLA::inla.stack.index(stk.full, "pred.y")$data
236 +
        pall <- round(out$summary.fitted.values[indpred, ], 3)  #get mean, sd, etc of predicted reponse at test.df location
237 +
        p <- c(p, as.numeric(pall[1]))  #first element is the mean of predictions
238 +
        
239 +
        
240 +
      } else {
241 +
        out <- INLA::inla(modform, family = family, data = datastk, control.predictor = list(A = INLA::inla.stack.A(stk.full), 
242 +
                                                                                             link = 1), control.inla = list(int.strategy = int.strategy), ...)
243 +
        indpred <- INLA::inla.stack.index(stk.full, "pred.y")$data
244 +
        pall <- round(out$summary.fitted.values[indpred, ], 3)  #get mean, sd, etc of predicted reponse at test.df location
245 +
        p <- c(p, as.numeric(pall[1]))  #first element is the mean of predictions
246 +
      }
247 +
    }  #end of i loop
167 248
    
168 -
    # Stupid way of getting correct column names without collisions (almost certainly).
169 -
    colnames(dataframe)[colnames(dataframe) == y] <- "tmp_y_column_name_eeeeee"
170 -
    colnames(dataframe)[colnames(dataframe) == long] <- "tmp_long_column_name_eeeeee"
171 -
    colnames(dataframe)[colnames(dataframe) == lat] <- "tmp_lat_column_name_eeeeee"
249 +
    if (sqroot == TRUE & any(test.df$y < 0) == FALSE & any(p < 0) == FALSE) {
250 +
      # compute the sqrt of predicted and observed for gamma function
251 +
      p.res <- sqrt(test.df$y) - sqrt(p)
252 +
    } else {
253 +
      p.res <- test.df$y - p
254 +
    }
172 255
    
173 -
    colnames(dataframe)[colnames(dataframe) == "tmp_y_column_name_eeeeee"] <- "y"
174 -
    colnames(dataframe)[colnames(dataframe) == "tmp_long_column_name_eeeeee"] <- "long"
175 -
    colnames(dataframe)[colnames(dataframe) == "tmp_lat_column_name_eeeeee"] <- "lat"
256 +
    p.rmse <- sqrt(mean(p.res^2, na.rm = TRUE))
176 257
    
177 -
    measurevar <- "y"  #replace response by 'y' in model formula
178 -
    for (i in 1:length(modform)) {
179 -
        modform[i] <- gsub(".*~", "", modform[i])
180 -
        modform[i] <- paste(measurevar, paste(modform[i]), sep = " ~ ")
258 +
    if (mae == TRUE) {
259 +
      p.mae <- mean(abs(p.res), na.rm = TRUE)
260 +
    } else {
261 +
      p.mae <- NA
181 262
    }
263 +
    if (ds == TRUE) {
264 +
      p.ds <- mean((test.df$y - (mean(p, na.rm = TRUE) * stats::sd(p, na.rm = TRUE)))^2 + 2 * log(stats::sd(p, na.rm = TRUE)), na.rm = TRUE)
265 +
    } else {
266 +
      p.ds <- NA
267 +
    }  #good for general predictive model choice; stable for size imbalance in categorical predictors)
182 268
    
183 -
    # start of SLOO process
184 -
    test_rows <- sample(nrow(dataframe), ss)
185 -
    test.df <- dataframe[test_rows, ]  #sample subset of size ss from the dataframe
186 -
    test.ntrials <- ntrials[test_rows]
187 -
    slooresults <- list()
188 -
    slooresultsprint<- list()
189 -
    familys <- family
190 -
    modforms <- modform
191 -
    result.list <- expand.grid(familys, modforms)  #create all possible combinations of family and modforms
192 -
    result.list <- lapply(apply(result.list, 1, identity), unlist)  #put into a list
193 -
    for (j in 1:(length(result.list))) {
194 -
        family <- as.character(result.list[[j]][1])  #select family for j element
195 -
        modform <- stats::as.formula(as.character(result.list[[j]][2]))  #select modform for j element
196 -
        spde <- INLA::inla.spde2.matern(mesh, alpha = 2)  #SPDE computation
197 -
        
198 -
        for (i in 1:nrow(test.df)) {
199 -
            if (i == 1) {
200 -
                p <- c()
201 -
            }
202 -
            # Training data (test point & buffer removed)
203 -
            train <- dataframe[sqrt((dataframe$long - test.df[i, ]$long)^2 + (dataframe$lat - test.df[i, ]$lat)^2) > rad, ]
204 -
            # predictions in test locations
205 -
            predptdf <- test.df[i, ]  #predict on test point (one point only)
206 -
            A.pred <- INLA::inla.spde.make.A(mesh = mesh, loc = cbind(test.df[i, ]$long, test.df[i, ]$lat))
207 -
            stk.pred <- INLA::inla.stack(data = list(y = NA), A = list(A.pred, 1, 1, 1), effects = list(spatial.field = 1:spde$n.spde, 
208 -
                iid = 1:length(predptdf$y), y.intercept = rep(1, length(predptdf$y)), covariate = subset(predptdf, select = -c(y, 
209 -
                  long, lat))), tag = "pred.y")  #covariate: remove resp and lat long to keep just covariates
210 -
            
211 -
            # model at response location (not prediction)
212 -
            A <- INLA::inla.spde.make.A(mesh = mesh, loc = cbind(train$long, train$lat))
213 -
            covariate_all = subset(train, select = -c(y, long, lat))  #remove resp and lat long to keep just covariates
214 -
            stk.y <- INLA::inla.stack(data = list(y = train$y), A = list(A, 1, 1), effects = list(list(spatial.field = 1:spde$n.spde), 
215 -
                list(iid = 1:length(train$y)), list(y.intercept = rep(1, length(train$y)), covariate = covariate_all)), tag = "est.y")
216 -
            # put all together prediction + at response level
217 -
            stk.full <- INLA::inla.stack(stk.y, stk.pred)
218 -
            datastk  <-  INLA::inla.stack.data(stk.full, spde = spde)
219 -
            
220 -
            if (family == "bernoulli") {
221 -
                out <- INLA::inla(modform, family = "binomial", Ntrials = 1, data = datastk, control.predictor = list(A = INLA::inla.stack.A(stk.full), 
222 -
                  link = 1), control.inla = list(int.strategy = int.strategy), ...)
223 -
                indpred <- INLA::inla.stack.index(stk.full, "pred.y")$data
224 -
                pall <- round(out$summary.fitted.values[indpred, ], 3)  #get mean, sd, etc of predicted reponse at test.df location
225 -
                p <- c(p, as.numeric(pall[1]))  #first element is the mean of predictions
226 -
                
227 -
                
228 -
            } else if (family == "binomial") {
229 -
                ntrials_sub <- c(ntrials[sqrt((dataframe$long - test.df[i, ]$long)^2 + (dataframe$lat - test.df[i, ]$lat)^2) > rad], test.ntrials[i])
230 -
                out <- INLA::inla(modform, family = "binomial", Ntrials = ntrials_sub, data = datastk, control.predictor = list(A = INLA::inla.stack.A(stk.full), 
231 -
                  link = 1), control.inla = list(int.strategy = int.strategy), ...)
232 -
                indpred <- INLA::inla.stack.index(stk.full, "pred.y")$data
233 -
                pall <- round(out$summary.fitted.values[indpred, ], 3)  #get mean, sd, etc of predicted reponse at test.df location
234 -
                p <- c(p, as.numeric(pall[1]))  #first element is the mean of predictions
235 -
                
236 -
                
237 -
            } else {
238 -
                out <- INLA::inla(modform, family = family, data = datastk, control.predictor = list(A = INLA::inla.stack.A(stk.full), 
239 -
                  link = 1), control.inla = list(int.strategy = int.strategy), ...)
240 -
                indpred <- INLA::inla.stack.index(stk.full, "pred.y")$data
241 -
                pall <- round(out$summary.fitted.values[indpred, ], 3)  #get mean, sd, etc of predicted reponse at test.df location
242 -
                p <- c(p, as.numeric(pall[1]))  #first element is the mean of predictions
243 -
            }
244 -
        }  #end of i loop
245 -
        
246 -
        if (sqroot == TRUE & any(test.df$y < 0) == FALSE & any(p < 0) == FALSE) {
247 -
            # compute the sqrt of predicted and observed for gamma function
248 -
            p.res <- sqrt(test.df$y) - sqrt(p)
249 -
        } else {
250 -
            p.res <- test.df$y - p
251 -
        }
252 -
        
253 -
        p.rmse <- sqrt(mean(p.res^2, na.rm = TRUE))
254 -
        
255 -
        if (mae == TRUE) {
256 -
            p.mae <- mean(abs(p.res), na.rm = TRUE)
257 -
        } else {
258 -
            p.mae <- NA
259 -
        }
260 -
        if (ds == TRUE) {
261 -
            p.ds <- mean((test.df$y - (mean(p, na.rm = TRUE) * stats::sd(p, na.rm = TRUE)))^2 + 2 * log(stats::sd(p, na.rm = TRUE)), na.rm = TRUE)
262 -
        } else {
263 -
            p.ds <- NA
264 -
        }  #good for general predictive model choice; stable for size imbalance in categorical predictors)
265 -
        
266 -
        slooresults[[j]] <- list(Observed_response = test.df$y, Predictions = p, Residuals = p.res, RMSE = p.rmse, 
267 -
            MAE = p.mae, DS = p.ds, family = family, ntrials = test.ntrials, test.df = test.df, Rownames_test = as.numeric(rownames(test.df)))
268 -
        
269 -
        slooresultsprint <- slooresults[[j]]
270 -
        slooresultsprint[[1]]<-round(as.numeric(slooresultsprint[[1]]),3)
271 -
        slooresultsprint[[2]]<-round(as.numeric(slooresultsprint[[2]]),3)
272 -
        slooresultsprint[[3]]<-round(as.numeric(slooresultsprint[[3]]),3)
273 -
        slooresultsprint[[4]]<-round(as.numeric(slooresultsprint[[4]]),3)
274 -
        slooresultsprint[[5]]<-round(as.numeric(slooresultsprint[[5]]),3)
275 -
        slooresultsprint[[6]]<-round(as.numeric(slooresultsprint[[6]]),3)
276 -
        if (print == TRUE) {
277 -
          message("\n")
278 -
          message("Summary of the Spatial leave-one-out analysis")
279 -
          message("#############################################", "\n")
280 -
          message("MODEL", "", j, "\n")
281 -
          print(slooresultsprint[-c(9:10)])
282 -
          message("End summary of the Spatial leave-one-out analysis")
283 -
          message("#################################################", "\n")
284 -
        }
285 -
    }  #end of j loop
269 +
    slooresults[[j]] <- list(Observed_response = test.df$y, Predictions = p, Residuals = p.res, RMSE = p.rmse, 
270 +
                             MAE = p.mae, DS = p.ds, family = family, ntrials = test.ntrials, test.df = test.df, Rownames_test = as.numeric(rownames(test.df)))
286 271
    
287 -
    # plot locations of observation and test points
288 -
    par(mfrow = grDevices::n2mfrow(length(result.list)), family = "mono", oma = c(0, 0, 2, 0))
289 -
    for (j in 1:(length(result.list))) {
290 -
        slooplot(alpha = alpha, df = slooresults[[j]], mae = mae, ds = ds, family = as.character(result.list[[j]][1]), ntrials = test.ntrials, 
291 -
            sqroot = sqroot)
292 -
        graphics::mtext(paste0("MODEL", "", j), side = 3, adj = 1, cex = 0.8, line = 1.6)
272 +
    slooresultsprint <- slooresults[[j]]
273 +
    slooresultsprint[[1]]<-round(as.numeric(slooresultsprint[[1]]),3)
274 +
    slooresultsprint[[2]]<-round(as.numeric(slooresultsprint[[2]]),3)
275 +
    slooresultsprint[[3]]<-round(as.numeric(slooresultsprint[[3]]),3)
276 +
    slooresultsprint[[4]]<-round(as.numeric(slooresultsprint[[4]]),3)
277 +
    slooresultsprint[[5]]<-round(as.numeric(slooresultsprint[[5]]),3)
278 +
    slooresultsprint[[6]]<-round(as.numeric(slooresultsprint[[6]]),3)
279 +
    if (print == TRUE) {
280 +
      message("\n")
281 +
      message("Summary of the Spatial leave-one-out analysis")
282 +
      message("#############################################", "\n")
283 +
      message("MODEL", "", j, "\n")
284 +
      print(slooresultsprint[-c(9:10)])
285 +
      message("End summary of the Spatial leave-one-out analysis")
286 +
      message("#################################################", "\n")
293 287
    }
294 -
    graphics::title(paste("SLOO with number of iterations =", "", ss), outer = TRUE)
295 -
    # end plot
296 -
    # plot locations of observation and test points
297 -
    par(mfrow = c(1, 1), family = "mono")
298 -
    sloopoint(points = cbind(dataframe$long, dataframe$lat), test = cbind(slooresults[[1]]$test.df$long, slooresults[[1]]$test.df$lat), 
299 -
        rad = rad)
300 -
    return(slooresults)#save the results as an object
301 -
}
288 +
  }  #end of j loop
289 +
  if (plot == TRUE) {
290 +
    
291 +
  # plot locations of observation and test points
292 +
  par(mfrow = grDevices::n2mfrow(length(result.list)), family = "mono", mar=c(5,4,4,2) + 0.1,oma = c(0, 0, 2, 0))
293 +
  for (j in 1:(length(result.list))) {
294 +
    slooplot(alpha = alpha, df = slooresults[[j]], mae = mae, ds = ds, family = as.character(result.list[[j]][1]), ntrials = test.ntrials, 
295 +
             sqroot = sqroot)
296 +
    graphics::mtext(paste0("MODEL", "", j), side = 3, adj = 1, cex = 0.8, line = 1.6)
297 +
  }
298 +
  graphics::title(paste("SLOO with number of iterations =", "", ss), outer = TRUE)
299 +
  # end plot
300 +
  # plot locations of observation and test points
301 +
  par(mfrow = c(1, 1), family = "mono", mar=c(5,4,4,2) + 0.1,oma = c(2, 2, 2, 2))
302 +
  sloopoint(points = cbind(dataframe$long, dataframe$lat), test = cbind(slooresults[[1]]$test.df$long, slooresults[[1]]$test.df$lat), 
303 +
            rad = rad)
304 +
  }
305 +
  return(slooresults)#save the results as an object
306 +
  }
307 +

@@ -0,0 +1,129 @@
Loading
1 +
#' A function to plot outputs of a spatial leave-one-out cross-validation in R using INLA.
2 +
#' 
3 +
#' This function plots the output of a spatial leave-one-out cross-validation (SLOO-CV) of one or several models running on INLA.
4 +
#' @param df dataframe (output of inla.sloo), which include the "Observed_response", "Predictions", "Residuals", RMSE, MAE (optional), and DS(optional)
5 +
#' @param alpha Numeric value (0,...,1) giving the threshold for computing confidence intervals (1-alpha) of rmse and mae estimation
6 +
#' @param mae If TRUE, compute the mean absolute error (mae) and the root mean square error (rmse). If FALSE, compute the rmse only.
7 +
#' @param ds If TRUE, compute the Dawid-Sebastiani score (ds). If FALSE, does not compute ds.
8 +
#' @param family Character string or list of Character string giving the family or families, respectively of the response in INLA format 
9 +
#' @param ntrials Numeric value (1,2,...) setting the number of trials for a Binomial family
10 +
#' @param sqroot If TRUE, compute the square root of the observed and predicted values to generate the rmse and/or mae. 
11 +
#' 
12 +
#' @export
13 +
#' @name slooplot.fun
14 +
#' @examples 
15 +
#' \dontrun{
16 +
#'  # SLOO function with one model formula (Bernoulli)
17 +
#'  df<-data.frame(Residuals=runif(10, 0.0, 1.0),RMSE=runif(10, 0.0, 2.0),MAE=runif(10, 0.0, 2.0),
18 +
#'  Observed_response= sample(c(0,1), replace=TRUE, size=10),Predictions=runif(10, 0.0, 1.0))
19 +
#'  alpha = 0.05
20 +
#'  family = bernoulli
21 +
#'  slooplot1<-slooplot.fun(df=df, alpha=0.05,mae=TRUE,ds=FALSE,family='bernoulli',sqroot=FALSE)
22 +
#'  }
23 +
24 +
slooplot.fun <- function(alpha = alpha, df = df, mae = mae, ds = ds, family = family, ntrials = ntrials, sqroot = sqroot) {
25 +
  # compute confidence intervals for rmse and mae confidence intervals
26 +
  rmse.LCI <- sqrt(length(df$Residuals)/stats::qchisq(1 - alpha/2, df = length(df$Residuals))) * df$RMSE
27 +
  rmse.HCI <- sqrt(length(df$Residuals)/stats::qchisq(alpha/2, df = length(df$Residuals))) * df$RMSE
28 +
  if (mae == TRUE) {
29 +
    mae.LCI <- sqrt(length(df$Residuals)/stats::qchisq(1 - alpha/2, df = length(df$Residuals))) * df$MAE
30 +
    mae.HCI <- sqrt(length(df$Residuals)/stats::qchisq(alpha/2, df = length(df$Residuals))) * df$MAE
31 +
  }
32 +
  if (family == "binomial") {
33 +
    df$Predictions <- df$Predictions*ntrials
34 +
  }
35 +
  
36 +
  if (sqroot == TRUE) {
37 +
    plotx <- sqrt(df$Observed_response)
38 +
    ploty <- sqrt(df$Predictions)
39 +
    xlabtext <- paste0("Sqrt.observed response", " ", "(", family, ")")
40 +
    ylabtext <- "Sqrt.prediction"
41 +
  } else {
42 +
    plotx <- df$Observed_response
43 +
    ploty <- df$Predictions
44 +
    xlabtext <- paste0("Observed response", " ", "(", family, ")")
45 +
    ylabtext <- "Prediction"
46 +
  }
47 +
  if (family == "gamma") {
48 +
    slooplot <- plot(plotx, ploty, pch = "+", xlab = xlabtext, ylab = ylabtext, pty = "s", asp = 1)
49 +
    abline(a = 0, b = 1, col = "grey")
50 +
    graphics::mtext(paste0("rmse mean sqrt", " ", "(", 1 - alpha, "% CI):", " ", round(df$RMSE, 2), " ", "(", round(rmse.LCI, 2), 
51 +
                           ";", round(rmse.HCI, 2), ")"), side = 3, adj = 0, cex = 0.8)
52 +
    if (mae == TRUE) {
53 +
      graphics::mtext(paste0(" ", "mae mean sqrt", " ", "(", 1 - alpha, "% CI):", " ", round(df$MAE, 2), " ", "(", round(mae.LCI, 
54 +
                                                                                                                         2), ";", round(mae.HCI, 2), ")"), side = 3, adj = 0, cex = 0.8, line = 0.8)
55 +
    }
56 +
    if (ds == TRUE) {
57 +
      graphics::mtext(paste0("dawid-sebastiani score:", " ", " ", " ", " ", " ", round(df$DS, 2)), side = 3, adj = 0, cex = 0.8, 
58 +
                      line = 1.6)
59 +
    }
60 +
    cat(slooplot)
61 +
    
62 +
  } else if (family == "poisson") {
63 +
    slooplot <- plot(plotx, ploty, pch = "+", xlab = xlabtext, ylab = ylabtext, pty = "s", asp = 1)
64 +
    abline(a = 0, b = 1, col = "grey")
65 +
    graphics::mtext(paste0("rmse mean", " ", "(", 1 - alpha, "% CI):", " ", round(df$RMSE, 2), " ", "(", round(rmse.LCI, 2), ";", 
66 +
                           round(rmse.HCI, 2), ")"), side = 3, adj = 0, cex = 0.8)
67 +
    if (mae == TRUE) {
68 +
      graphics::mtext(paste0(" ", "mae mean", " ", "(", 1 - alpha, "% CI):", " ", round(df$MAE, 2), " ", "(", round(mae.LCI, 2), 
69 +
                             ";", round(mae.HCI, 2), ")"), side = 3, adj = 0, cex = 0.8, line = 0.8)
70 +
    }
71 +
    if (ds == TRUE) {
72 +
      graphics::mtext(paste0("dawid-sebastiani score:", " ", " ", " ", " ", " ", round(df$DS, 2)), side = 3, adj = 0, cex = 0.8, 
73 +
                      line = 1.6)
74 +
    }
75 +
    cat(slooplot)
76 +
    
77 +
  } else if (family == "bernoulli") {
78 +
    
79 +
    yjitter <- jitter(plotx, factor = 0.01)
80 +
    slooplot <- plot(yjitter, ploty, pch = "+", xlab = xlabtext, ylab = ylabtext, xaxt = "n", pty = "s", asp = 1)
81 +
    # abline(a = 0, b = 1, col = 'grey')
82 +
    graphics::axis(1, at = c(0, 1), labels = c(0, 1))
83 +
    graphics::mtext(paste0("rmse mean", " ", "(", 1 - alpha, "% CI):", " ", round(df$RMSE, 2), " ", "(", round(rmse.LCI, 2), ";", 
84 +
                           round(rmse.HCI, 2), ")"), side = 3, adj = 0, cex = 0.8)
85 +
    if (mae == TRUE) {
86 +
      graphics::mtext(paste0(" ", "mae mean", " ", "(", 1 - alpha, "% CI):", " ", round(df$MAE, 2), " ", "(", round(mae.LCI, 2), 
87 +
                             ";", round(mae.HCI, 2), ")"), side = 3, adj = 0, cex = 0.8, line = 0.8)
88 +
    }
89 +
    if (ds == TRUE) {
90 +
      graphics::mtext(paste0("dawid-sebastiani score:", " ", " ", " ", " ", " ", round(df$DS, 2)), side = 3, adj = 0, cex = 0.8, 
91 +
                      line = 1.6)
92 +
    }
93 +
    cat(slooplot)
94 +
    
95 +
  } else if (family == "binomial") {
96 +
    
97 +
    slooplot <- plot(plotx, ploty, pch = "+", xlab = xlabtext, ylab = ylabtext, pty = "s", asp = 1)
98 +
    abline(a = 0, b = 1, col = "grey")
99 +
    graphics::axis(1, at = c(0, 1), labels = c(0, 1))
100 +
    graphics::mtext(paste0("rmse mean", " ", "(", 1 - alpha, "% CI):", " ", round(df$RMSE, 2), " ", "(", round(rmse.LCI, 2), ";", 
101 +
                           round(rmse.HCI, 2), ")"), side = 3, adj = 0, cex = 0.8)
102 +
    if (mae == TRUE) {
103 +
      graphics::mtext(paste0(" ", "mae mean", " ", "(", 1 - alpha, "% CI):", " ", round(df$MAE, 2), " ", "(", round(mae.LCI, 2), 
104 +
                             ";", round(mae.HCI, 2), ")"), side = 3, adj = 0, cex = 0.8, line = 0.8)
105 +
    }
106 +
    if (ds == TRUE) {
107 +
      graphics::mtext(paste0("dawid-sebastiani score:", " ", " ", " ", " ", " ", round(df$DS, 2)), side = 3, adj = 0, cex = 0.8, 
108 +
                      line = 1.6)
109 +
    }
110 +
    cat(slooplot)
111 +
    
112 +
  } else {
113 +
    slooplot <- plot(plotx, ploty, pch = "+", xlab = xlabtext, ylab = ylabtext, pty = "s", asp = 1)
114 +
    abline(a = 0, b = 1, col = "grey")
115 +
    graphics::mtext(paste0("rmse mean", " ", "(", 1 - alpha, "% CI):", " ", round(df$RMSE, 2), " ", "(", round(rmse.LCI, 2), ";", 
116 +
                           round(rmse.HCI, 2), ")"), side = 3, adj = 0, cex = 0.8)
117 +
    if (mae == TRUE) {
118 +
      graphics::mtext(paste0(" ", "mae mean", " ", "(", 1 - alpha, "% CI):", " ", round(df$MAE, 2), " ", "(", round(mae.LCI, 2), 
119 +
                             ";", round(mae.HCI, 2), ")"), side = 3, adj = 0, cex = 0.8, line = 0.8)
120 +
    }
121 +
    if (ds == TRUE) {
122 +
      graphics::mtext(paste0("dawid-sebastiani score:", " ", " ", " ", " ", " ", round(df$DS, 2)), side = 3, adj = 0, cex = 0.8, 
123 +
                      line = 1.6)
124 +
    }
125 +
    cat(slooplot)
126 +
  }
127 +
}
128 +
129 +
  

@@ -1,16 +1,16 @@
Loading
1 1
#' A function to plot outputs of a spatial leave-one-out cross-validation in R using INLA.
2 2
#' 
3 3
#' This function plots the output of a spatial leave-one-out cross-validation (SLOO-CV) of one or several models running on INLA.
4 -
#' @param df = dataframe (output of inla.sloo), which include the "Observed_response", "Predictions", "Residuals", RMSE, MAE (optional), and DS(optional)
5 -
#' @param alpha = Numeric value (0,...,1) giving the threshold for computing confidence intervals (1-alpha) of rmse and mae estimation
6 -
#' @param mae = If TRUE, compute the mean absolute error (mae) and the root mean square error (rmse). If FALSE, compute the rmse only.
7 -
#' @param ds = If TRUE, compute the Dawid-Sebastiani score (ds). If FALSE, does not compute ds.
8 -
#' @param family = Character string or list of Character string giving the family or families, respectively of the response in INLA format 
9 -
#' @param ntrials = Numeric value (1,2,...) setting the number of trials for a Binomial family
10 -
#' @param sqroot = If TRUE, compute the square root of the observed and predicted values to generate the rmse and/or mae. 
4 +
#' @param df dataframe (output of inla.sloo), which include the "Observed_response", "Predictions", "Residuals", RMSE, MAE (optional), and DS(optional)
5 +
#' @param alpha Numeric value (0,...,1) giving the threshold for computing confidence intervals (1-alpha) of rmse and mae estimation
6 +
#' @param mae If TRUE, compute the mean absolute error (mae) and the root mean square error (rmse). If FALSE, compute the rmse only.
7 +
#' @param ds If TRUE, compute the Dawid-Sebastiani score (ds). If FALSE, does not compute ds.
8 +
#' @param family Character string or list of Character string giving the family or families, respectively of the response in INLA format 
9 +
#' @param ntrials Numeric value (1,2,...) setting the number of trials for a Binomial family
10 +
#' @param sqroot If TRUE, compute the square root of the observed and predicted values to generate the rmse and/or mae. 
11 11
#' 
12 +
#' @export
12 13
#' @name slooplot
13 -
#' 
14 14
#' @examples 
15 15
#' \dontrun{
16 16
#'  # SLOO function with one model formula (Bernoulli)
@@ -21,105 +21,108 @@
Loading
21 21
#'  slooplot1<-slooplot(df=df, alpha=0.05,mae=TRUE,ds=FALSE,family='bernoulli',sqroot=FALSE)
22 22
#'  }
23 23
24 -
slooplot <- function(alpha, df, mae, ds, family, ntrials, sqroot) {
25 -
    # compute confidence intervals for rmse and mae confidence intervals
26 -
    rmse.LCI <- sqrt(length(df$Residuals)/stats::qchisq(1 - alpha/2, df = length(df$Residuals))) * df$RMSE
27 -
    rmse.HCI <- sqrt(length(df$Residuals)/stats::qchisq(alpha/2, df = length(df$Residuals))) * df$RMSE
24 +
slooplot <- function(alpha = alpha, df = df, mae = mae, ds = ds, family = family, ntrials = ntrials, sqroot = sqroot) {
25 +
  # compute confidence intervals for rmse and mae confidence intervals
26 +
  rmse.LCI <- sqrt(length(df$Residuals)/stats::qchisq(1 - alpha/2, df = length(df$Residuals))) * df$RMSE
27 +
  rmse.HCI <- sqrt(length(df$Residuals)/stats::qchisq(alpha/2, df = length(df$Residuals))) * df$RMSE
28 +
  if (mae == TRUE) {
29 +
    mae.LCI <- sqrt(length(df$Residuals)/stats::qchisq(1 - alpha/2, df = length(df$Residuals))) * df$MAE
30 +
    mae.HCI <- sqrt(length(df$Residuals)/stats::qchisq(alpha/2, df = length(df$Residuals))) * df$MAE
31 +
  }
32 +
  if (family == "binomial") {
33 +
    df$Predictions <- df$Predictions*ntrials
34 +
  }
35 +
  
36 +
  if (sqroot == TRUE) {
37 +
    plotx <- sqrt(df$Observed_response)
38 +
    ploty <- sqrt(df$Predictions)
39 +
    xlabtext <- paste0("Sqrt.observed response", " ", "(", family, ")")
40 +
    ylabtext <- "Sqrt.prediction"
41 +
  } else {
42 +
    plotx <- df$Observed_response
43 +
    ploty <- df$Predictions
44 +
    xlabtext <- paste0("Observed response", " ", "(", family, ")")
45 +
    ylabtext <- "Prediction"
46 +
  }
47 +
  if (family == "gamma") {
48 +
    slooplot <- plot(plotx, ploty, pch = "+", xlab = xlabtext, ylab = ylabtext, pty = "s", asp = 1)
49 +
    abline(a = 0, b = 1, col = "grey")
50 +
    graphics::mtext(paste0("rmse mean sqrt", " ", "(", 1 - alpha, "% CI):", " ", round(df$RMSE, 2), " ", "(", round(rmse.LCI, 2), 
51 +
                           ";", round(rmse.HCI, 2), ")"), side = 3, adj = 0, cex = 0.8)
52 +
    if (mae == TRUE) {
53 +
      graphics::mtext(paste0(" ", "mae mean sqrt", " ", "(", 1 - alpha, "% CI):", " ", round(df$MAE, 2), " ", "(", round(mae.LCI, 
54 +
                                                                                                                         2), ";", round(mae.HCI, 2), ")"), side = 3, adj = 0, cex = 0.8, line = 0.8)
55 +
    }
56 +
    if (ds == TRUE) {
57 +
      graphics::mtext(paste0("dawid-sebastiani score:", " ", " ", " ", " ", " ", round(df$DS, 2)), side = 3, adj = 0, cex = 0.8, 
58 +
                      line = 1.6)
59 +
    }
60 +
    cat(slooplot)
61 +
    
62 +
  } else if (family == "poisson") {
63 +
    slooplot <- plot(plotx, ploty, pch = "+", xlab = xlabtext, ylab = ylabtext, pty = "s", asp = 1)
64 +
    abline(a = 0, b = 1, col = "grey")
65 +
    graphics::mtext(paste0("rmse mean", " ", "(", 1 - alpha, "% CI):", " ", round(df$RMSE, 2), " ", "(", round(rmse.LCI, 2), ";", 
66 +
                           round(rmse.HCI, 2), ")"), side = 3, adj = 0, cex = 0.8)
67 +
    if (mae == TRUE) {
68 +
      graphics::mtext(paste0(" ", "mae mean", " ", "(", 1 - alpha, "% CI):", " ", round(df$MAE, 2), " ", "(", round(mae.LCI, 2), 
69 +
                             ";", round(mae.HCI, 2), ")"), side = 3, adj = 0, cex = 0.8, line = 0.8)
70 +
    }
71 +
    if (ds == TRUE) {
72 +
      graphics::mtext(paste0("dawid-sebastiani score:", " ", " ", " ", " ", " ", round(df$DS, 2)), side = 3, adj = 0, cex = 0.8, 
73 +
                      line = 1.6)
74 +
    }
75 +
    cat(slooplot)
76 +
    
77 +
  } else if (family == "bernoulli") {
78 +
    
79 +
    yjitter <- jitter(plotx, factor = 0.01)
80 +
    slooplot <- plot(yjitter, ploty, pch = "+", xlab = xlabtext, ylab = ylabtext, xaxt = "n", pty = "s", asp = 1)
81 +
    # abline(a = 0, b = 1, col = 'grey')
82 +
    graphics::axis(1, at = c(0, 1), labels = c(0, 1))
83 +
    graphics::mtext(paste0("rmse mean", " ", "(", 1 - alpha, "% CI):", " ", round(df$RMSE, 2), " ", "(", round(rmse.LCI, 2), ";", 
84 +
                           round(rmse.HCI, 2), ")"), side = 3, adj = 0, cex = 0.8)
85 +
    if (mae == TRUE) {
86 +
      graphics::mtext(paste0(" ", "mae mean", " ", "(", 1 - alpha, "% CI):", " ", round(df$MAE, 2), " ", "(", round(mae.LCI, 2), 
87 +
                             ";", round(mae.HCI, 2), ")"), side = 3, adj = 0, cex = 0.8, line = 0.8)
88 +
    }
89 +
    if (ds == TRUE) {
90 +
      graphics::mtext(paste0("dawid-sebastiani score:", " ", " ", " ", " ", " ", round(df$DS, 2)), side = 3, adj = 0, cex = 0.8, 
91 +
                      line = 1.6)
92 +
    }
93 +
    cat(slooplot)
94 +
    
95 +
  } else if (family == "binomial") {
96 +
    
97 +
    slooplot <- plot(plotx, ploty, pch = "+", xlab = xlabtext, ylab = ylabtext, pty = "s", asp = 1)
98 +
    abline(a = 0, b = 1, col = "grey")
99 +
    graphics::axis(1, at = c(0, 1), labels = c(0, 1))
100 +
    graphics::mtext(paste0("rmse mean", " ", "(", 1 - alpha, "% CI):", " ", round(df$RMSE, 2), " ", "(", round(rmse.LCI, 2), ";", 
101 +
                           round(rmse.HCI, 2), ")"), side = 3, adj = 0, cex = 0.8)
28 102
    if (mae == TRUE) {
29 -
        mae.LCI <- sqrt(length(df$Residuals)/stats::qchisq(1 - alpha/2, df = length(df$Residuals))) * df$MAE
30 -
        mae.HCI <- sqrt(length(df$Residuals)/stats::qchisq(alpha/2, df = length(df$Residuals))) * df$MAE
103 +
      graphics::mtext(paste0(" ", "mae mean", " ", "(", 1 - alpha, "% CI):", " ", round(df$MAE, 2), " ", "(", round(mae.LCI, 2), 
104 +
                             ";", round(mae.HCI, 2), ")"), side = 3, adj = 0, cex = 0.8, line = 0.8)
105 +
    }
106 +
    if (ds == TRUE) {
107 +
      graphics::mtext(paste0("dawid-sebastiani score:", " ", " ", " ", " ", " ", round(df$DS, 2)), side = 3, adj = 0, cex = 0.8, 
108 +
                      line = 1.6)
31 109
    }
32 -
    if (sqroot == TRUE) {
33 -
        plotx <- sqrt(df$Observed_response)
34 -
        ploty <- sqrt(df$Predictions)
35 -
        xlabtext <- paste0("Sqrt.observed response", " ", "(", family, ")")
36 -
        ylabtext <- "Sqrt.prediction"
37 -
    } else {
38 -
        plotx <- df$Observed_response
39 -
        ploty <- df$Predictions
40 -
        xlabtext <- paste0("Observed response", " ", "(", family, ")")
41 -
        ylabtext <- "Prediction"
110 +
    cat(slooplot)
111 +
    
112 +
  } else {
113 +
    slooplot <- plot(plotx, ploty, pch = "+", xlab = xlabtext, ylab = ylabtext, pty = "s", asp = 1)
114 +
    abline(a = 0, b = 1, col = "grey")
115 +
    graphics::mtext(paste0("rmse mean", " ", "(", 1 - alpha, "% CI):", " ", round(df$RMSE, 2), " ", "(", round(rmse.LCI, 2), ";", 
116 +
                           round(rmse.HCI, 2), ")"), side = 3, adj = 0, cex = 0.8)
117 +
    if (mae == TRUE) {
118 +
      graphics::mtext(paste0(" ", "mae mean", " ", "(", 1 - alpha, "% CI):", " ", round(df$MAE, 2), " ", "(", round(mae.LCI, 2), 
119 +
                             ";", round(mae.HCI, 2), ")"), side = 3, adj = 0, cex = 0.8, line = 0.8)
42 120
    }
43 -
    if (family == "gamma") {
44 -
        slooplot <- plot(plotx, ploty, pch = "+", xlab = xlabtext, ylab = ylabtext, pty = "s", asp = 1)
45 -
        abline(a = 0, b = 1, col = "grey")
46 -
        graphics::mtext(paste0("rmse mean sqrt", " ", "(", 1 - alpha, "% CI):", " ", round(df$RMSE, 2), " ", "(", round(rmse.LCI, 2), 
47 -
            ";", round(rmse.HCI, 2), ")"), side = 3, adj = 0, cex = 0.8)
48 -
        if (mae == TRUE) {
49 -
            graphics::mtext(paste0(" ", "mae mean sqrt", " ", "(", 1 - alpha, "% CI):", " ", round(df$MAE, 2), " ", "(", round(mae.LCI, 
50 -
                2), ";", round(mae.HCI, 2), ")"), side = 3, adj = 0, cex = 0.8, line = 0.8)
51 -
        }
52 -
        if (ds == TRUE) {
53 -
            graphics::mtext(paste0("dawid-sebastiani score:", " ", " ", " ", " ", " ", round(df$DS, 2)), side = 3, adj = 0, cex = 0.8, 
54 -
                line = 1.6)
55 -
        }
56 -
        cat(slooplot)
57 -
        
58 -
    } else if (family == "poisson") {
59 -
        slooplot <- plot(plotx, ploty, pch = "+", xlab = xlabtext, ylab = ylabtext, pty = "s", asp = 1)
60 -
        abline(a = 0, b = 1, col = "grey")
61 -
        graphics::mtext(paste0("rmse mean", " ", "(", 1 - alpha, "% CI):", " ", round(df$RMSE, 2), " ", "(", round(rmse.LCI, 2), ";", 
62 -
            round(rmse.HCI, 2), ")"), side = 3, adj = 0, cex = 0.8)
63 -
        if (mae == TRUE) {
64 -
            graphics::mtext(paste0(" ", "mae mean", " ", "(", 1 - alpha, "% CI):", " ", round(df$MAE, 2), " ", "(", round(mae.LCI, 2), 
65 -
                ";", round(mae.HCI, 2), ")"), side = 3, adj = 0, cex = 0.8, line = 0.8)
66 -
        }
67 -
        if (ds == TRUE) {
68 -
            graphics::mtext(paste0("dawid-sebastiani score:", " ", " ", " ", " ", " ", round(df$DS, 2)), side = 3, adj = 0, cex = 0.8, 
69 -
                line = 1.6)
70 -
        }
71 -
        cat(slooplot)
72 -
        
73 -
    } else if (family == "bernoulli") {
74 -
        
75 -
        yjitter <- jitter(plotx, factor = 0.01)
76 -
        slooplot <- plot(yjitter, ploty, pch = "+", xlab = xlabtext, ylab = ylabtext, xaxt = "n", pty = "s", asp = 1)
77 -
        # abline(a = 0, b = 1, col = 'grey')
78 -
        graphics::axis(1, at = c(0, 1), labels = c(0, 1))
79 -
        graphics::mtext(paste0("rmse mean", " ", "(", 1 - alpha, "% CI):", " ", round(df$RMSE, 2), " ", "(", round(rmse.LCI, 2), ";", 
80 -
            round(rmse.HCI, 2), ")"), side = 3, adj = 0, cex = 0.8)
81 -
        if (mae == TRUE) {
82 -
            graphics::mtext(paste0(" ", "mae mean", " ", "(", 1 - alpha, "% CI):", " ", round(df$MAE, 2), " ", "(", round(mae.LCI, 2), 
83 -
                ";", round(mae.HCI, 2), ")"), side = 3, adj = 0, cex = 0.8, line = 0.8)
84 -
        }
85 -
        if (ds == TRUE) {
86 -
            graphics::mtext(paste0("dawid-sebastiani score:", " ", " ", " ", " ", " ", round(df$DS, 2)), side = 3, adj = 0, cex = 0.8, 
87 -
                line = 1.6)
88 -
        }
89 -
        cat(slooplot)
90 -
        
91 -
    } else if (family == "binomial") {
92 -
        
93 -
        yjitter <- jitter(plotx, factor = 0.01)
94 -
        slooplot <- plot(yjitter, ploty, pch = "+", xlab = xlabtext, ylab = ylabtext, xaxt = "n", pty = "s", asp = 1)
95 -
        abline(a = 0, b = 1, col = "grey")
96 -
        graphics::axis(1, at = c(0, 1), labels = c(0, 1))
97 -
        graphics::mtext(paste0("rmse mean", " ", "(", 1 - alpha, "% CI):", " ", round(df$RMSE, 2), " ", "(", round(rmse.LCI, 2), ";", 
98 -
            round(rmse.HCI, 2), ")"), side = 3, adj = 0, cex = 0.8)
99 -
        if (mae == TRUE) {
100 -
            graphics::mtext(paste0(" ", "mae mean", " ", "(", 1 - alpha, "% CI):", " ", round(df$MAE, 2), " ", "(", round(mae.LCI, 2), 
101 -
                ";", round(mae.HCI, 2), ")"), side = 3, adj = 0, cex = 0.8, line = 0.8)
102 -
        }
103 -
        if (ds == TRUE) {
104 -
            graphics::mtext(paste0("dawid-sebastiani score:", " ", " ", " ", " ", " ", round(df$DS, 2)), side = 3, adj = 0, cex = 0.8, 
105 -
                line = 1.6)
106 -
        }
107 -
        cat(slooplot)
108 -
        
109 -
    } else {
110 -
        slooplot <- plot(plotx, ploty, pch = "+", xlab = xlabtext, ylab = ylabtext, pty = "s", asp = 1)
111 -
        abline(a = 0, b = 1, col = "grey")
112 -
        graphics::mtext(paste0("rmse mean", " ", "(", 1 - alpha, "% CI):", " ", round(df$RMSE, 2), " ", "(", round(rmse.LCI, 2), ";", 
113 -
            round(rmse.HCI, 2), ")"), side = 3, adj = 0, cex = 0.8)
114 -
        if (mae == TRUE) {
115 -
            graphics::mtext(paste0(" ", "mae mean", " ", "(", 1 - alpha, "% CI):", " ", round(df$MAE, 2), " ", "(", round(mae.LCI, 2), 
116 -
                ";", round(mae.HCI, 2), ")"), side = 3, adj = 0, cex = 0.8, line = 0.8)
117 -
        }
118 -
        if (ds == TRUE) {
119 -
            graphics::mtext(paste0("dawid-sebastiani score:", " ", " ", " ", " ", " ", round(df$DS, 2)), side = 3, adj = 0, cex = 0.8, 
120 -
                line = 1.6)
121 -
        }
122 -
        cat(slooplot)
121 +
    if (ds == TRUE) {
122 +
      graphics::mtext(paste0("dawid-sebastiani score:", " ", " ", " ", " ", " ", round(df$DS, 2)), side = 3, adj = 0, cex = 0.8, 
123 +
                      line = 1.6)
123 124
    }
125 +
    cat(slooplot)
126 +
  }
124 127
}
125 128
Files Coverage
R 90.49%
Project Totals (14 files) 90.49%
Notifications are pending CI completion. Waiting for GitHub's status webhook to queue notifications. Push notifications now.
Sunburst
The inner-most circle is the entire project, moving away from the center are folders then, finally, a single file. The size and color of each slice is representing the number of statements and the coverage, respectively.
Icicle
The top section represents the entire project. Proceeding with folders and finally individual files. The size and color of each slice is representing the number of statements and the coverage, respectively.
Grid
Each block represents a single file in the project. The size and color of each block is represented by the number of statements and the coverage, respectively.
Loading