[R] extension to nlme self start SSmicmen?

[David Airey]

.

Just to close the query, here are the data and my R code to use nlme  
to test for a strain difference in receptor binding using nlme and the  
self start function SSmicmen.

Cheers,

Dave

"rownames","y","x","strain","pool"
"1",1,54.0423393249512,0.169009417295456,"dba","5_10"
"2",2,96.1945190429688,0.158872991800308,"dba","5_10"
"3",3,155.323593139648,0.669624865055084,"dba","5_10"
"4",4,153.497665405273,0.681416213512421,"dba","5_10"
"5",5,245.413208007812,1.43854534626007,"dba","5_10"
"6",6,219.874313354492,1.46481728553772,"dba","5_10"
"7",7,376.499328613281,2.15823173522949,"dba","5_10"
"8",8,275.310150146484,2.24408102035522,"dba","5_10"
"9",9,272.352142333984,3.06637287139893,"dba","5_10"
"10",10,258.903350830078,3.01941442489624,"dba","5_10"
"11",11,366.927612304688,4.61828088760376,"dba","5_10"
"12",12,312.222229003906,4.53119039535522,"dba","5_10"
"13",13,424.75,6.34457635879517,"dba","5_10"
"14",14,340.815551757812,6.15901756286621,"dba","5_10"
"15",15,369.316772460938,8.0274305343628,"dba","5_10"
"16",16,357.210906982422,8.27318668365479,"dba","5_10"
"17",17,83.4424362182617,0.319814652204514,"dba","1_6"
"18",18,85.219223022461,0.315470457077026,"dba","1_6"
"19",19,141.1943359375,0.686587870121002,"dba","1_6"
"20",20,138.046035766602,0.65224814414978,"dba","1_6"
"21",21,200.457916259766,1.35704016685486,"dba","1_6"
"22",22,195.697128295898,1.37710618972778,"dba","1_6"
"23",23,227.243957519531,2.02066588401794,"dba","1_6"
"24",24,228.378860473633,2.06659007072449,"dba","1_6"
"25",25,303.666809082031,2.84481954574585,"dba","1_6"
"26",26,233.546936035156,2.74655818939209,"dba","1_6"
"27",27,252.836318969727,4.49436807632446,"dba","1_6"
"28",28,238.709274291992,4.36507701873779,"dba","1_6"
"29",29,265.972564697266,5.97635507583618,"dba","1_6"
"30",30,267.110992431641,5.85864877700806,"dba","1_6"
"31",31,307.06640625,7.56363725662231,"dba","1_6"
"32",32,312.790069580078,7.44117259979248,"dba","1_6"
"33",33,72.9625701904297,0.303885966539383,"dba","12_8"
"34",34,76.9071884155273,0.324158817529678,"dba","12_8"
"35",35,123.434700012207,0.68782901763916,"dba","12_8"
"36",36,102.190269470215,0.663625717163086,"dba","12_8"
"37",37,175.760848999023,1.31111598014832,"dba","12_8"
"38",38,157.757446289062,1.28505086898804,"dba","12_8"
"39",39,241.185928344727,2.08934545516968,"dba","12_8"
"40",40,163.552764892578,1.98529183864594,"dba","12_8"
"41",41,204.231719970703,2.89943218231201,"dba","12_8"
"42",42,179.495147705078,2.78358721733093,"dba","12_8"
"43",43,206.535339355469,4.26909112930298,"dba","12_8"
"44",44,338.019897460938,4.11042499542236,"dba","12_8"
"45",45,261.55615234375,5.70411968231201,"dba","12_8"
"46",46,209.903076171875,5.56076145172119,"dba","12_8"
"47",47,203.267913818359,7.45896291732788,"dba","12_8"
"48",48,203.768753051758,7.4076600074768,"dba","12_8"
"49",49,50.6933288574219,0.151632696390152,"dba","7_2"
"50",50,57.050724029541,0.157838672399521,"dba","7_2"
"51",51,130.287200927734,0.641491115093231,"dba","7_2"
"52",52,147.611511230469,0.645835280418396,"dba","7_2"
"53",53,180.577438354492,1.36035001277924,"dba","7_2"
"54",54,181.326904296875,1.34773123264313,"dba","7_2"
"55",55,203.960754394531,2.02666497230530,"dba","7_2"
"56",56,212.268936157227,1.94702160358429,"dba","7_2"
"57",57,237.200775146484,2.75938391685486,"dba","7_2"
"58",58,224.804031372070,2.9393572807312,"dba","7_2"
"59",59,261.538757324219,4.02767848968506,"dba","7_2"
"60",60,246.376846313477,4.08829021453857,"dba","7_2"
"61",61,332.598663330078,5.55889987945557,"dba","7_2"
"62",62,283.046020507812,5.4238166809082,"dba","7_2"
"63",63,293.742919921875,7.18093538284302,"dba","7_2"
"64",64,261.284545898438,6.96331214904785,"dba","7_2"
"65",65,62.1605453491211,0.168802559375763,"dba","11_9"
"66",66,57.884651184082,0.157838672399521,"dba","11_9"
"67",67,189.111709594727,0.709549963474274,"dba","11_9"
"68",68,176.258209228516,0.733339548110962,"dba","11_9"
"69",69,230.409515380859,1.51632690429688,"dba","11_9"
"70",70,218.730087280273,1.45075035095215,"dba","11_9"
"71",71,285.564056396484,2.32579302787781,"dba","11_9"
"72",72,265.167175292969,2.39964413642883,"dba","11_9"
"73",73,289.183013916016,3.11540007591248,"dba","11_9"
"74",74,317.505859375,3.029137134552,"dba","11_9"
"75",75,438.351287841797,4.74446868896484,"dba","11_9"
"76",76,414.60400390625,4.58456182479858,"dba","11_9"
"77",77,416.385955810547,6.29844522476196,"dba","11_9"
"78",78,421.841491699219,6.11992025375366,"dba","11_9"
"79",79,469.797760009766,7.46703052520752,"dba","11_9"
"80",80,363.786895751953,7.76222848892212,"dba","11_9"
"81",81,39.6947517395020,0.133842229843140,"dba","3_4"
"82",82,40.8877220153809,0.151218950748444,"dba","3_4"
"83",83,119.681869506836,0.624528110027313,"dba","3_4"
"84",84,123.664993286133,0.648938238620758,"dba","3_4"
"85",85,164.695404052734,1.35062730312347,"dba","3_4"
"86",86,167.642990112305,1.39737904071808,"dba","3_4"
"87",87,204.152923583984,2.00204801559448,"dba","3_4"
"88",88,204.188491821289,2.10630846023560,"dba","3_4"
"89",89,225.391983032227,2.8406822681427,"dba","3_4"
"90",90,238.616790771484,2.90936160087585,"dba","3_4"
"91",91,258.135070800781,4.23992300033569,"dba","3_4"
"92",92,280.171173095703,4.25109386444092,"dba","3_4"
"93",93,346.150909423828,5.78045320510864,"dba","3_4"
"94",94,345.513641357422,5.96870136260986,"dba","3_4"
"95",95,420.243377685547,7.73554277420044,"dba","3_4"
"96",96,246.498519897461,7.43207025527954,"dba","3_4"
"97",97,47.404598236084,0.152460157871246,"c57","12_4"
"98",98,43.5518913269043,0.155563145875931,"c57","12_4"
"99",99,139.432891845703,0.7014821767807,"c57","12_4"
"100",100,135.885025024414,0.70313709974289,"c57","12_4"
"101",101,186.337890625,1.44599246978760,"c57","12_4"
"102",102,212.822143554688,1.40275752544403,"c57","12_4"
"103",103,234.293060302734,2.08686304092407,"c57","12_4"
"104",104,229.746109008789,2.00080680847168,"c57","12_4"
"105",105,238.368408203125,2.77076148986816,"c57","12_4"
"106",106,227.190399169922,2.87295341491699,"c57","12_4"
"107",107,255.688491821289,4.46706199645996,"c57","12_4"
"108",108,265.674682617188,4.45009899139404,"c57","12_4"
"109",109,297.111419677734,5.97180414199829,"c57","12_4"
"110",110,291.05224609375,5.79493379592896,"c57","12_4"
"111",111,362.946868896484,7.53591728210449,"c57","12_4"
"112",112,279.840484619141,7.72499227523804,"c57","12_4"
"113",113,43.4696846008301,0.145219847559929,"c57","2_8"
"114",114,44.1865348815918,0.129704907536507,"c57","2_8"
"115",115,134.279449462891,0.630734086036682,"c57","2_8"
"116",116,138.789947509766,0.652454972267151,"c57","2_8"
"117",117,208.447219848633,1.33304369449615,"c57","2_8"
"118",118,202.162582397461,1.37296879291534,"c57","2_8"
"119",119,253.445007324219,1.92261147499084,"c57","2_8"
"120",120,255.973876953125,2.01218438148499,"c57","2_8"
"121",121,276.047576904297,2.96128511428833,"c57","2_8"
"122",122,292.277160644531,2.83799290657043,"c57","2_8"
"123",123,361.648315429688,4.35411310195923,"c57","2_8"
"124",124,343.133483886719,4.28233051300049,"c57","2_8"
"125",125,420.476837158203,6.02352046966553,"c57","2_8"
"126",126,330.344818115234,5.90788269042969,"c57","2_8"
"127",127,369.262390136719,7.36670064926147,"c57","2_8"
"128",128,353.220764160156,7.60397577285767,"c57","2_8"
"129",129,24.6443214416504,0.0690932050347328,"c57","3_6"
"130",130,27.0475311279297,0.064335286617279,"c57","3_6"
"131",131,154.345153808594,0.64542156457901,"c57","3_6"
"132",132,163.389480590820,0.611288666725159,"c57","3_6"
"133",133,218.992340087891,1.33842217922211,"c57","3_6"
"134",134,223.767669677734,1.32249355316162,"c57","3_6"
"135",135,306.3505859375,1.92571449279785,"c57","3_6"
"136",136,259.520812988281,1.94888341426849,"c57","3_6"
"137",137,293.181671142578,2.75607419013977,"c57","3_6"
"138",138,286.219024658203,2.67560338973999,"c57","3_6"
"139",139,381.964324951172,3.78337001800537,"c57","3_6"
"140",140,353.642944335938,4.02747201919556,"c57","3_6"
"141",141,419.319305419922,5.51752662658691,"c57","3_6"
"142",142,428.754119873047,5.46808576583862,"c57","3_6"
"143",143,423.506134033203,7.13976907730103,"c57","3_6"
"144",144,378.537811279297,7.22851419448853,"c57","3_6"
"145",145,43.4618949890137,0.155149415135384,"c57","11_1"
"146",146,43.9843521118164,0.159700453281403,"c57","11_1"
"147",147,134.337509155273,0.665073812007904,"c57","11_1"
"148",148,136.9423828125,0.636940062046051,"c57","11_1"
"149",149,186.966949462891,1.32456219196320,"c57","11_1"
"150",150,184.947372436523,1.34090459346771,"c57","11_1"
"151",151,228.087814331055,2.16816115379333,"c57","11_1"
"152",152,225.084869384766,2.02211403846741,"c57","11_1"
"153",153,264.287872314453,2.81254839897156,"c57","11_1"
"154",154,241.713577270508,2.76621055603027,"c57","11_1"
"155",155,298.817077636719,4.46044206619263,"c57","11_1"
"156",156,319.769958496094,4.31377410888672,"c57","11_1"
"157",157,294.957427978516,6.49393367767334,"c57","11_1"
"158",158,274.323547363281,6.65590953826904,"c57","11_1"
"159",159,342.263977050781,7.58453035354614,"c57","11_1"
"160",160,350.169586181641,7.35739183425903,"c57","11_1"
"161",161,81.609748840332,0.314642995595932,"c57","10_5"
"162",162,85.1363830566406,0.321883320808411,"c57","10_5"
"163",163,131.298370361328,0.664246320724487,"c57","10_5"
"164",164,129.688034057617,0.681209325790405,"c57","10_5"
"165",165,187.991760253906,1.43233931064606,"c57","10_5"
"166",166,226.305999755859,1.36097061634064,"c57","10_5"
"167",167,209.508758544922,2.19877743721008,"c57","10_5"
"168",168,220.805053710938,2.1075496673584,"c57","10_5"
"169",169,240.582992553711,2.95445847511292,"c57","10_5"
"170",170,268.121887207031,2.90708613395691,"c57","10_5"
"171",171,251.843154907227,4.4126558303833,"c57","10_5"
"172",172,303.686584472656,4.54918766021729,"c57","10_5"
"173",173,227.569854736328,6.05661916732788,"c57","10_5"
"174",174,359.158538818359,6.13419389724731,"c57","10_5"
"175",175,289.635559082031,7.8149790763855,"c57","10_5"
"176",176,267.569030761719,7.61349201202393,"c57","10_5"
"177",177,68.1292266845703,0.182455703616142,"c57","9_7"
"178",178,61.7317924499512,0.193005859851837,"c57","9_7"
"179",179,161.687805175781,0.713480412960052,"c57","9_7"
"180",180,170.710861206055,0.69010454416275,"c57","9_7"
"181",181,231.013900756836,1.36655604839325,"c57","9_7"
"182",182,243.102920532227,1.42427158355713,"c57","9_7"
"183",183,306.300842285156,2.15305995941162,"c57","9_7"
"184",184,285.980926513672,2.03204345703125,"c57","9_7"
"185",185,316.344787597656,2.77427840232849,"c57","9_7"
"186",186,372.712951660156,2.89529490470886,"c57","9_7"
"187",187,325.51611328125,4.61041975021362,"c57","9_7"
"188",188,327.375091552734,4.49560928344727,"c57","9_7"
"189",189,332.587463378906,5.62240791320801,"c57","9_7"
"190",190,456.245758056641,5.79493379592896,"c57","9_7"
"191",191,320.459655761719,7.48068380355835,"c57","9_7"
"192",192,412.822357177734,7.26781892776489,"c57","9_7"


rm(list=ls())
library(nlme)
load("nlme.rdata")
is.data.frame(nlme)
nlme$pool = factor(nlme$pool)
nlme$strain = factor(nlme$strain)
is.factor(nlme$pool)
is.factor(nlme$strain)
uptake = groupedData(y ~ x | pool, data = nlme)
attach(uptake)

plot(uptake)
plot(uptake, outer=~1)

fm1uptake.nls = nls(y ~ (bmax*x) / (kd + x), data = uptake,  
start=list(bmax=300, kd=2), trace=T)
summary(fm1uptake.nls)

plot(fm1uptake.nls)

plot(fm1uptake.nls, pool ~ resid(.), abline=0)

fm1uptake.lis = nlsList(y ~ SSmicmen(x, bmax, kd), data = uptake)

fm1uptake.lis.noSS = nlsList(y ~ (bmax*x) / (kd + x), data = uptake,  
start = c(bmax = 300, kd = 2))

summary(fm1uptake.lis)
summary(fm1uptake.lis.noSS)

plot(intervals(fm1uptake.lis), layout = c(2,1))
plot(fm1uptake.lis, pool ~ resid(.), abline = 0)

fm1uptake.nlme = nlme(fm1uptake.lis)

# fmuptake.nlme = nlme(y ~ SSmicmen(x, bmax, kd), data = uptake, fixed  
= bmax + kd ~ 1, random = bmax + kd ~ 1, start = fixef(fm1uptake.lis))


summary(fm1uptake.nlme)
summary(fm1uptake.nls)
pairs(fm1uptake.nlme)
fm2uptake.nlme = update(fm1uptake.nlme, random = bmax ~ 1)
anova(fm1uptake.nlme, fm2uptake.nlme)
summary(fm2uptake.nlme)
plot(fm2uptake.nlme)
plot(augPred(fm2uptake.nlme, level = 0:1), layout = c(2,6)) # awesome  
plot!
qqnorm(fm2uptake.nlme, abline = c(0,1))

# but what about strain differences??

plot(uptake, outer = ~strain, layout = c(2,1))

# try to look at strain subset

st = coef(nls(y ~ (bmax*x) / (kd + x), data = uptake,   
start=list(bmax=300, kd=2)))
st

fm1uptake.c57.nlme = nlme(y ~ SSmicmen(x, bmax, kd), data = uptake,  
subset = strain == "c57", fixed = bmax + kd ~ 1, random = bmax + kd ~  
1, start = list(fixed = st))

fm1uptake.dba.nlme = update(fm1uptake.c57.nlme, subset = strain ==  
"dba")

fm1uptake.c57.nlme
fm1uptake.dba.nlme

fm2uptake.c57.nlme = update(fm1uptake.c57.nlme, random = bmax ~ 1)
fm2uptake.dba.nlme = update(fm1uptake.dba.nlme, random = bmax ~ 1)

fm2uptake.c57.nlme
fm2uptake.dba.nlme

# trying to follow MASS4 page 291 to allow strain fixed effects

c1 = c(400, 0, 1, 0)
fm3uptake.nlme = nlme(y ~ SSmicmen(x, bmax, kd), data = uptake, fixed  
= bmax + kd ~ strain, random = bmax ~ 1, start = list(fixed = c1))
summary(fm3uptake.nlme)

fm4uptake.nlme = update(fm3uptake.nlme, fixed = bmax + kd ~ 1, start =  
list(fixed = st))

anova(fm4uptake.nlme, fm3uptake.nlme) # finally!

# fun!