I denne opgave skal I se på data fra artiklen Prediction of full load electrical power output of a base load operated combined cycle power plant using machine learning methods, hvor man ønsker at prædiktere Full Load Electrical Power Output (PE målt i mega watt) time for time baseret på fire forklarende variable. De fire forklarende variable er Ambient Temperature (AT målt i grader celsius), Vacuum (Exhaust Steam Pressure) (V målt i cm Hg), Atmospheric Pressure (AP målt i mbar) og Relative Humidity (RH målt i procent). Datafilen PowerOutput.csv har søjlerne AT, V, AP, RH, PE, men der er også tilføjet følgende søjler, $$\mathit{AT2}=\mathit{AT}^{\,2},\quad \mathit{AT3}=\mathit{AT}^{\,3},\quad \mathit{V2}=\mathit{V}^{\,2},\quad \mathit{AP2}=\mathit{AP}^{\,2},\quad \mathit{ATRH}=\mathit{AT}\cdot\mathit{RH}.$$

1. Indlæs data som en datatabel.
   Opskriv regressionsmodellen, hvor middelværdien af PE afhænger lineært af de fire forklarende variable AT, V, AP og RH.
2. Lav en figur, hvor residualerne i modellen afsættes mod AT. Hvad ser du i denne figur (to ting) ?
3. Betragt nu regressionsmodellen med de seks forklarende variable AT, V, AP, RH, AT2 og AT3. Lav igen en figur, hvor residualerne i modellen afsættes mod AT.
4. Betragt nu regressionsmodellen med alle ni forklarende variable AT, V, AP, RH, AT2, AT3, V2, AP2 og ATRH.
   Reducer den fulde multiple regressionsmodel ved successivt at fjerne led i modellen (backward selektion). Lav en tabel, som for hvert trin angiver den variabel der fjernes i modellen og den tilhørende $p$-værdi.
   Hvor meget stiger power output med når luftfugtigheden falder fra 80 til 60 procent ?
5. Opskriv slutmodellen ved backward selektionsproceduren, og lav et $F$-test for reduktion fra den fulde model til slutmodellen.
   Lav et qqplot af residualerne i slutmodellen. Angiv Resultat nummer, for hvordan et qqplot vurderes
6. Angiv et 95%-konfidensinterval for regressionsparameteren, der ganges på den forklarende variabel $RH$ (luftfugtighed) i slutmodellen.
7. Det oprindelige datasæt har 9568 målinger, hvoraf jeg har givet jer 5000 ovenfor til at vælge en model. I filen PowerOutputTest.csv er de resterende 4568 målinger.
   Indlæs data fra PowerOutputTest.csv som en datatabel. Benyt den estimerede model ovenfor til at prædiktere PE for de nye data fra PowerOutputTest.csv.
   Lav en figur, hvor forskel mellem den målte og prædikterede værdi af PE afsættes mod AT.
   Udregn kvadratroden af den gennemsnitlige værdi af den kvadrerede forskel mellem den målte og prædikterede værdi af PE. Sammenlign med skøn over spredningen i din slutmodel ovenfor.
8. Betragt den første måling i datafilen PowerOutputTest.csv. Lav både et 95%-konfidensinterval for middelværdien af PE og et 95%-prædiktionsinterval.

Luk igen

På et raffinaderi går råolien igennem flere processer, og undervejs skal man vurdere, hvad det nuværende produkt bedst kan anvendes til. Dette kræver mange tidskrævende analyser, og man er derfor interesseret i simple alternative analysemetoder. I skal i denne opgave se på anvendelsen af NIR-spektre til bestemmelse af cetantallet for dieselolie. Data er fra artiklen IPA: A deep CNN based on Inception for Petroleum Analysis. Filen TrainY124.txt indeholder data for 124 prøver med cetantal varierende mellem 19.0 og 71.1. Filen har 124 rækker og 3891 søjler, hvor den sidste søjle er cetantallet bestemt ved en kompliceret analysemetode. I skal i opgaven etablere en model ved henholdsvis backward selektion og ved ridge regression. I artiklen omtales partial least squares, som giver en model, der er af samme type (lineær), som de to modeller I skal betragte. Der betragtes også to ikke-lineære og noget mere komplicerede modeller i artiklen. Hvert NIR-spekter er målt ved 3890 bølgelængder (fra 4500 til 12000 $\text{cm}^{-1}$).

1. Indlæs data som en $124\times 3891$ matriks Traen med brug af indlæsningskommandoen np.loadtxt.
   Placer NIR-spektrene i en $124\times 3890$ matrix $T$ og cetantallet (søjle 3891) i vektoren cetan (python: T=Traen[:,0:3890]).
   Angiv, hvor mange olieprøver der har et cetantal under 27.
2. Kør forward selektionsalgoritmen, hvor der inkluderes op til 20 led. Illustrer proceduren med en figur som i afsnit 9.4 baseret på output fra forward.
   Et kald til cvForward giver følgende crossvalidation prædiktionsspredninger (tager meget lang tid at køre).
   $$\begin{array}{lccccccc} \hline \text{Antal variable} & 1 & 2 & 3 & 4 & 5 & 6 & 7 \\ \hline \text{CV-prædiktionsspredning} & 6.10 & 2.46 & 2.25 & 2.24 & 2.14 & 2.08 & 2.05 \\ \hline &&&&&&& \\ \text{Antal variable} & 8 & 9 & 10 & 11 & 12 & 13 & 14 \\ \hline \text{CV-prædiktionsspredning} & 2.07 & 2.04 & 2.09 & 2.07 & 2.12 & 2.11 & 2.14 \\ \hline \end{array}$$ Hvor mange variable, fundet ved forward selektion, vil du medtage i en model til beskrivelse af cetantallet ?
   Synes du at cetantallet er velbestemt ud fra NIR-spektret.
3. Du skal nu bruge ridge regression til at beskrive cetantallet ud fra NIR-spektrene. Med alle 3890 forklarende variable inkluderet tager det lang tid at køre cvRidge. Jeg har derfor lavet et mindre datasæt hvor der kun medtages 350 bølgelængder, nemlig søjlerne 501-750 og 3751-3850 blandt de 3890 søjler. Det reducerede datasæt med spektrene er i filen Train124subset.txt. Indlæs filen og brug denne delmængde af spektret til at lave ridge regression.
   Lav en tabel med krydsvalideringsspredning $s_{\text{cv}}(\lambda)$ for $\lambda$-værdierne $0.001,0.01,0.08,0.1,0.2,1,10,100,200,400,1000$.
   Lav også en figur, hvor $s_{\text{cv}}(\lambda)$ afsættes mod $\log(\lambda)$.
   Vælg en værdi $\hat\lambda$ af $\lambda$, som du vil bruge til at beskrive data med. Lav en figur med det observerede cetantal afsat mod det forventede cetantal, og indsæt identitetslinjen i denne figur. Du kan eventuelt lade dig inspirere af koden i den skjulte kode Ridge regression i python.
4. Hvilken model foretrækker du, modellen fra forward selektion eller modellen fra ridge regression ?

Kommentar omkring data. I artiklen omtales 249 prøver, der deles op i et træningssæt med 174 prøver og et testsæt med 75 prøver. Imidlertid ser det ud til, at 70 ud af de 75 i testsættet er indeholdt i de 174 prøver i træningssættet. Jeg har derfor set bort fra testsættet i artiklen, og bruger kun træningsættet fra artiklen i denne og den følgende opgave. De 174 prøver i artiklens træningssæt har jeg delt op på et nyt træningsæt med 124 prøver og et nyt testsæt med 50 prøver. Jeg har gjort dette på en sådan måde, at fordelingen af cetantallet i træningssættet og testsættet ligner hinanden. Sammenligning med ridge regression baseret på alle 3890 variable i spektrene: Følgende tabel viser cross-validation prædiktionsspredning $s_{\text{cv}}(\lambda)$ for nogle få værdier af regulariseringsparameteren $\lambda.$ $$\begin{array}{lcccccc} \hline \lambda & 0.1 & 0.2 & 0.9 & 1 & 1.2 & 2 \\ \hline \text{CV-prædiktionsspredning} & 2.105 & 2.041 & 1.976 & 1.975 & 1.976 & 1.981 \\ \hline \end{array}$$

Luk igen

Denne opgave er en fortsættelse af den foregående opgave. Udover datasættet med 124 olieprøver er der også et andet datasæt med 50 olieprøver. Vi kan bruge dette datasæt til et uafhængigt tjek af, hvor god modellen er til at prædiktere. Det nye datasæt er i filen TestY50.txt, som har samme struktur som TrainY124.txt (men altså 50 rækker i stedet for 124 rækker).

1. Indlæs som i foregående opgave data fra filen TrainY124.txt og dan en matriks $T$ med spektrene og en vektor cetan med cetantallene for de 124 prøver.
   Indlæs dernæst data fra filen TestY50.txt og dan en matrik Ttest med spektrene og en vektor cetantest med cetantallene for de 50 nye prøver.
2. Betragt først forward selektionsmodellen fra foregående opgave, hvor der medtages 10 variable. Estimer modellen baseret på $T$ og cetan, og lav prædikterede værdier for prøverne med spektre i Ttest. Dette kan gennemføres under brug af koden i det skjulte punkt Bruge slutmodellen til prædiktion.
   Udregn testspredningen som kvadratroden af den gennemsnitlige værdi af den kvadrerede afstand mellem prædiktion og den sande værdi, det vil sige
   $$s_{\text{test}}=\sqrt{\frac{1}{50} \sum_i (\text{cetantest}_i-\text{prædiktion}_i)^2}.$$ Overrasker værdien dig ?
3. Betragt nu ridge regression, hvor regulariseringsparameteren $\lambda$ vælges til $\lambda=0.1.$ Benyt funktionen ridge til at finde $\hat\alpha(\lambda),$ $\hat\beta(\lambda)$ og spredningsskøn $s_M(0.1)$ baseret på spektrene i $T$ og respons i cetan.
   Angiv den første koordinat i vektoren $\hat\beta(\lambda).$
4. Udregn de prædikterede værdier for de nye data med spektre i Ttest (se koden i det skjulte punkt Ridge prædiktion i python)
   Lav en figur hvor de observerede værdier for de nye data afsættes mod de prædikterede værdier. Indsæt identitetslinjen i figuren.
   Hvad er den største abolutte prædiktionsfejl, og hvor optræder den ?
5. Udregn testspredningen som kvadratroden af den gennemsnitlige værdi af den kvadrerede afstand mellem prædiktion og den sande værdi.
   Sammenlign denne værdi med cross-validation prædiktionsspredningen $s_{\text{cv}}(0.1)=1.98.$

Som nævnt i starten af foregående opgave betragtes i artiklen IPA: A deep CNN based on Inception for Petroleum Analysis partial least squares og to mere avancerede metoder. Desværre kan vi ikke bruge de testspredniger der nævnes i artiklen, men vi kan lave en ny beregning for partial least squares baseret på de data I har betragtet i denne og den foregående opgave. I partial least squares skal man vælge et antal "komponenter", og vælges 8 får man en cross-validation prædiktionsspredning på 2.14. Benyttes 8 komponenter bliver testspredningen fra de 50 prøver i TestY50.txt 1.54. Overordnet set, giver de tre metoder forward selektion, ridge regression og partial least squares lige gode resultater.

Luk igen

Til bestemmelse af mængden af protein i en prøve bruges ofte et Bradford assay. Lys med en bølgelængde på 595 nm sendes gennem prøven og absorbansen måles. Ud fra en kalibreringskurve laves der et skøn over proteinmængden. Imidlertid er respons som funktion af proteinmængden kun lineær i et mindre område, hvilket øger usikkerheden i metoden. I artiklen Linearization of the Bradford Protein Assay fra 2010 (baseret på artikel fra 1996) foreslås at bruge forholdet $a_{590}/a_{450}$ mellem absorbans ved bølgelængderne 590 nm og 450 nm, hvor teori og praksis peger på en lineær sammenhæng i et større område. I skal i denne opgave se på data fra artiklen, hvor det er muligt at lave et test for linearitet. For hver af 6 proteinmængder er der lavet 3 målinger af absorbans ved både 590 og 450 nm, svarende til brønd A, B og C i tabellen nedenfor. Data er i i filen TreGentagBradford.csv, der har tre søjler, hvor første søjle indeholder proteinmængden BSA (Bovine serum albumin, $\mu\text{g}$), den anden søjle er en variabel gruppe, der koder for BSA-gruppe (med værdierne G1,G2,G3,G4,G5,G6), og den tredje søjle indeholder absorbansratio. Data er gengivet i den følgende tabel, hvor søjlen F indeholder forholdet mellem måling ved 590 nm og 450 nm. $$\begin{array}{c|ccc|ccc|ccc}\hline && \text{590nm} &&& \text{450nm} &&& \text{F} & \\ \text{BSA} & \text{BrønA} & \text{BrønB} & \text{BrønC} & \text{BrønA} & \text{BrønB} & \text{BrønC} & \text{BrønA} & \text{BrønB} & \text{BrønC} \\ \hline 0 & 0.290 & 0.298 & 0.306 & 0.629 & 0.654 & 0.672 & 0.461 & 0.456 & 0.455 \\ 1 & 0.491 & 0.457 & 0.462 & 0.600 & 0.593 & 0.590 & 0.818 & 0.771 & 0.783 \\ 2 & 0.580 & 0.591 & 0.601 & 0.516 & 0.534 & 0.539 & 1.124 & 1.107 & 1.115 \\ 3 & 0.683 & 0.717 & 0.751 & 0.483 & 0.495 & 0.495 & 1.414 & 1.448 & 1.517 \\ 4 & 0.802 & 0.853 & 0.840 & 0.446 & 0.459 & 0.444 & 1.798 & 1.858 & 1.892 \\ 5 & 1.226 & 0.965 & 0.971 & 0.677 & 0.394 & 0.403 & 1.811 & 2.449 & 2.409 \\ \hline \end{array}$$ Et relevant spørgsmål for data af denne type er, om man skal tage gennemsnit for hver bølgelændge af de tre målinger, før man beregner forhold, eller om man skal lave forholdet for hver brønd og tage gennemsnit af disse. Formodentligt giver det en meget lille forskel, men da respons er lineær for forholdet, vil det være bedst at betragte forholdet for hver brønd og så tage gennemsnit. Et andet spørgsmål, I skal overveje nedenfor, er, om man skal basere kalibreringskurven på gennemsnit af de tre forhold for hver proteinmængde, eller om man skal basere beregningerne på alle 18 datapunkter.

1. Opstil normalfordelingsmodellen, hvor hver BSA-gruppe har sin egen middelværdi af absorbansratio, og alle grupperne har den samme varians.
   Analyser modellen i python, og angiv skøn over spredningen i modellen.
2. Opstil regressionsmodellen, hvor middelværdien af absorbansratio afhænger lineært af BSA.
   Analyser modellen i python, og angiv skøn over spredningen i modellen.
   Hvis man kun betragter data bestående af 6 proteinmængder og 6 absorbansratioer, som hver er gennemsnit af tre målinger, vil en regressionsmodel give et skøn over spredningen på 0.02665. Kan du forstå dette i sammenligning med det skøn, du fik ovenfor (svaret kan godt være "nej") ?
3. Forklar, at regressionsmodellen er en undermodel af modellen, hvor hver BSA-gruppe har sin egen middelværdi af absorbansratio (den ensidede variansanalysemodel).
   Lav et test for reduktion fra den ensidede variansanalysemodel til regressionsmodellen. Angiv den $F$-fordeling, der bruges. Er $p$-værdien stor nok, til at du vil acceptere regressionsmodellen ?
4. Beregn et 95%-konfidensinterval for hældningen i den lineære sammenhæng mellem middelværdien af Absorbansratio og BSA.
   Hvis man igen betragter regressionsmodellen baseret på de 6 gennemsnit, bliver konfidensintervallet $[0.335,0.371].$ Sammenlign dette med det konfidensinterval du lige har fundet.
5. Betragt en ny prøve med ukendt indhold af proteinmængden BSA, hvor absorbansratio er målt tre gange med værdierne 0.979, 0.977, 0.984. Angiv et 95%-konfidensinterval for det ukendte indhold af BSA.
   Hvis man igen betragter regressionsmodellen baseret på de 6 gennemsnit, og en ny måling som er 0.980 (gennemsnittet af de tre målinger ovenfor), bliver konfidensintervallet $[1.31,1.78].$
   Ved fremtidig brug af kalibreringskurven til bestemmelse af proteinindholdet vil du da vælge at bruge kurven og spredningskøn baseret på alle 18 datapunkter, eller vil du bruge kurven og spredningsskøn baseret på de 6 gennemsnit ?

Luk igen

I denne sidste opgave er det meningen, at I hovedsageligt skal bruge metoderne fra de 6 første øvelser (undtagelsen er det sidste spørgsmål). I artiklen Use of Multiple Linear Regression Models for Setting Water Quality Criteria for Copper: A Complementary Approach to the Biotic Ligand Model betragtes forskellige modeller til at beskrive giftigheden af kobber i vandløb. Giftigheden, under givne forhold, findes ved at lave en række forsøg med forskellige koncentrationer af kobber i det vand hvor vandlopper (Daphnia magna) opholder sig i et givet tidsrum. Som beskrivelse af giftigheden bruges den koncentration af kobber under hvilken 50 procent af vandlopperne dør. Denne kaldes acute EC50 (Half maximal effective concentration). Som modeller for EC50-koncentrationen har man på den ene side en simplificeret model, hvor man kun bruger vandets hårdhed som forklarende variabel i en regressionsmodel. På den anden side har man en kompliceret multipel regressionsmodel med 10 forklarende variable, som ifølge forfatterne ikke har vundet indpas i vandkvalitetsprogrammer. Forfatterne betragter i stedet en mere simpel multipel regressionsmodel med tre forklarende variable, henholdsvis DOC (dissolved organic carbon, mg/L), hårdhed (mg/L) og pH. I opgaven her skal I betragte data for EC50-koncentrationen af kobber ($\mu$g/L Cu). Den bedste skala at betragte data på er ved at logaritmetransformere værdierne (dog ikke pH, som allerede er på en log-skala). Filen Dmagna.csv indeholder variablene DOC, hardness, pH og EC50. Der er data for 302 målinger. Lad logEC50 være logaritmen til acute EC50, og for nemhed i notationen lad $D$ være logaritmen til DOC og $H$ være logaritmen til hardness. En analyse af den multiple regression af logEC50 på $D$, $H$ og pH viser, at modellen kan reduceres til $$\mathrm{LogEC50}\sim N(\alpha+D+\frac{1}{2}H+\mathit{pH},\sigma^2).$$ Du kan eventuelt køre en multipel regressionsanalyse og se, at parametertabellen ikke strider mod $\beta_D=1,$ $\beta_H=\frac{1}{2}$ og $\beta_{\text{pH}}=1.$ Denne analyse giver anledning til at indføre et indeks, her kaldet giftIndex, på formen $$\text{giftIndex}=\text{logEC50}-D- \frac{1}{2}H-\mathit{pH}.$$ Fra data ses, at dette indeks svinger omkring værdien $-6.86.$ Opgaven her går ud på at lave nogle (lidt usædvanlige) undersøgelser for at se, om den multiple regressionsmodel kan forbedres ved at inkludere flere variable. Specifikt skal I se, om produktet $H\cdot\mathit{pH}$ indeholder information om giftindex ved at dele data op i grupper efter denne variabel.

1. Indlæs data og dan variablene logEC50, $D$, $H$ og pH. Dan dernæst giftIndex.
   Tjek at gennemsnit af målingerne af giftIndeks er $-6.86.$
2. Konstruer en kategorisk variabel gruppe, der deler op i tre grupper, $A$, $B$ og $C$, efter om $H\cdot\mathit{pH}$ er mindre end 25, er mellem 25 og 35 eller over 35. I python (pandas) kan dette gøres med kommandoen
   $$\text{gruppe=pd.cut(H*pH,bins=[0,25,35,100],labels=['A','B','C'])}$$
   Dan nu to datasæt, giftA og giftC, med værdierne af giftIndex hørende til gruppe $A$ og til gruppe $C$. Der er 60 observationer i gruppe $A$ og 96 i gruppe $C$.
3. For de 60 prøver i gruppe $A$ er der 14, der har en værdi af giftIndex under $-7$ (konstater, at dette er rigtigt). Opstil en statistisk model til beskrivelse af observationen 14, og lav et 95%-konfidensinterval for sandsynligheden for, at værdien af giftIndex er under $-7$ i gruppe $A$.
4. For de 96 prøver i gruppe $C$ er der 29, der har en værdi af giftIndex under $-7.$ Undersøg, om der er samme frekvens af prøver med en værdi af giftIndex under $-7$ i de to grupper $A$ og $C$.
5. Opstil en statistisk model til beskrivelse af værdierne af giftIndex i giftA og giftC.
   Lav et test for hypotesen, at der er samme middelværdi af giftIndex i de to grupper $A$ og $C$.
6. Betragt nu giftIndex for alle 3 grupper dannet ud fra faktoren gruppe. Opstil en statistisk model for data, og undersøg først, om der er samme varians for de tre grupper, og dernæst, om der er samme middelværdi for de tre grupper.
7. Ovenstående analyse viser, at der er information i produktet $H\cdot\mathit{pH}$, der kan bruges til at beskrive giftIndex, og dermed kan bruges til at beskrive logEC50. I dette delspørgsmål skal I analysere en multipel regressionsmodel til beskrivelse af logEC50, hvor I som forklarende variable ud over $D$, $H$ og pH også bruger alle produkter af disse, det vil sige, variablene
   $$\begin{aligned} & \mathit{DD}=D*D,\enspace \mathit{HH}=H*H,\enspace \mathit{ppH}=\mathit{pH}*\mathit{pH}, \\ & \mathit{DH}=D*H,\enspace \mathit{DpH}=D*\mathit{pH},\enspace \mathit{HpH}=H*\mathit{pH}. \end{aligned}$$ Opskriv den multiple regressionsmodel, hvor alle 9 forklarende variable inddrages.
   Reducer modellen ved brug af backward selektion.
8. Lav et test for reduktion fra startmodel til slutmodellen fremkommet ved backward selektion.
   Angiv et 95%-konfidensinterval for regressionskoefficienten hørende til variablen $H\cdot\mathit{pH}$ i slutmodellen.
   I din slutmodel, hvor meget stiger middelværdien af logEC50 hvis dissolved organic carbon stiger fra 2 til 10 ?
9. Overvej, hvilke figurer du ville lave til kontrol af slutmodellen (du skal ikke lave figurerne).
10. Sammenlign endelig skøn over spredning i slutmodellen med skøn over spredning i den multiple regressionsmodel hvor kun de tre forklarende variable $D$, $H$ og pH inddrages. Er det rimeligt at bruge denne sidste simple model ?

Luk igen

Betragt data fra opgave 7.1 og slutmodellen ved backward selektion du fandt frem til i opgaven. Analyse af slutmodellen giver $$\text{Cov}_s(\hat\alpha,\hat\beta_{\text{AT}})=-5.981,\enspace \text{Cov}_s(\hat\alpha,\hat\beta_{\text{AT2}})=0.6301,\enspace \text{Cov}_s(\hat\beta_{\text{AT}},\hat\beta_{\text{AT2}})=-0.001295.$$

1. Lav et approksimativt 95%-konfidensinterval for parameteren $\theta=\alpha/1000+30\big(\beta_{\text{AT}}^2+\beta_{\text{AT2}}^2\big).$

Luk igen