-2
In the R, I surrounded the function Mostest of vegan package to verify whether the relationship between two variables had a hump shape. So far all right, I checked that in fact the curve had this format and generated the graph with the Plot function. Now I need to get the R² of quadratic adjustment and the equation of the curve, how can I do? It has how to get the equation based on the data?
I did not find this information giving a Summary. This is the script:
mod <- MOStest(PET, Riqueza, family=gaussian)
summary(mod)
plot(mod)
As requested, in order to allow replication, I reproduce here the result of the dput function on the "mod" object generated in the above script. I reiterate that this is not a mistake, the script has worked perfectly, my problem is how to get two pieces of information in the results, R2 and the curve equation. Follow the result of dput:
dput(mod)
structure(list(isHump = TRUE, isBracketed = TRUE, hump = c(min = 1087,
hump = 1258.31714127466, max = 1750), family = structure(list(
family = "gaussian", link = "identity", linkfun = function (mu)
mu, linkinv = function (eta)
eta, variance = function (mu)
rep.int(1, length(mu)), dev.resids = function (y, mu, wt)
wt * ((y - mu)^2), aic = function (y, n, mu, wt, dev)
{
nobs <- length(y)
nobs * (log(dev/nobs * 2 * pi) + 1) + 2 - sum(log(wt))
}, mu.eta = function (eta)
rep.int(1, length(eta)), initialize = expression({
n <- rep.int(1, nobs)
if (is.null(etastart) && is.null(start) && is.null(mustart) &&
((family$link == "inverse" && any(y == 0)) || (family$link ==
"log" && any(y <= 0))))
stop("cannot find valid starting values: please specify some")
mustart <- y
}), validmu = function (mu)
TRUE, valideta = function (eta)
TRUE), class = "family"), coefficients = structure(list(`min/max` = c(1087,
1750, NA), F = c(1.47951715715753, 8.32257715767822, NA), `Pr(>F)` = c(0.240469622558875,
0.0102866435223026, 0.248282640795972)), row.names = c("hump at min",
"hump at max", "Combined"), class = "data.frame"), mod = structure(list(
coefficients = c(`(Intercept)` = -610.411556247032, x = 1.27173798958452,
`I(x^2)` = -0.000505332856030344), residuals = c(`1` = -19.0903389236728,
`2` = 36.2323157794438, `3` = 32.4967819365999, `4` = -28.6873895737368,
`5` = 16.1646964784381, `6` = 10.7498764032539, `7` = 52.5163341100767,
`8` = -4.88200106442309, `9` = 23.1179989355769, `10` = -29.4749022142456,
`11` = 26.47737929651, `12` = 33.9294745606495, `13` = -18.8467120784045,
`14` = 24.957850571315, `15` = -17.8286205965798, `16` = -9.1830694159712,
`17` = -10.6595101557635, `18` = -84.994718570696, `19` = -47.4473915454337,
`20` = 14.451946067051), fitted.values = c(`1` = 132.090338923673,
`2` = 165.767684220556, `3` = 176.5032180634, `4` = 187.687389573737,
`5` = 188.835303521562, `6` = 141.250123596746, `7` = 189.483665889923,
`8` = 174.882001064423, `9` = 174.882001064423, `10` = 181.474902214246,
`11` = 183.52262070349, `12` = 187.07052543935, `13` = 186.846712078404,
`14` = 179.042149428685, `15` = 184.82862059658, `16` = 185.183069415971,
`17` = 189.659510155764, `18` = 187.994718570696, `19` = 133.447391545434,
`20` = 67.548053932949), effects = c(`(Intercept)` = -759.815898754429,
x = 106.275663892188, `I(x^2)` = -78.0307813641536, -19.2473885023306,
16.6345964780407, 6.14597586854161, 57.896572949065, 17.8854812803235,
45.8854812803235, -33.4235587888267, 23.1930390208033, 44.3212539743377,
-8.13125523356534, 44.136777943802, -4.52502693883779, -11.7480403942007,
-6.23860669678687, -76.0710038322058, -51.2562947790223,
21.0033178425398), R = structure(c(-4.47213595499958, 0,
0, -5926.92178116094, -833.327186643997, 0, -8010229.64428267,
-2307488.00713431, 154414.62084442), .Dim = c(3L, 3L), .Dimnames = list(
c("(Intercept)", "x", "I(x^2)"), c("(Intercept)", "x",
"I(x^2)"))), rank = 3L, qr = structure(list(qr = structure(c(-4.47213595499958,
0.223606797749979, 0.223606797749979, 0.223606797749979,
0.223606797749979, 0.223606797749979, 0.223606797749979,
0.223606797749979, 0.223606797749979, 0.223606797749979,
0.223606797749979, 0.223606797749979, 0.223606797749979,
0.223606797749979, 0.223606797749979, 0.223606797749979,
0.223606797749979, 0.223606797749979, 0.223606797749979,
0.223606797749979, -5926.92178116094, -833.327186643997,
0.0542778401641913, -0.215724151377843, -0.0897232219915606,
0.231879150156285, -0.16532377962333, -0.34532510731802,
-0.34532510731802, 0.0134775392200616, -0.00692261125200326,
-0.226524231039525, -0.230124257593419, -0.314124877184274,
-0.257724461173271, -0.0261227528727702, -0.152123682259053,
-0.209724107121354, 0.260679362587436, 0.450280761092509,
-8010229.64428267, -2307488.00713431, 154414.62084442, 0.0613632736741654,
0.225096497148415, -0.00436787233684219, 0.143992246106897,
-0.256023572028568, -0.256023572028568, 0.25282586685686,
0.254940910739639, 0.0406845488494826, 0.0335585020297847,
-0.165808644983827, -0.0249472138599037, 0.253512174952593,
0.161857618246654, 0.0723981291840941, -0.0703012237387962,
-0.690588342533872), .Dim = c(20L, 3L), .Dimnames = list(
c("1", "2", "3", "4", "5", "6", "7", "8", "9", "10",
"11", "12", "13", "14", "15", "16", "17", "18", "19",
"20"), c("(Intercept)", "x", "I(x^2)"))), rank = 3L,
qraux = c(1.22360679774998, 1.12147833583688, 1.23736627251662
), pivot = 1:3, tol = 1e-11), class = "qr"), family = structure(list(
family = "gaussian", link = "identity", linkfun = function (mu)
mu, linkinv = function (eta)
eta, variance = function (mu)
rep.int(1, length(mu)), dev.resids = function (y, mu,
wt)
wt * ((y - mu)^2), aic = function (y, n, mu, wt, dev)
{
nobs <- length(y)
nobs * (log(dev/nobs * 2 * pi) + 1) + 2 - sum(log(wt))
}, mu.eta = function (eta)
rep.int(1, length(eta)), initialize = expression({
n <- rep.int(1, nobs)
if (is.null(etastart) && is.null(start) && is.null(mustart) &&
((family$link == "inverse" && any(y == 0)) ||
(family$link == "log" && any(y <= 0))))
stop("cannot find valid starting values: please specify some")
mustart <- y
}), validmu = function (mu)
TRUE, valideta = function (eta)
TRUE), class = "family"), linear.predictors = c(`1` = 132.090338923673,
`2` = 165.767684220556, `3` = 176.5032180634, `4` = 187.687389573737,
`5` = 188.835303521562, `6` = 141.250123596746, `7` = 189.483665889923,
`8` = 174.882001064423, `9` = 174.882001064423, `10` = 181.474902214246,
`11` = 183.52262070349, `12` = 187.07052543935, `13` = 186.846712078404,
`14` = 179.042149428685, `15` = 184.82862059658, `16` = 185.183069415971,
`17` = 189.659510155764, `18` = 187.994718570696, `19` = 133.447391545434,
`20` = 67.548053932949), deviance = 21148.4804239745, aic = 204.029362543648,
null.deviance = 38531.8, iter = 2L, weights = c(`1` = 1,
`2` = 1, `3` = 1, `4` = 1, `5` = 1, `6` = 1, `7` = 1, `8` = 1,
`9` = 1, `10` = 1, `11` = 1, `12` = 1, `13` = 1, `14` = 1,
`15` = 1, `16` = 1, `17` = 1, `18` = 1, `19` = 1, `20` = 1
), prior.weights = c(`1` = 1, `2` = 1, `3` = 1, `4` = 1,
`5` = 1, `6` = 1, `7` = 1, `8` = 1, `9` = 1, `10` = 1, `11` = 1,
`12` = 1, `13` = 1, `14` = 1, `15` = 1, `16` = 1, `17` = 1,
`18` = 1, `19` = 1, `20` = 1), df.residual = 17L, df.null = 19L,
y = c(`1` = 113L, `2` = 202L, `3` = 209L, `4` = 159L, `5` = 205L,
`6` = 152L, `7` = 242L, `8` = 170L, `9` = 198L, `10` = 152L,
`11` = 210L, `12` = 221L, `13` = 168L, `14` = 204L, `15` = 167L,
`16` = 176L, `17` = 179L, `18` = 103L, `19` = 86L, `20` = 82L
), converged = TRUE, boundary = FALSE, model = structure(list(
y = c(113L, 202L, 209L, 159L, 205L, 152L, 242L, 170L,
198L, 152L, 210L, 221L, 168L, 204L, 167L, 176L, 179L,
103L, 86L, 82L), x = c(1596L, 1476L, 1420L, 1195L, 1300L,
1568L, 1237L, 1087L, 1087L, 1386L, 1369L, 1186L, 1183L,
1113L, 1160L, 1353L, 1248L, 1200L, 1592L, 1750L), `I(x^2)` = structure(c(2547216,
2178576, 2016400, 1428025, 1690000, 2458624, 1530169,
1181569, 1181569, 1920996, 1874161, 1406596, 1399489,
1238769, 1345600, 1830609, 1557504, 1440000, 2534464,
3062500), class = "AsIs")), terms = y ~ x + I(x^2), row.names = c(NA,
20L), class = "data.frame"), call = glm(formula = y ~ x +
I(x^2), family = ..1), formula = y ~ x + I(x^2), terms = y ~
x + I(x^2), data = <environment>, offset = NULL, control = list(
epsilon = 1e-08, maxit = 25, trace = FALSE), method = "glm.fit",
contrasts = NULL, xlevels = structure(list(), .Names = character(0))), class = c("glm",
"lm"))), class = "MOStest")
Welcome to Stackoverflow! Unfortunately, this question cannot be reproduced by anyone trying to answer it. Please take a look at this link (mainly in the use of function
dput
) and see how to ask a reproducible question in R. So, people who wish to help you will be able to do this in the best possible way.– Marcus Nunes
The shared data does not allow the reproduction of the object. I am receiving the error
Error: unexpected ')' in " "
when I try to replicate on my computer. I recommend re-copying the output ofdput(mod)
and test it in a new session ofR
, without any object in memory. Rotaterm(list = ls())
to wipe the memory ofR
and try to play your own example to ensure it works for other users.– Marcus Nunes
Marcus, thank you for your answer. I couldn’t get the output to work, as the analysis runs with a small database, it’s easier for me to post the data here to allow playback: library(vegan) X <- c(1596, 1476, 1420, 1195, 1300, 1568, 1237, 1087, 1386, 1369, 1186, 1183, 1113, 1160, 1353, 1248, 1200, 1592, 1750) Y <- c(113, 202, 209, 159, 205, 152, 242, 170, 198, 152, 210, 221, 168, 204, 167, 176, 179, 103, 86, 82) modXY <- Mostest(X, Y, family=Gaussian) ##To join the 4 charts into a single picture op <- par(mfrow=c(2,2), mar=c(4,4,1,1)+. 1) Plot(modXY)
– Roberto
Marcus, when I copied and pasted the above script, I don’t know why he put it together by ignoring the paragraph entries, so it got a disorganized aspect, and it’s not like that in the original script. I don’t know why it happened, I tried to tidy up but I couldn’t make the line changes.
– Roberto