Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

R - ggplot2 extrapolated regression lines in linear region

Tags:

r

ggplot2

I have the following data frame (attached). I have plotted CR vs A for various Collimator / Head combinations.

(p <- ggplot(df,aes(x=A,y=CR,col=Head))+geom_point()+geom_line() +facet_grid(Collimator~Head, scales="fixed")+scale_x_continuous("Activity [MBq]", expand = c(0,0))+ylim(0,80000)+ ylab("Count Rate [cps]") + theme_bw()+theme(legend.position = "none"))

enter image description here

In an ideal world the above plot would be linear. In reality CR will begin to fall off with A due to increased system deadtime. What I would like to add to each facet is a straight line fit which passes through the first 2 data points only - this is before dead time kicks in.

Is there a simple way to do this in ggplot2. Can I use geom_smooth(method = "lm") with some additional options?

structure(list(A0 = c(76L, 274L, 786L, 1060L, 1294L, 2092L, 2639L, 
3437L, 4223L, 76L, 274L, 786L, 1060L, 1294L, 2092L, 2639L, 3437L, 
4223L, 76L, 274L, 786L, 1060L, 1294L, 2092L, 2639L, 3437L, 4223L, 
76L, 274L, 786L, 1060L, 1294L, 2092L, 2639L, 3437L, 4223L, 76L, 
274L, 786L, 1060L, 1294L, 2092L, 2639L, 3437L, 4223L, 76L, 274L, 
786L, 1060L, 1294L, 2092L, 2639L, 3437L, 4223L), T = c(85L, 87L, 
88L, 89L, 89L, 90L, 91L, 92L, 93L, 97L, 98L, 99L, 100L, 101L, 
102L, 103L, 103L, 104L, 306L, 308L, 310L, 311L, 313L, 315L, 316L, 
317L, 321L, 328L, 330L, 331L, 332L, 336L, 338L, 340L, 341L, 342L, 
352L, 354L, 357L, 358L, 361L, 363L, 364L, 366L, 368L, 376L, 378L, 
379L, 380L, 385L, 386L, 388L, 389L, 390L), A = c(64.8860944957, 
233.0628375794, 667.3247389509, 898.2821229937, 1096.5821388243, 
1769.5416286837, 2228.0796200189, 2896.4301555482, 3552.1951906822, 
63.4538798403, 228.3428223318, 653.8100019998, 880.0900106808, 
1072.3775535078, 1730.4829756141, 2178.8997717016, 2837.7713207042, 
3480.2557273313, 43.0160516527, 154.5083447781, 441.5789350636, 
594.4069507774, 722.9307781622, 1164.4170158664, 1466.1502053416, 
1905.9469999005, 2324.4554021136, 41.2913626085, 148.3134747414, 
424.6633733747, 571.6370071823, 692.6559941892, 1115.653739112, 
1402.1405840581, 1822.7365994889, 2235.4125025879, 39.4886520314, 
141.8383609947, 404.6153215597, 544.6504360348, 661.1848647194, 
1064.9635212237, 1340.925513839, 1739.9197611259, 2129.880289335, 
37.7646447282, 135.6459396877, 388.3926422248, 522.8131775262, 
632.318659446, 1020.3651426265, 1282.3829893435, 1667.0556688926, 
2044.484697239), Counts = c(102830L, 328231L, 784020L, 1010212L, 
1160531L, 1582051L, 1760850L, 1888034L, 1897347L, 99780L, 317952L, 
749548L, 965314L, 1106831L, 1488386L, 1672990L, 1793667L, 1789803L, 
129507L, 453800L, 1053106L, 1327867L, 1473197L, 1900706L, 2075742L, 
1991265L, 1756820L, 121230L, 424329L, 994864L, 1237568L, 1374478L, 
1734922L, 1921046L, 1878514L, 1664225L, 213389L, 712467L, 1498082L, 
1777791L, 1882367L, 1824631L, 1525162L, 1250229L, 1072038L, 193591L, 
651249L, 1354850L, 1594421L, 1653835L, 1669993L, 1436444L, 1144518L, 
1015859L), CR = c(3428L, 10941L, 26134L, 33674L, 38684L, 52735L, 
58695L, 62934L, 63245L, 3326L, 10598L, 24985L, 32177L, 36894L, 
49613L, 55766L, 59789L, 59660L, 4317L, 15127L, 35104L, 44262L, 
49107L, 63357L, 69191L, 66376L, 58561L, 4041L, 14144L, 33162L, 
41252L, 45816L, 57831L, 64035L, 62617L, 55474L, 7113L, 23749L, 
49936L, 59260L, 62746L, 60821L, 50839L, 41674L, 35735L, 6453L, 
21708L, 45162L, 53147L, 55128L, 55666L, 47881L, 38151L, 33862L
), Head = structure(c(2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L), .Label = c("H1", 
"H2"), class = "factor"), Collimator = structure(c(1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 
3L, 3L, 3L), .Label = c("HRGP", "LEGP", "LEHS"), class = "factor")), .Names = c("A0", 
"T", "A", "Counts", "CR", "Head", "Collimator"), row.names = c(NA, 
-54L), class = "data.frame")
like image 951
moadeep Avatar asked Dec 09 '22 10:12

moadeep


1 Answers

This should do it:

library(ggplot2)
(p <- ggplot(df,aes(x=A,y=CR,col=Head))+
    geom_point()+geom_line() +
    facet_grid(Collimator~Head, scales="fixed")+
    scale_x_continuous("Activity [MBq]",
                       expand = c(0,0))+ylim(0,80000)+
  ylab("Count Rate [cps]") + theme_bw()+theme(legend.position = "none"))

 library(plyr)
 subdf <- ddply(df,c("Collimator","Head"),
                function(x) x[1:2,])
 p + geom_smooth(method="lm",data=subdf,colour="gray",se=FALSE,
                 fullrange=TRUE)

enter image description here

like image 110
Ben Bolker Avatar answered Jan 07 '23 00:01

Ben Bolker