Tutorial from Amidlige and MBuynak

Tutorial 9

@@ -0,0 +1,80 @@
+#Import Packages
+import pandas
+import numpy
+from scipy.optimize import minimize
+from scipy.stats import norm
+import re
+import os
+from plotnine import *
+
+os.chdir("/Users/madelinebuynak/Desktop")
+
+#load data
+data = pandas.read_csv("ponzr1.csv",header=0,sep='\t')
+
+subset=data.loc[data.mutation.isin(['WT','I231N']),:]
+
+#two likelihood functions
+def nllike(p,obs):
+    B0=p(0)
+    B1=p(1)
+    sigma=p(2)
+
+def null (p,obs):
+    B0=p[0]
+    sigma=p[1]
+
+    expected=B0
+    nll=-1*norm(expected,sigma).logpdf(obs.y).sum()
+    return nll
+
+def alter (p,obs):
+    B0=p[0]
+    B1=p[1]
+    sigma=p[2]
+
+expected=B0+B1*obs.x
+nll=-1*norm(expected,sigma).logpdf(obs.y).sum()
+return nll
+
+initialGuess=numpy.array([1,1,1])
+fitNull=minimize(nllike,initialGuess,method="Nelder-Mead",options={'disp':True},args=subset)
+fitAlter=minimize(nllike,initialGuess,method="Nelder-Mead",options={'disp':True},args=subset)
+
+
+subset=data.loc[data.mutation.sin(['WT','V456D']),:]
+
+
+def null (p,obs):
+    B0=p[0]
+    sigma=p[1]
+
+    expected=B0
+    nll=-1*norm(expected,sigma).logpdf(obs.y).sum()
+    return nll
+
+def alter (p,obs):
+    B0=p[0]
+    B1=p[1]
+    sigma=p[2]
+
+    expected=B0+B1*obs.x
+    nll=-1*norm(expected,sigma).logpdf(obs.y).sum()
+    return nll
+
+initialGuess=numpy.array([1,1,1])
+fitNull=minimize(null,initialGuess,method="Nelder-Mead",options={'disp':True},args=subset)
+fitAlter=minimize(alter,initialGuess,method="Nelder-Mead",options={'disp':True},args=subset)
+
+D=2*(fitNull.fun-fitAlter.fun)
+1-ch2.cdf(x=D,df=1)
+
+#M124K: p-value ~ 0.72 (no effect of treatment)
+
+#V456D: p-value ~ 5.6e-6 (effect of treatment)
+
+#I213N: p-value ~ 0.88 (no effect of treatment)
+
+
+
+

Tutorial 9 Question 3

@@ -0,0 +1,72 @@
+#Import Packages
+import pandas
+import numpy
+from scipy.optimize import minimize
+from scipy.stats import norm
+import re
+import os
+from plotnine import *
+
+os.chdir("/Users/madelinebuynak/Desktop")
+
+#load data
+data=pandas.read_csv("leafDecomp.csv",header=0,sep='\t')
+
+subset=data.loc[data.mutation.isin(['WT','I231N']),:]
+
+
+#plot data
+plot = ggplot(data, aes(y = 'Ms', x = 'decomp'))
+plot + geom_jitter(colour='black')
+
+#quadratic
+
+#define 3 custom likelihood functions
+def nllike(p,obs):
+    B0=p(0)
+    B1=p(1)
+    sigma=p(2)
+
+def null (p,obs):
+    B0=p[0]
+    sigma=p[1]
+
+    expected=B0
+    nll=-1*norm(expected,sigma).logpdf(obs.y).sum()
+    return nll
+
+def alter (p,obs):
+    B0=p[0]
+    B1=p[1]
+    sigma=p[2]
+
+    expected=B0+B1*obs.x
+    nll=-1*norm(expected,sigma).logpdf(obs.y).sum()
+    return nll
+
+def alter (p,obs):
+    B0=p[0]
+    B1=p[1]
+    sigma=p[3]
+
+    expected=B0+B1*obs.x+B2*obs.x^2
+    nll=-1*norm(expected,sigma).logpdf(obs.y).sum()
+    return nll
+
+#Estimate the Parameters
+#Initial Guess for Model 1 (mean~590)
+#Initial Guess for Model 2 (intercept~200 slope~6.33)
+
+initialGuess=numpy.array([1,1,1])
+fitNull=minimize(nllike,initialGuess,method="Nelder-Mead",options={'disp':True},args=subset)
+fitAlter=minimize(nllike,initialGuess,method="Nelder-Mead",options={'disp':True},args=subset)
+
+#Compare Models
+#likelihood test
+D=2*(fitNull.fun-fitAlter.fun)
+1-ch2.cdf(x=D,df=1)
+
+
+#quadratic model is the best fit, but linear model is better than constant model
+
+

tutorial9question2.py

@@ -0,0 +1,26 @@
+#Import Packages
+import pandas
+import numpy
+from scipy.optimize import minimize
+from scipy.stats import norm
+import re
+import os
+from plotnine import *
+
+data = pandas.read_csv("mMarinumGrowth.csv",header=0,sep='\t')
+
+def Monod(p,obs):
+    uMax=p[0]
+    K=p[1]
+    sigma=p[2]
+
+    expected=uMax* (x/(x+K))
+    nll=-1+norm(expected,sigma).logpdf(obs.y).sum
+    return nll
+
+    initialGuess=numpy.array([1,1,1])
+    fit-minimize(Monod,initialGuess,method="Monod Equation",options={'disp':True},args=data)
+
+    print(fit.x)
+
+