How to obtain the R² and plot curve equation from the object created by the Mostest function

Asked

Viewed 63 times

-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")
r
  • 3

    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.

  • 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 of dput(mod) and test it in a new session of R, without any object in memory. Rotate rm(list = ls()) to wipe the memory of R and try to play your own example to ensure it works for other users.

  • 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)

  • 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.

1 answer

1

Since your question does not specifically depend on your measurements, I will use generic simulated data to facilitate playback by other users:

set.seed(86)
x <- sample(1:20, 100, replace = TRUE)  + rnorm(100, 0, 10)
y <- (-6 + 1.2*x - .5*x^2) + rnorm(100, 0, 20)

Check out Mostest’s help, she makes a call to glm with a quadratic function (y = a + b.x + c.x²).

library(vegan)

most <- MOStest(x, y)

> summary(most)
            Length Class      Mode
isHump        1     -none-     logical
isBracketed   1     -none-     logical
hump          3     -none-     numeric
family       11     family     list
coefficients  3     data.frame list
mod          30     glm        list

> formula(most$mod) # ou most$mod$formula
y ~ x + I(x^2)
<environment: 0x56193a7abfc0>

> coef(most$mod) # ou most$mod$coefficients
(Intercept)           x      I(x^2)
 -5.5864792   1.0073143  -0.4913561

As for R², it does not exist in the object (you can check with str(most)). This is because it only serves as a measure of fit quality when using ordinary minimum squares, which is not the case with Generalized Linear Models.

You can calculate a "pseudo-R²"; in the case of a GLM with Gaussian family and link identity, an equivalent measure is:

> with(summary(most$mod), 1 - deviance/null.deviance)
[1] 0.9764036

But remember that this is an approximation, and that for any error distribution model other than gaussian(link = "identity") she will be skewed.

The package rsq implements several generalized and partial R² measures. You can read more about them at this link (in English). But it is best to simply use Akaike Information Criterion (AIC), which is robust and comparable to any model set to the same data set.

Browser other questions tagged r

You are not signed in. Login or sign up in order to post.