Disabled.R

monsterlist <- c() #This will hold the number of disabled in each sim run

for (h in 1:100){ #Repeat simulation 100 times


disabled <- c() #This will hold the disabled within each sim run

population <- as.vector(c(rep(0,10000))) #The population of 10,000


for(j in 1:30) { # Number of Years

infected <- sample(1:length(population),floor(0.1*length(population))) # infect x% of the population every year


population[infected] <- population[infected]-1 #infected get an increasingly-negative counter

# The next 2 paragraphs control who gets LC. One should be commented out.
# The first puts everyone in the same risk category.
# The second creates low, medium, and high risk categories.


# for (k in 1:length(population)){ #Use this paragraph if risk to everyone of LC is equal. 
#Use next paragraph if LC is based on low risk, medium risk, and high risk groups.
#   #if (runif(1) <= pweibull(abs(population[k]),scale=6,shape=2)){ #Use this one if successive Covid infections increase LC risk
#   if (runif(1) <= 0.1) { #Use this one if risk of LC is constant. Comment out one of these.
#     population[k] <- population[k]+ 100 #positive number means Long Covid
#   }
#   
# }



for (p in 1:length(population)){
  if (p <= 3000) { #The first 3000 in the population are "low risk" for LC. 
    if(runif(1) <= 0.05) { #Only a 5% chance of getting LC with each infection.  
    population[p] <- population[p] + 100
    }
    
  } else if (p >= 7000) { #High risk group for LC.
      if(runif(1) <= 0.25) { #This group has a 25% chance of getting LC with each infection.
        population[p] <- population[p] + 100
      }
  } else { #This is the "medium" risk group of getting LC with each infection.
      if(runif(1) <= 0.1) { #This group has a 10% chance of getting LC with each infection.
        population[p] <- population[p] + 100
      }
    }
}

for (n in 1:length(population)){ # give them a 50% chance to recover from LC this year
  if (population[n] > 0){
    if(runif(1)>=0.5){
      population[n] <- population[n] - 100
      
    }
  }
}

print(j)

del_list <- c()

for (m in 1:length(population)){
  if (population[m] > 0){
    if(runif(1)>=0.9){
      disabled <- c(disabled,population[m])
      del_list <- c(del_list,m)
      
    }
  }
}
if (length(del_list)>0){

population <- population[-del_list]
}






}



monsterlist <- c(monsterlist,length(disabled))

}

sum(monsterlist)/10000 #gives the % (not the decimal) average over 100 simulations

Hello,
you created a really nice model. One suggestion for improvement:
You can simply create the population vector with the code 'population <- rep(0,10000)'
There's no need to use the as.vector-function (or even the c), because the rep-function creates a vector by default.
Best Wishes