From 81621ff50c4c404fccd0b2cee765bf5cd7f87ad1 Mon Sep 17 00:00:00 2001 From: SNSTATS <65264926+snstats@users.noreply.github.com> Date: Wed, 28 Feb 2024 21:52:47 +1100 Subject: [PATCH 1/2] looking into dsum variogram --- R/vario for dsum - seeing how it could work.r | 156 ++++++++++++++++++ 1 file changed, 156 insertions(+) create mode 100644 R/vario for dsum - seeing how it could work.r diff --git a/R/vario for dsum - seeing how it could work.r b/R/vario for dsum - seeing how it could work.r new file mode 100644 index 0000000..0007399 --- /dev/null +++ b/R/vario for dsum - seeing how it could work.r @@ -0,0 +1,156 @@ +library(ggplot2) +library(asreml) +library(biometryassist) +library(tidyverse) +library(lubridate) +library(viridis) + +theme_set(theme_bw()) + + + +datn <- read.csv("C:/Users/sharo/Documents/SNSTATS/Clients/Daniel Menadue/Objective 2 Hart/HartNDVI.csv") +str(datn) +datn$Treatment <- factor(datn$Treatment, levels = c("Control", + "Balance&Grow", "GABA", "Kelpak", "RemoveFlorets", "Seasol")) + +datn <- datn %>% mutate(across(c(1:5), factor)) + +datn$Date <- dmy(datn$Date) + +datn <- datn %>% arrange(Date, Row, Bay) + +datn$facDate <- factor(datn$Date) +datn$plot <- factor(paste(datn$Row, datn$Bay, sep = "")) + +dat.asr <- asreml(fixed = NDVI ~ facDate + Variety + Variety:facDate + Treatment + + Variety:Treatment + Treatment:facDate + + Variety:Treatment:facDate, +random = ~ Block, +residual = ~ dsum(~ ar1(Row):id(Bay)|facDate), data = datn) + + + +dat.asr2d <- asreml(fixed = Biomass ~ Variety + Treatment + + Variety:Treatment, +random = ~ Block, +residual = ~ ar1(Row):id(Bay), data = dat) +dat.asr2d <- update(dat.asr2d) + +Row <- "Row" +Coulmn <- "Bay" +model.obj <- dat.asr2d + +######################################################################### + +# vario_df <- function(model.obj, Row = NA, Column = NA) { +# # The 'z' value for the variogram is the residuals +# # Need to be able to pull out the x/y from the model object + +# This works for 2 dimensions + +dims <- unlist(strsplit(names(model.obj$R.param[1]), ":")) + +dims <- unlist(strsplit(names(dat.asr2d$R.param[1]), ":")) + + +# This will flag that a dsum is used in the model - should do this first then +# do the 2D case + +# now to identify the terms in the residual - the lasst will be the factor +# associated with the last term in the res.terms vector + + +# Once dsum is identified - the dat.asr$formulae$residual will need to be broken +# down into it's components. +# Seperate data set for each |term and row column within each |term + +# pulling out the residual formulae as a character +res.call <- as.character(dat.asr$formulae$residual)[2] + +# identifying if dsum is present +# if present pull out residual terms +if(any(grepl("dsum", res.call))){ +res.terms <- unlist(strsplit(res.call, "[~:| ()]+"))[-1] +asr.terms <- c("id", "ar1", "ar2", "cor") + +res.terms <- res.terms[!is.element(res.terms, asr.terms)] +} +res.terms + +# Break up the data and the residuals into all levels of grouping factor +# Calculate the gamas for each level and create a separate variogram for each +# level + +# Maybe have the option to combine into a single graph or keep as separate plots + + +# This is for 2D + + if(missing(Row) | is.na(Row) | is.null(Row)) { + Row <- as.numeric(model.obj$mf[[dims[1]]]) + } + if(missing(Column) | is.na(Column) | is.null(Column)) { + Column <- as.numeric(model.obj$mf[[dims[2]]]) + } + + nrows <- max(Row) + ncols <- max(Column) + + Resid <- residuals(model.obj)#[order(Column, Row)], nrow = nrows) + # Resid <- matrix(residuals(model.obj)[order(Column, Row)], nrow = nrows) + + vario <- expand.grid(Row = 0:(nrows-1), Column = 0:(ncols-1)) + + # Ignore the 0, 0 case (gamma=0, counted row*cols times) + gammas <- rep(0, nrows*ncols) + nps <- rep(nrows*ncols, nrows*ncols) + + for (index in 2:nrow(vario)) { + i <- vario[index, 'Row'] + j <- vario[index, 'Column'] + + gamma <- 0 + np <- 0 + for (val_index in 1:nrow(vario)) { + # val <- vals[val_index, ] + + # Deliberate double-counting so that offset handling is easy + # (so e.g. we compute distance from (1,1)->(2,3), and then again + # later from (2,3)->(1,1)). + for (offset in unique(list(c(i, j), c(-i, j), c(i, -j), c(-i, -j)))) { + row <- Row[val_index] + offset[1] + col <- Column[val_index] + offset[2] + + if (0 < row && row <= nrows && 0 < col && col <= ncols && !is.na(Resid[val_index])) { + other <- Resid[Row == row & Column == col] + + if (!is.na(other)) { + gamma <- gamma + (Resid[val_index] - other)^2 + np <- np + 1 + } + } + } + } + # Since we double-counted precisely, halve to get the correct answer. + np <- np / 2 + gamma <- gamma / 2 + + if (np > 0) { + gamma <- gamma / (2*np) + } + + gammas[index] <- gamma + nps[index] <- np + } + nps[1] <- nps[1]-sum(is.na(Resid)) + vario <- cbind(vario, data.frame(gamma = gammas, np = nps)) + colnames(vario) <- c(dims, "gamma", "np") + class(vario) <- c("variogram", "data.frame") + return(vario) +} + + +vario_df(dat.asr, Row = "Row", Column = "Bay") + + From c3955cb4e4ea8ef97d7b326c6ce1047bded63a68 Mon Sep 17 00:00:00 2001 From: SNSTATS <65264926+snstats@users.noreply.github.com> Date: Wed, 28 Feb 2024 21:56:24 +1100 Subject: [PATCH 2/2] fixing position of files in folder --- ...it could work.r => vario for dsum - seeing how it could work.r | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename R/vario for dsum - seeing how it could work.r => vario for dsum - seeing how it could work.r (100%) diff --git a/R/vario for dsum - seeing how it could work.r b/vario for dsum - seeing how it could work.r similarity index 100% rename from R/vario for dsum - seeing how it could work.r rename to vario for dsum - seeing how it could work.r