Formålet med dette afsnit er at gøre jer bekendt med den statistiske model, der ligger bag avancerede programmer i python til estimation af ikke-lineære sammenhænge, specifikt pythons funktion curve$\text{\textunderscore}$fit. I skal ikke selv analysere data med disse programmer, da dette ligger uden for indholdet i denne webbog. I nogle situationer er en lineær funktion ikke tilstrækkelig til at beskrive sammenhængen i data. Et eksempel er inden for studiet af iontransport. I artiklen Unbiased Simulations Reveal the Inward-Facing Conformation of the Human Serotonin Transporter and Na+ Ion Release, skrevet af forskere fra Aarhus Universitet, omtales et eksperiment med transport af Na+ i celler fra pattedyr. Forskellige koncentrationer af NA+ er indsat i cellerne, og raten hvormed ionerne transporteres måles. Der omtales to eksperimenter, hvor der bruges forskellige proteiner til transport af Na+, henholdsvis en vildtype hSERT (human serotonin transporter) og en muteret type Asp437Asn. Koncentrationen måles i mM og raten i cpm (counts per minute). Data der anvendes nedenfor er aflæst fra figur 6 i artiklen og er vist i nedenstående figuren. Michaelis-Menten kinetik forbinder en reaktionsrate $v$ med koncentrationen $T$ af et substrat gennem ligningen $$v(T)=\frac{\alpha T}{\beta+T}. \tag{7.8.1}$$ Der er således to parametre i denne model. Parameteren $\alpha$ angiver den maksimale reaktionsrate, og parameteren $\beta$ angiver den koncentration, under hvilken reaktionsraten vil være det halve af den maksimale værdi $\alpha$. Relationen (7.8.1) med skønnede parameterværdier er indtegnet i figuren nedenfor, hvor den venstre kurve er for hSERT eksperimentet og den højre kurve er for Asp437Asn eksperimentet. Den generelle model vi vil betragte beskriver den ikke-lineære sammenhæng gennem parametre i vektoren $\theta.$ Typisk vil der også være en forklarende variabel $t,$ således at den $i$'te observation af respons er knyttet til værdien $t_i$ af den forklarende variabel. Vi kan skrive dette generelt ved at lade middelværdien af det $i$'te respons $X_i$ være $E(X_i)=g_i(\theta)$. Den lineære regressionsmodel vi har betragtet i dette kapitel svarer til $g_i(\alpha,\beta)=\alpha+\beta t_i,$ og Michaelis-Menten kinetik omtalt ovenfor svarer til $g_i(\alpha,\beta)=\frac{\alpha t_i}{\beta+t_i}.$ Hvor vi i den lineære regressionsmodel forlanger at variansen af respons $X_i$ er den samme for alle $i,$ vil vi her tillade at variansen skalerer med kendte vægte, $\text{Var}(X_i)=\sigma^2w_i,$, hvor $(w_1,\ldots,w_n)$ er de kendte vægte. I ovenstående figur er vægtene illustreret gennem en errorbar ved hvert målepunkt. Som det sidste element i modellen siger vi at data er normalfordelte. Skøn over $\theta$ findes ved vægtet mindste kvadraters metode svarende til at vi finder skønnet $\hat\theta$ ved at minimere $$D(\theta)=\sum_{i=1}^n\frac{1}{w_i}\big(x_i-g_i(\theta)\big)^2.$$ Når skønnet $\hat\theta$ er fundet bruger vi som skøn over $\sigma^2$ værdien $$s^2=\frac{D(\hat\theta)}{n-d}= \frac{1}{n-d}\sum_{i=1}^n\frac{1}{w_i}\big(x_i-g_i(\hat\theta)\big)^2,$$ hvor $d$ er antallet af parametre i $\theta.$ Der gælder approksimativt, at $$\begin{aligned} s^2 & \sim \sigma^2\chi^2(n-d)/(n-d), \\ t & =\frac{\hat\theta_j-\theta_j}{\text{std}_s(\hat\theta_j)} \sim t(n-d), \end{aligned}$$ hvor det sidste bruges til at lave et approksimativt 95%-konfidensinterval på formen $$\hat\theta_j\pm t_0\text{std}_s(\hat\theta_j),\quad t_0=t_{\text{inv}}(0.975,n-d).$$ Som en del af output fra curve$\text{\textunderscore}$fit kan man få en kovariansmatriks for $\hat\theta.$ Dette er en $d\times d$ matriks, hvor det $j$'te diagonalelement er $\text{std}_s(\hat\theta_i)^2$ og det $(j,k)$'te element er $\text{Cov}_s(\hat\theta_j,\hat\theta_k).$ Med denne matriks til rådighed kan vi bruge ophobningsloven til beregning af usikkerheden på afledede parametre. Formlen for matricen er $$s^2 J^{-1},\quad J_{jk}=\sum_{i=1}^n\frac{1}{w_i} \frac{\partial g_i}{\partial \theta_j}(\hat\theta) \frac{\partial g_i}{\partial \theta_k}(\hat\theta).$$ Som nævnt er fordelingsresultaterne (og konfidensintervallerne baseret på disse) approksimative, og erfaringsmæssigt ved man, at approksimationerne ikke altid er gode. Nedenfor analyseres de to datasæt vist i figuren tidligere i dette afsnit omkring Michaelis-Menten kinetik. Modellen er således $$X_i\sim N\Big(\frac{\alpha t_i}{\beta+t_i},\frac{\sigma^2}{w_i}\Big), \enspace i=1,\ldots,n, \tag{7.8.2}$$ hvor $t_i$ er koncentrationen af Na+ ioner, og $X_i$ er den målte transportrate. Pythons funktion curve$\text{\textunderscore}$fit skal som input have en funktion, der definerer den ikke-lineære sammenhæng, en vektor med værdierne af den forklarende variabel, en vektor med responsværdierne, eventuelt en vektor med vægte (tildeles sigma), og eventuelt en vektor med startværdier for søgning efter skøn over parametrene (tildeles p0). Ofte vil søgerutinen ikke kunne finde passende skøn over parametrene medmindre den hjælpes på vej med startværdier. Output fra curve$\text{\textunderscore}$fit består af parameterskøn (popt i kode nedenfor) og kovariansmatriks beskrevet ovenfor (pcov i kode nedenfor).

Se opstartskoden (til/fra)

Kørsel af koden ovenfor giver følgende skøn og 95%-konfidensintervaller for $\alpha$ og $\beta$. $$\begin{array}{lcccccccc} \hline && &\alpha & && & \beta & \\ \text{Protein} && \text{Skøn} & \text{Nedre} & \text{Øvre} && \text{Skøn} & \text{Nedre} & \text{Øvre} \\ \hline \text{hSERT} && 21678 & 20453 & 22902 && 2.32 & 1.06 & 3.59 \\ \text{Asp437Asn} && 72553 & 33437 & 111669 && 367 & 117 & 616 \\ \hline \end{array}$$ For begge datasæt får vi et meget bredt konfidensinterval (den øvre grænse er flere gange større end den nedre grænse). For ASP datasættet (til højre i figur) er det tydeligt, at forklaringen på dette er, at vi ikke er i stand til at fastlægge den øvre grænse for raten, og dermed kan vi ikke fastlægge hvilken koncentration, der giver en rate, der er halvdelen af maksimum. For hSert datasættet har vi god viden om den øvre grænse, men i det område hvor raten cirka er halvdelen af den øvre grænse, er der stor usikkerhed på ratemålingen, hvorfor vi har svært ved at fastlægge, for hvilket koncentration raten er lig med halvdelen af den øvre grænse. ForegåendeNæste   

• Øvelse 5:
• Kemi
• Fysik