minor edits and including script that generated numbers

_posts/2024-03-18-Hypothesis-why-do-neutralizing-antibody-titers-max-out.md

@@ -72,7 +72,7 @@ Equating terms, we see first that we should have
 
 $$\frac{\mu_0}{\log_2(\text{max titer})}=\gamma \tag{4}$$
 
-relating our model derived from antibody response and our model derived from susceptibility to infection. The remarkable thing is $\gamma=0.46$ and $\frac{\mu_0}{\log_2(\text{max titer})}=0.48$!  I've been lazy and not showing uncertainty intervals, but this is damn near exactly the same, and very much not significantly different!
+relating our model derived from antibody response and our model derived from susceptibility to infection. The remarkable thing is $\gamma=0.46$ and $\frac{\mu_0}{\log_2(\text{max titer})}=0.48$!  I've been lazy and not showing uncertainty intervals, but this is damn near exactly the same, and very much not significantly different with the intervals that I have!
 
 **I want you to share in my excitement. We get the exact same model, with the exact same numbers, from two very different places.  The first is from just directly measuring adaptive immune response.  The second is from fitting a model to human challenge studies for the probability a person gets infected from a range of oral doses,[^4] taking seriously the biological interpretation of the model in terms of the human infectiousness per cell culture infectious unit, and *finding empirically via this agreement* that pre-challenge neutralizing antibodies have *exactly the same effect* on both getting infected and building new responses to prevent the next infection. And the *assumption that ties them together is that this all is determined by what happens during the viral replication and propagation process at the cellular level*!**
 
@@ -105,11 +105,13 @@ For attribution, please cite this work as:
 
 `Famulare (2024, Mar 18). Hypothesis---why do neutralizing antibody titers max out? Retrieved from https://famulare.github.io/2024/03/18/Hypothesis-why-do-neutralizing-antibody-titers-max-out.html.`
 
+The script used to generate various numbers in this post can be found [here](/assets/2024-03-18-Hypothesis-why-do-neutralizing-antibody-titers-max-out/linkages_among_polio_immune_model_components.R).
+
 ____
 [^1]: but unfortunately the numbers aren't directly comparable because it's based on hemagglutinin assay inhibition titers and not neutralizing antibody titers. Different reference standard.
 
-[^2]: To my eternal shame, we fit the mean and variance separately in the paper even though we knew they should've been coupled. In my defence, I was very frazzled at the time by intense overwork on the [Seattle Flu Study](https://seattleflu.org/history) and Wes had no help debugging what I now suspect was a subtle convergence error in the fitting script when he tried to do it right. The figures and parameters mentioned in this post differ from the paper because I redid the analysis correctly here.
+[^2]: To my eternal shame, we fit the mean and variance separately in the paper even though we knew they should've been coupled. In my defence, I was very frazzled at the time by intense overwork on the [Seattle Flu Study](https://seattleflu.org/history) and Wes had no help debugging what I now suspect was a subtle convergence error in the fitting script when he tried to do it right. The figures and parameters mentioned in this post differ from the paper because I redid the analysis correctly here. 
 
 [^3]: except you might notice that at very high titer, some of the non-responder dots are inside the range of expected responses. At very high titer, responses are small enough that you can't tell from seroresponse if someone was infected in the interval or not. This should be accounted for in how I fit the model, but I was too lazy to do so. Which also means the max titer is likely biased a bit higher than the real value. Or closer to $2^{14}$ :-). 
 
-[^4]: As an aside, I started my deep dive because I couldn't make sense of the susceptibility dose response part of the original model in the [Behrend paper](https://www.ijidonline.com/article/S1201-9712(13)00297-X/fulltext), and so I re-did it my own way (and found the original was wrong but for reasonable reasons). In going down that road, I since redeveloped a bunch of other stuff, but I periodically find myself surprised that Matthew already thought and wrote about it and I forgot!
+[^4]: As an aside, I started my deep dive because I couldn't make sense of the susceptibility dose response part of the original model in the [Behrend paper](https://www.ijidonline.com/article/S1201-9712(13)00297-X/fulltext), and so I re-did it my own way (and found the original was wrong but for reasonable reasons). In going down that road, I've since redeveloped a bunch of other stuff, but I periodically find myself surprised at everything Matthew thought about, wrote about, got right, and I forgot!

assets/2024-03-18-Hypothesis-why-do-neutralizing-antibody-titers-max-out/linkages_among_polio_immune_model_components.R

@@ -0,0 +1,144 @@
+# scratch script for 2024-03-18-Hypothesis-why-do-neutralizing-antibody-titers-max-out.md
+
+library(tidyverse)
+
+# polio dose response model
+alpha=0.44
+beta=2.3
+gamma=0.46
+
+rbar = function(NAb=1){
+  1/(1+beta/alpha*NAb^gamma)
+}
+
+p = function(NAb=1,dose=1){
+  1-(1+dose/beta)^(-(alpha/NAb^gamma))
+}
+
+dens = function(x,NAb=1){
+  dbeta(x,shape1=alpha/(NAb^gamma),shape2=beta)
+}
+
+rand = function(dose,NAb=1){
+  rbeta(dose,shape1=alpha/(NAb^gamma),shape2=beta)
+}
+
+1*rbar(NAb=2^14 )
+
+plot(seq(0,14),1/rbar(NAb=2^(0:14)))
+
+plot(0:14,p(dose=1000,NAb=2^(0:14)))
+
+
+
+x=seq(0,1,by=1e-4)
+plot(x,dens(x,NAb=2^14))
+
+
+p(NAb=2^14,dose=1e6)
+
+
+# median from the shedding duration model
+log(30.3)-log(1.16)*log2(2^14)
+
+
+
+## LOL did I misremember 2^14 for the polio model?!
+
+temp_file <- tempfile(fileext = ".xlsx")
+download.file(url = 'https://github.com/famulare/cessationStability/raw/master/data/Louisiana1957/householdSeroconversion.xlsx', destfile = temp_file, mode = "wb", quiet = TRUE)
+
+df=readxl::read_xlsx(temp_file,sheet='prePostAntibodyLineList')[,1:2] |>
+  mutate(seroconverted = (postExposureTiter/preExposureTiter)>=4)
+df
+
+ggplot(df)+
+  geom_jitter(aes(x=preExposureTiter,y=postExposureTiter/preExposureTiter,color=seroconverted)) +
+  scale_x_continuous(trans='log2',breaks=2^(3*(0:5)))+
+  scale_y_continuous(trans='log2',breaks=2^(3*(0:5)))
+
+responseModelLik = function(params=c(4.8,3.3,-0.31/4.4,-0.09,0.36),data=df){
+  Msamp=params[4]
+  Ssamp=params[5]
+
+  logL=0
+
+  for(k in 1:nrow(data)){
+
+    x=log2(data$postExposureTiter[k])-log2(data$preExposureTiter[k])
+
+    if(df$seroconverted[k]==TRUE){
+
+      Msc=params[1]*(1+params[3]*log2(data$preExposureTiter[k])) + Msamp
+      Ssc=sqrt( (params[2]*(1+params[3]*log2(data$preExposureTiter[k])))^2 + Ssamp)
+
+      # df$logL[k]=dnorm(x=x,mean=Msc,sd=Ssc,log=TRUE)
+
+      logL = logL + dnorm(x=x,mean=Msc,sd=Ssc,log=TRUE)
+
+    } else {
+      logL = logL + dnorm(x=x,mean=Msamp,sd=Ssamp,log=TRUE)
+
+      # df$logL[k]=dnorm(x=x,mean=Msamp,sd=Ssamp,log=TRUE)
+
+    }
+
+  }
+  return(logL)
+  # return(df)
+}
+
+
+res=stats4::mle(function(x=c(4.8,3.3,-0.31/4.4,-0.09,0.36)){responseModelLik(params=x,data=df)},
+          method='L-BFGS-B',
+          lower=c(0,0,-1/12.34,-Inf,0.01),
+          upper=c(Inf,Inf,0,Inf,Inf),
+          control=list(fnscale=-1,trace=0))
+
+res
+-1/res@coef[3]
+# profile likelihood
+# res_ci=stats4::confint(res)
+
+-res@coef[1]*res@coef[3]
+gamma
+
+
+C=2^(res@coef[1])/900*beta/alpha
+C
+
+1/C
+
+
+# index of dispersion, not unlike Poisson process...
+res@coef[2]^2/res@coef[1]
+
+
+
+
+# matthew reference 
+# table 5 https://www.ijidonline.com/action/showFullTableHTML?isHtml=true&tableId=tbl0025&pii=S1201-9712%2813%2900297-X 
+# mean of max titers = 13.5
+mean(c(11.1,11.2,15.3,12.9,15.5,15.2))
+
+
+pred=data.frame(preExposureTiter=2^seq(0,-1/res@coef[3],by=0.01)) |>
+  mutate(median=exp(res@coef[1]*(1+res@coef[3]*log2(preExposureTiter)))) |>
+  mutate(lower=exp(res@coef[1]*(1+res@coef[3]*log2(preExposureTiter))-(res@coef[2]*(1+res@coef[3]*log2(preExposureTiter))))) |>
+  mutate(upper=exp(res@coef[1]*(1+res@coef[3]*log2(preExposureTiter))+(res@coef[2]*(1+res@coef[3]*log2(preExposureTiter))))) |>
+  mutate(lower2=exp(res@coef[1]*(1+res@coef[3]*log2(preExposureTiter))-2*(res@coef[2]*(1+res@coef[3]*log2(preExposureTiter))))) |>
+  mutate(upper2=exp(res@coef[1]*(1+res@coef[3]*log2(preExposureTiter))+2*(res@coef[2]*(1+res@coef[3]*log2(preExposureTiter)))))
+
+ggplot(df)+
+  geom_jitter(aes(x=log2(preExposureTiter),y=postExposureTiter/preExposureTiter,color=seroconverted)) +
+  geom_line(data=pred,aes(x=log2(preExposureTiter),y=median)) +
+  geom_line(data=pred,aes(x=log2(preExposureTiter),y=median)) +
+  geom_ribbon(data=pred,aes(x=log2(preExposureTiter),ymin=lower,ymax=upper),alpha=0.1)+
+  geom_ribbon(data=pred,aes(x=log2(preExposureTiter),ymin=lower2,ymax=upper2),alpha=0.1)+
+  scale_x_continuous(breaks=(3*(0:5)))+
+  scale_y_continuous(trans='log2',breaks=2^(2*(0:5))) 
+
+
+## venom
+delta_mu = 1/3*log(8.4e6/15e3)
+delta_mu