Mathematica file for Midterm Take Home

Thanks to Donna for the format of this page.

Your Name Here
Spring, 2018

Your job:
	1. You are to complete this file, by incorporating two additional “missing” models for kill rate (two others are already determined for you, the linear regression model and the General Power Model): you are to add in the models 
		a. Hollings Type III
		b. The “Rachel Model” -- a square root model
		
		
Note that the independent variable is the ratio of moose to wolves.
	2. At the bottom of this file you will compare all four models, using the residual plots and adjusted RSquared values. 
	3. Note the difference between the adjusted RSquared values produced by the NonlinearModelFit algorithms, and our values discussed in class. Which seem more reasonable? Why?
	4. Once you have these models for kill rate, you will turn your attention to InsightMaker, and try each of these models in a compartmental model you will clone from me, at
	https://insightmaker.com/insight/102161/Final-Midterm-Student-version-of-A-More-Realistic-Model-of-Isle-Royale-Predator-Prey-Interactions

Clear["Global`*"]

mooseWolf=Import["http://ceadserv1.nku.edu/longa/classes/2018spring/mat375/docs/WolfMoose/isleroyale_graph_data_28Dec2011.xlsx"];
Dimensions[mooseWolf]

moose=mooseWolf[[1,Table[i,{i,15,55,1}],3]];
wolves=mooseWolf[[1,Table[i,{i,15,55,1}],2]];
kills=mooseWolf[[1,Table[i,{i,15,55,1}],14]];
ybar=Mean[kills];

ratio=moose/wolves;
data=Transpose[{ratio,kills}];
p1=ListPlot[data,PlotRange->{{0,150},{0,1.5}},AxesLabel->{"moose per wolf","kill rate"}]

Use non-linear regression to fit a Hollings Type III model to the Kill Rate

i. Hollings Type III model: f(x)=((ax^2)/((b^2)+(x^2)))

Clear[a,b]
Manipulate[
Module[
{},
p2=Plot[((a (x^2))/((b^2)+(x^2))),{x,0,150}];
Show[p1,p2]
]
,{{a,.986198},0,150,1}
,{{b,Sqrt[601.468]},0,50,.01}
]

NLM1=NonlinearModelFit[data, ((a (x^2))/((b^2)+(x^2))),{{a,1.5},{b,50}},x]
NLM1["ParameterConfidenceIntervalTable"]
NLM1["ParameterTable"]
rsq1=NLM1["AdjustedRSquared"]
yhat=Table[NLM1[i], {i,ratio}];
residuals=kills-yhat;
rsquared1=1-residuals.residuals/((kills-ybar).(kills-ybar))
n=Length[ratio];
p=2; (* number of parameters estimated *)
rsquaredadj1=1-(n-1)/(n-p) (1-rsquared1)

Show[p1, Plot[NLM1[x],{x,0,150}]]
r1=ListPlot[NLM1["FitResiduals"],PlotLabel->"Hollings Type III Model",AxesLabel->{"Fitted Value", "Residual"},PlotRange->{Automatic,{-.5,.5}}];

Use non-linear regression to fit a root model to the Kill Rate:

i. Root model: f(x)=a x^b

Clear[a,b]
Manipulate[
Module[
{},
p2=Plot[a x^b,{x,0,150}];
Show[p1,p2]
]
,{{a,.09},0,1,.01}
,{{b,.5},0,1,.01}
]

NLM2=NonlinearModelFit[data,a x^b,{{a,.09},{b,0.5}},x]
NLM2["ParameterConfidenceIntervalTable"]
NLM2["ParameterTable"]
rsq2=NLM2["AdjustedRSquared"]
yhat=Table[NLM2[i], {i,ratio}];
residuals=kills-yhat;
rsquared2=1-residuals.residuals/((kills-ybar).(kills-ybar))
n=Length[ratio];
p=2; (* number of parameters estimated *)
rsquaredadj2=1-(n-1)/(n-p) (1-rsquared2)

Show[p1, Plot[NLM2[x],{x,0,150}]]
r2=ListPlot[NLM2["FitResiduals"],PlotLabel->"Power Model",AxesLabel->{"Fitted Value", "Residual"},PlotRange->{Automatic,{-.5,.5}}];

Use linear regression to fit a square root model to the Kill Rate:

i. Square Root model: f(x)=a Sqrt[x]

Clear[a,b]
Manipulate[
Module[
{},
p2=Plot[a Sqrt[x],{x,0,150}];
Show[p1,p2]
]
,{{a,.09},0,1,.01}
]

ii. I actually don’t want an intercept, but what do I find if I do the linear regression with the intercept term?

LM=LinearModelFit[data,{Sqrt[x]},x]
LM["ParameterConfidenceIntervalTable"]
LM["ParameterTable"]
rsq3=LM["AdjustedRSquared"]
LM["ANOVATable"]
yhat=Table[LM[i], {i,ratio}];
residuals=kills-yhat;
rsquared3=1-residuals.residuals/((kills-ybar).(kills-ybar))
n=Length[ratio];
p=2; (* number of parameters estimated *)
rsquaredadj3=1-(n-1)/(n-p) (1-rsquared3)

ii. I actually don’t want an intercept:

LM=LinearModelFit[data,{Sqrt[x]},x,IncludeConstantBasis->False]
LM["ParameterConfidenceIntervalTable"]
LM["ParameterTable"]
rsq3=LM["AdjustedRSquared"]
LM["ANOVATable"]
yhat=Table[LM[i], {i,ratio}];
residuals=kills-yhat;
rsquared3=1-residuals.residuals/((kills-ybar).(kills-ybar))
n=Length[ratio];
p=1; (* number of parameters estimated *)
rsquaredadj3=1-(n-1)/(n-p) (1-rsquared3)

What if we use non-linear regression to fit a square root model to the Kill Rate? Same thing!

NLMsq=NonlinearModelFit[data,a Sqrt[x],{{a,.09}},x]
NLMsq["ParameterConfidenceIntervalTable"]
NLMsq["ParameterTable"]
rsq3=NLMsq["AdjustedRSquared"]
yhat=Table[NLMsq[i], {i,ratio}];
residuals=kills-yhat;
rsquared3=1-residuals.residuals/((kills-ybar).(kills-ybar))
n=Length[ratio];
p=2; (* number of parameters estimated *)
rsquaredadj3=1-(n-1)/(n-p) (1-rsquared3)

Show[p1, Plot[LM[x],{x,0,150}]]
r3=ListPlot[LM["FitResiduals"],PlotLabel->"Power (Root) Model",AxesLabel->{"Fitted Value", "Residual"},PlotRange->{Automatic,{-.5,.5}}];

Use linear regression to fit model to the Kill Rate:

LM2=LinearModelFit[data, {x},x]
LM2["ParameterConfidenceIntervalTable"]
LM2["ParameterTable"]
rsq4=LM2["AdjustedRSquared"]
yhat=Table[LM2[i], {i,ratio}];
residuals=kills-yhat;
rsquared4=1-residuals.residuals/((kills-ybar).(kills-ybar))
n=Length[ratio];
p=2; (* number of parameters estimated *)
rsquaredadj4=1-(n-1)/(n-p) (1-rsquared4)

Show[p1, Plot[LM2[x],{x,0,150}]]
r4=ListPlot[LM2["FitResiduals"],PlotLabel->"Linear Model",AxesLabel->{"Fitted Value", "Residual"},PlotRange->{Automatic,{-.5,.5}}];

Now let’s graph the residuals:

GraphicsGrid[{{r1,r2},{r3,r4}},Frame->True]

Grid[{{"Model","Nominal Adjusted (R^2)","(R^2)","Adjusted (R^2)"},{"Hollings", rsq1,rsquared1, rsquaredadj1},{"Power", rsq2,rsquared2,rsquaredadj2},{"Square Root",rsq3,rsquared3,rsquaredadj3},{"Linear", rsq4,rsquared4,rsquaredadj4}},Frame->All]

This is about the wolves and moose of Isle Royale, Michigan:
http://www.isleroyalewolf.org/data/data/home.html

Andy Long
Spring, 2018

Thanks to:
https://mathematica.stackexchange.com/questions/74214/the-lotka-volterra-predator-prey-model

Perhaps helpful stuff.... from other sources: http://reference.wolfram.com/language/tutorial/NDSolveDelayDifferentialEquations.html

Now let’s make some moose and wolves!
Starting parameters from Wang, et al: MBS6971.pdf

Clear["Global `*"]

Now let’s get the actual wolf/moose data, extracted from a much more extensive data set, available from 
http://www.isleroyalewolf.org/data/data/home.html

IsleRoyale=Import["http://ceadserv1.nku.edu/longa//classes/mat375/data/isleroyale_graph_data_28Dec2011.csv"];

A look at the actual data:

year=Rest[IsleRoyale[[All,1]]];
minYear=Min[year]
year=year-minYear;
wolves=Rest[IsleRoyale[[All,2]]];
moose=Rest[IsleRoyale[[All,3]]];
ListPlot[Transpose[{year+minYear,moose}]]
ListPlot[Transpose[{year+minYear,wolves}]]
(* Do you see a ring-around-the-rosey? *)
ListPlot[Transpose[{moose,wolves}]]

eta=0.9861984930612291;
dee=25;
kappa=2000;
HollingTypeIII[x_]:=eta *x ^2/(x^2 + dee^2);
MoosePerWolfPerYear[x_]:=12*HollingTypeIII[x];

MoosePerWolfPerYear[x_]:=12*LM[x]; (* Square Root Model *)
MoosePerWolfPerYear[x_]:=12*NLM1[x]; (* Hollings III Model *)
MoosePerWolfPerYear[x_]:=12*NLM2[x]; (* Power Model *)
MoosePerWolfPerYear[x_]:=12*LM2[x]; (* Linear Model *)

p2=Plot[HollingTypeIII[x],{x,0, 150}];
Show[p1,p2]
Manipulate[
lv=
NDSolve[
{(*  moose: *)
m'[t]==bm m[t](1-m[t]/kappa)-MoosePerWolfPerYear[ m[t]/w[t]]w[t] ,
(* wolves: *)
w'[t]==bw*w[t]*(wpm*MoosePerWolfPerYear[ m[t]/w[t]])-deathRateWolf w[t],
m[0]==m0,
w[0]==w0
},
{m[t],w[t]},
{t,0,tend}];
plotMoose=ListPlot[Transpose[{year,moose/m0}],PlotMarkers->{{▲,10}},PlotStyle->{Blue}];
plotWolves=ListPlot[Transpose[{year,wolves/w0}],PlotMarkers->{{▲,10}},PlotStyle->{Red}];
p2=Plot[Evaluate[{m[t]/m0,w[t]/w0}/.lv],{t,0,tend},
PlotRange->All,AxesOrigin->{0,0},PlotLabel->Evaluate[{m[tend]}/.lv]];
Show[p2,plotMoose,plotWolves],
{{lv,{}},None},
{{tend,60},2,1000,10},
{{m0,563},10,3000,10},
{{w0,20},5,50,1},
{{bm,0.31},0.01,1.7,.1},
{{bw,0.769},.001,1},
{{wpm,.051},.001,.1,.01},
{{deathRateWolf,.39},.1,1,.01}
]