0

I have some data I'm trying plot a fit for in ggplot. I'm using nls with a self starting four parameter logisitic regression (SSfpl). It can fit the formula fine using nls and it gives me coefficients but when i try to pass it to geom_smooth it gives me a computation failed in 'stat_smooth()': singular matrix 'a' in solve.

My data looks like so:

dput(SODat)
structure(list(Treatment = structure(1:6, .Label = c(".0032uM", 
".016uM", ".08uM", ".4uM", "2uM", "10uM"), class = "factor"), 
    Peak = c(1028.06247583425, 1466.14079105711, 1007.51610918954, 
    1171.03478670369, 1204.37248245019, 1042.44546803658), Area1 = c(902.715126326407, 
    818.917775151066, 1003.80586032706, 1268.82827952821, 1479.87729958222, 
    1321.77972824089), SS1 = c(131.573938723728, 90.0065951864569, 
    111.723312483588, 147.027334638628, 195.159502773361, 173.819528362056
    ), Conc_nm = c(3.2, 16, 80, 400, 2000, 10000), logconc_nm = c(0.505149978319906, 
    1.20411998265592, 1.90308998699194, 2.60205999132796, 3.30102999566398, 
    4)), .Names = c("Treatment", "Peak", "Area1", "SS1", "Conc_nm", 
"logconc_nm"), class = "data.frame", row.names = 2:7)

My code for it looks as follows. 

    val=NULL # EC50
vallog=NULL# logEC50
allDR=NULL
vertline=NULL
for (i in 2:4){
  currentfit=tryCatch(nls(SODat[,i] ~ SSfpl(logconc_nm,A,B,xmid,scal),dat=testdat),error=function(e) 0)
  if(typeof (currentfit)=='list'){
    vallog[i]= summary(currentfit)$coefficients[3]
    val[i]=10^summary(currentfit)$coefficients[3]}
  else{
    vallog[i]=0
    val[i]=0
  }
}  

vallog=vallog[-c(1:2)]
val=val[-c(1:2)]

My plotting code looks likes this:

    ## SS1
DRPlot_SS=ggplot(data = subset(testdatMelt,testdatMelt$variable == 'SS1'),aes(logconc_nm,value))
DRPlot_SS = DRPlot_SS +  geom_point()
DRPlot_SS = DRPlot_SS +  scale_x_log10(breaks=round(testdat$logconc_nm,2))

if((vertline$vallog[2])!= 0){
  DRPlot_SS = DRPlot_SS +  geom_smooth(method = 'nls', formula = y ~ SSfpl(x,A,B,xmid,scal),se=FALSE)
  DRPlot_SS = DRPlot_SS +  geom_vline(data=vertline,color='red',aes(xintercept = vallog[2]),alpha=.5)
  DRPlot_SS = DRPlot_SS + geom_text(data=vertline,aes(x=(vallog[2])-.5,y=max(testdatMelt$value[testdatMelt$variable=='SS1']),color='red',
                                                          label = paste('EC50 =',round(val[2],2),'nM')))}
if(vertline$vallog[2] == 0){DRPlot_area=DRPlot_area + geom_text(label='No Fit',size=10,col='red',aes(x=median(testdatMelt$logconc_nm),
                                                                                     y=max(testdatMelt$value[testdatMelt$variable=='SS1'])))}

DRPlot_SS = DRPlot_SS +  theme(legend.position='none')
DRPlot_SS = DRPlot_SSylim(0,max(testdatMelt$value[testdatMelt$variable=='SS1'])+100)

****Question**** Does anybody know why it would be able to fit in nls and not plot in ggplot? In this example the fit is probably pretty weak. So that's really my only explanation.

Ted Mosby
  • 1,426
  • 1
  • 16
  • 41
  • 1
    in general, when you use `$` notation within `aes()` mappings, you're asking for trouble, especially when you have multiple data sets or are subsetting across layers. The *best* solution would be to invest more time in learning how to use ggplot2 differently, but if time doesn't allow that, I'd first suggest adding a `inherit.aes = FALSE` to the smooth, vline and text layers. That will result in those layers not looking for the mappings that you established in the original `ggplot()` call. – Matt74 Feb 15 '16 at 18:30
  • Hmm. interesting. I will give it a shot! thanks! – Ted Mosby Feb 15 '16 at 18:34
  • So the error is gone with the inherit.aes = FALSE, but the smooth still isn't plotting. any further thoughts? – Ted Mosby Feb 15 '16 at 18:42
  • afraid not, sorry. But you might get more useful help if you gave a reproducible example. We don't have access to `testdatMelt`. – Matt74 Feb 15 '16 at 19:39
  • shoot you're right. I'll get around to it. For now, it only seems to be a problem on a small subset of data. my guess is that it's just too weak of a fit to plot. – Ted Mosby Feb 15 '16 at 20:00
  • 2
    If the model runs outside of ggplot, why dont you get the fitted values from that and plot them using geom_line? – user2957945 Feb 15 '16 at 22:53
  • Would geom_line fit a four parameter logistic regression? I thought it wouldn't which is why I pass the formula through geom_smooth. – Ted Mosby Feb 16 '16 at 13:48
  • If you mak predictions on small intervals of x , it should look smooth enough – user2957945 Feb 18 '16 at 10:06

0 Answers0