0

I have a dataset on fecal egg counts (FEC) in sheep and would like to test the effect of different treatments. I have 120 farms treated with combination of different drugs. I am using the "eggcount" package that is based on the Bayesian approach to simulate and produce the outputs for each farm and treatment (Drug).

I want to make a loop for two variables; i) Farm and ii) Drug. This means the for loop can generate the output for each farm and drug, similar to the subset function.

setwd("C:/Data")
dat <- read.csv("McMaster_UTC.csv", header=T, sep="'", na.string=NA)
dat$EPG1 <- as.numeric(dat$EPG1)
dat$EPG2 <- as.numeric(dat$EPG2)

# subset for Fram & Drug: if I want to run the data for each Farm and Drug. I assume with a loop I don't #need the subset function.
dat1 <-subset(dat, Farm==1 & Drug=="CLO", select=1:9)

# A model to compare the two groups (unpaired)
model <- fecr_stan(egg.1$EPG1/50, egg.1$EPG2/50, rawCounts=FALSE, paired=FALSE, indEfficacy=FALSE)

#This is the for loop that I have in my mind, but not sure how to fit these in the loop:
out <- vectro()
for(i in 1:ncol(dat)) {
  for(j in 1:ncol(dat)) {
....
....
set.seed(1234)
model <- fecr_stan(egg.1$EPG1/50, egg.1$EPG2/50, rawCounts=FALSE, paired=FALSE, indEfficacy=FALSE)

    model[i] <- out[i]

Here is a subset of my data:

Sample data for 2 farms with two drugs

Ahmad
  • 3
  • 3

1 Answers1

0

If I understand what you want correctly, you just want to loop over each combination of Farm and Drug. Step 1 would be to work out what the levels of each of these factors are. Given you have read your CSV file without stringsAsFactors set to FALSE, we can get the levels using the levels function:

Farms = levels(dat$Farm)
Drugs = levels(dat$Drug)

Then to loop over each, you can just do:

set.seed(1234) ## Only put this inside your loop if you want to use the same
               ## sequence of random numbers for each model.  
for(farm in Farms){
  for(drug in Drugs){
    dat1 <-subset(dat, Farm == farm & Drug == drug, select=1:9)
    model <- fecr_stan(dat1$EPG1/50, dat1$EPG2/50, 
                       rawCounts = TRUE, paired=FALSE, 
                       indEfficacy=FALSE)

   ## Note you'll have to do something to save or store the model output here
  }
}

Thanks James for your codes/comments, much appreciated. First; there are a couple of errors in my codes: egg.1$EPG1/50 and egg.1$EPG1/50 should be dat$EPG1/50 & dat$EPG1/50. I must say that I am not good with loops, just in case if I am asking silly questions here. Yes, I want to compare the difference between dat$EPG1/50 & dat$EPG2/50 for each Drug within each Farm, and save the output in a vector. When you say "loop over each", what do you mean by that? Are you suggesting that I put model <- fecr_stan(dat$EPG1/50, dat$EPG2/50, rawCounts=TRUE, paired=FALSE, indEfficacy=FALSE), and it will work? Should I put the setseed(12345) outside the loop? I appreciate having your thoughts! I attached a sample of my data as an image, if you want to have a look. I didn't know how I can paste a sample of my data here.

James Curran
  • 1,274
  • 7
  • 23
  • @Ahmad - does it do what you want to do? – James Curran Apr 10 '20 at 03:04
  • @Ahmad - I've altered the code to include some of the steps. `for` loops in R are more of the "iterate over elements of a list/vector" type (and are sometimes called `foreach` loops in other languages), then traditional "increment a counter" type `for` loops. So `for(drug in Drugs)` will step through each of the levels of the factor which are stored in `Drugs` etc. – James Curran Apr 10 '20 at 03:13
  • Thanks, James, I used your codes in the loop. I created a vector ```out <- NA``` (outside the loop) to save the output and put this ```out[farm, drug] <- model``` inside the loop as you suggested. It worked ok (no error)- but didn't produce an output to be saved in the vector ```out```. Your thoughts? I will keep working on it to see how I can make it work. – Ahmad Apr 10 '20 at 08:00
  • I also noticed that ```Farms = levels(dat$Farm)``` didn't work as expcted and generated ```NULL``` ?? – Ahmad Apr 10 '20 at 08:13
  • Hi @Ahmad. Try typing `is.factor(dat$Farm)` after you have read the data in. If it returns `FALSE`, then you have two options. Either turn it into a factor `dat$Farm = factor(dat$Farm)`, or instead of using the levels,just use the unique values of `Farm`, i.e. `Farms = sort(unique(dat$Farm))` – James Curran Apr 10 '20 at 15:08
  • Thanks James, I noticed even with Fram as a numeric variable in the data, I had to convert it to a factor. I will try this to see if this will generate an output. Thanks again! – Ahmad Apr 12 '20 at 05:16