Doppelter Logarithmus und Ermittlung der Regressionsgeraden(gleichung)

[All-Right]

Hallo zusammen,

für eine Auswertung bezüglich der Ermittlung des sogenannten Mass Median Aerodynamic Diameter (MMAD; die mittlere Teilchengröße für die Angabe von Pulvern zur Inhalation) werden die Pulver-Teilchengrößen aus sieben verschiedenen Fraktionen (Schalen) aus einem Next Generation Impactor (NGI; wird in der pharmazeutischen Technologie verwendet) gesammelt und deren Grenzdurchmesser gegen den prozentualen (Wirkstoff)Anteil aufgetragen. Im studentischen Praktikum geschieht die graphische Auftragung auf doppelt-logarithmischem Papier und es wird anschließend eine Regressionsgerade eingezeichnet, wobei dann der Schnittpunkt bei 50% der Regressionsgeraden mit dem Schnittpunkt auf der x-Achseder der zu ermittelnden Größe des MMAD entspricht. Dies finde ich ziemlich ungenau, da jeder die Gerade zwischen den einzelnen Punkten etwas anders einzeichnen würde. Von dem her muss es ja auch eine genauere mathematische Herangehensweise geben.

Im Anhang habe ich eine CSV-formatierte Excel-Tabelle mit den Werten angehängt, die wir aus einem Versuch erhalten haben. Anschließend habe ich mit R die einzelnen Grenzdurchmesser [µm] der Schalen gegen den Wirkstoffanteil [%] aufgetragen und schon einmal eine Horizontale bei 50% einzeichnen lassen:

Code: Alles auswählen

NGI <- read.csv2(file.choose()) 
NGI 	# Anmerkung: d0 ist keine Schale, nur Aufsummierung aller Fraktionen
plot(NGI$Grenzdurchmesser, NGI$Wirkstoffanteil, log="xy", xlim=c(0.1,10), ylim=c(1,100), 
	xlab="Grenzdurchmesser [µm]", ylab="Wirkstoffanteil [%]")
abline(h=50)

Da eine Kollegin von mir das ganze in Excel ausgerechnet hat und eine Potenz-Regressionsgerade/-kurve (power law regression) zur Ermittlung des MMAD verwendete, habe ich analog dasselbe über das "händische" Rückrechnen von x aus y über die Ausgabe der Koeffizienten (Intercept) versucht:

Code: Alles auswählen

 summary(lm(log(NGI$Wirkstoffanteil)~log(NGI$Grenzdurchmesser)))
# Intercept estimate: 		2.63466
# log(NGI$Grenzdurchmesser):	0.94268

# ln(y) = 2.63466 + 0.94268*ln(x)
# y = e^2.63466 + e^0.94268*ln(x)
# y = 13.93857 * x^(0.94268)		--> x = (y/13.93857)^(1/0.94268)
(50/13.93857)^(1/0.94268)		# --> 3.876892 µm MMAD bei 50%
[attachment=0]NGI_Attempt5.csv[/attachment]
abline(v=(50/13.93857)^(1/0.94268))

Als nächstes wollte ich anhand dieses Ergebnisses meine Regressionsgerade einzeichnen. Die Steigung 0.94268 passt perfekt, jedoch ist der Ordinatenabschnitt (Intercept) verschoben:

Code: Alles auswählen

abline(log(2.63466),0.94268)

Ich weiß nicht, wie ich diesen genau ermitteln kann. Er scheint um log(2.64466) + x verschoben zu sein. Kann mir jemand weiterhelfen?
Vielen Dank schon mal für eure Hilfe!

P.S. Falls dieser Beitrag diesem Unterforum falsch zugeordnet sein sollte, bitte ich um Entschuldigung.

[bigben]

Hallo All-Right,

willkommen im Forum. Ich muss sagen, dass mich die Vermischung aus pharmazeutischem und rechnerischem Inhalt etwas durcheinander bringt und auch diese Beschreibung kann man bestimmt nochmal leichter darstellen:

es wird anschließend eine Regressionsgerade eingezeichnet, wobei dann der Schnittpunkt bei 50% der Regressionsgeraden mit dem Schnittpunkt auf der x-Achseder der zu ermittelnden Größe des MMAD entspricht.

Beschreib nochmal in Ruhe, was gemeint ist, denn 50% einer Regressionsgeraden gibt es nicht.

Aufgefallen ist mir aber dieser Rechenschritt:

Code: Alles auswählen

# ln(y) = 2.63466 + 0.94268*ln(x)
# y = e^2.63466 + e^0.94268*ln(x)

Das scheint mir falsch zu sein.

Code: Alles auswählen

ln(y) = 2.63466 + 0.94268*ln(x)   | auf beiden Seiten die Exponentialfunktion
y = e^(2.63466 + 0.94268*ln(x))   | Wenn wir jetzt das "Plus" im Exponent auflösen entsteht ein Produkt, keine Summe
y = e^(2.63466) * e^(0.94268*ln(x)) 

Um auf Deine konkreten Daten zurück zu kommen:
Am einfachsten wird es sein, beide Spalten logarithmisch zu transformieren und dann damit einfach linear zu rechnen:

Code: Alles auswählen

NGI <- read.csv2("http://forum.r-statistik.de/download/file.php?id=1830")

NGI$GD.log <- log(NGI$Grenzdurchmesser)
NGI$WA.log <- log(NGI$Wirkstoffanteil)

regr <- lm(WA.log ~ GD.log, data = NGI)
coefficients(regr)

So, jetzt suchen wir die Stelle, wo diese Regression den Wert log(50) ergibt:

Code: Alles auswählen

3.912023 = 2.6346619 + 0.9426844*GD.log  
0.9426844*GD.log = 3.912023 - 2.6346619
GD.log = (3.912023 - 2.6346619)/0.9426844
GD.log = 1.355025
GD = e^1.355025
GD = 3.876858

Und das passt als Schnittpunkt auch ganz gut zu Deinen Daten:

Code: Alles auswählen

plot(NGI$Grenzdurchmesser, NGI$Wirkstoffanteil, log="xy", xlim=c(0.1,10), ylim=c(1,100), 
     xlab="Grenzdurchmesser [µm]", ylab="Wirkstoffanteil [%]")

abline(h = 50)
abline(v = 3.88)

Klärt das das Problem?

Im base-Grafiksystem würde ich auf die automatischen logarithmischen Achsen in der Regel verzichten und so wie oben vorgemacht, die Logarithmierung händisch vornehmen. Das ist im Zweifel immer einfacher. Hier mal im ganzen Stück für Deine konkreten Daten:

Code: Alles auswählen

NGI <- read.csv2("http://forum.r-statistik.de/download/file.php?id=1830")

NGI$GD.log <- log(NGI$Grenzdurchmesser)
NGI$WA.log <- log(NGI$Wirkstoffanteil)

regr <- lm(WA.log ~ GD.log, data = NGI)
coefficients(regr)

# Koordinatensystem (Beschriftung, Grenzen etc. kannst Du noch anpassen)
plot(WA.log ~ GD.log, data = NGI, xlab  = "Grenzdur...", ylab = "Wirksst...",
     xaxt = "n", yaxt = "n")
     
# Achseneinteilung
abline(h = log(seq(10, 100, 10)), col = "lightgrey")
abline(v = log(seq(.5, 10, .5)), col = "lightgrey")
axis(1, at = log(c(.25, .5, 1, 2, 4, 8)), labels = c(.25, .5, 1, 2, 4, 8))
axis(2, at = log(c(2.5, 5, 10, 20, 40, 80)), labels = c(2.5, 5, 10, 20, 40, 80))
axis(3, at = 1.355, labels = "3.877")
axis(4, at = log(50), labels = "50 %")

# Geraden einzeichnen
abline(regr, col = "darkgrey", lwd = 2)
abline(h = log(50), lty = 2, lwd = 1.5)
abline(v = 1.355025, lty = 2, lwd = 1.5)

# Punkte zum Schluss nochmal deutlich nachzeichnen
points(WA.log ~ GD.log, data = NGI, pch = 15)

Rplot01.pdf

(5.06 KiB) 12-mal heruntergeladen

LG,
Bernhard

[All-Right]

Hallo Bernhard,
vielen lieben Dank für Deine schnelle Antwort. Ja, das hilft ungemein und ist letztlich meine gesuchte Antwort
Es tut mir leid, dass ich mich etwas ungünstig ausgedrückt habe. Ich meinte natürlich nicht 50% einer Regressionsgeraden, ich meinte, dass der Schnittpunkt des 50% Anteils der Teilchen mit der Regressionsgerade den MMAD auf der x-Achse markieren
Und auch Danke für die Korrektur der Rechnung

Das ist eine wirklich schöne Herangehensweise, erst alles zu logarithmieren und dann nach Deinem Vorschlag linear weiterzurechnen. Auch Dein Vorschlag zur Darstellung mit den geänderten Achsen finde ich sehr schön!

Herzlichen Dank nochmals und hab noch einen schönen Tag und eine schöne Woche!
Flo