I have a data.frame(x) and I am doing a GAM analysis.
x = data[, c("asthma3", "TargetGroup2012", "age3", "SmokingGroup_Kai", "bmi3", "PA_Score")]
x = na.omit(x)
x = data.frame(scale(x))
g = gam(formula=asthma3~s(age3)+s(bmi3)+SmokingGroup_Kai+PA_Score+TargetGroup2012, data=x)
plot(g, se=T)
The produced plot has the function s and its SE. I would like to save in a data.frame the value of s and its SE for each variable.
With p = plot(g, se=T) and str(p) I can see the values S(x) and SE but I do not know how to access to them..
> p = plot(g, se=T)
> p
Call:
gam(formula = asthma3 ~ s(age3) + s(bmi3) + SmokingGroup_Kai +
PA_Score + TargetGroup2012, data = x)
Degrees of Freedom: 3185 total; 3174 Residual
Residual Deviance: 258.4571
> str(p)
List of 35
$ smooth.frame :'data.frame': 3186 obs. of 2 variables:
..$ s(age3):Class 'smooth' atomic [1:3186] 55 49 58 51 52 62 62 54 49 54 ...
.. .. ..- attr(*, "spar")= num 1
.. .. ..- attr(*, "df")= num 4
.. .. ..- attr(*, "call")= language gam.s(data[["s(age3)"]], z, w, spar = 1, df = 4)
..$ s(bmi3):Class 'smooth' atomic [1:3186] 20.9 20.2 30.5 34.1 23.2 ...
.. .. ..- attr(*, "spar")= num 1
.. .. ..- attr(*, "df")= num 4
.. .. ..- attr(*, "call")= language gam.s(data[["s(bmi3)"]], z, w, spar = 1, df = 4)
$ coefficients : Named num [1:6] -0.051222 -0.001759 0.005935 0.010658 0.000831 ...
..- attr(*, "names")= chr [1:6] "(Intercept)" "s(age3)" "s(bmi3)" "SmokingGroup_Kai" ...
$ residuals : Named num [1:3186] -0.0705 -0.0954 -0.1249 -0.171 -0.1079 ...
..- attr(*, "names")= chr [1:3186] "1" "2" "3" "4" ...
$ fitted.values : Named num [1:3186] 0.0705 0.0954 0.1249 0.171 0.1079 ...
..- attr(*, "names")= chr [1:3186] "1" "2" "3" "4" ...
$ effects : Named num [1:3186] -5.13778 0.00287 -1.58832 -0.97786 -0.04935 ...
..- attr(*, "names")= chr [1:3186] "(Intercept)" "s(age3)" "s(bmi3)" "SmokingGroup_Kai" ...
$ weights : Named num [1:3186] 1 1 1 1 1 1 1 1 1 1 ...
..- attr(*, "names")= chr [1:3186] "1" "2" "3" "4" ...
$ rank : int 6
$ assign : int [1:6] 0 1 2 3 4 5
$ qr :List of 5
..$ qr : num [1:3186, 1:6] -56.4447 0.0177 0.0177 0.0177 0.0177 ...
.. ..- attr(*, "dimnames")=List of 2
.. .. ..$ : chr [1:3186] "1" "2" "3" "4" ...
.. .. ..$ : chr [1:6] "(Intercept)" "s(age3)" "s(bmi3)" "SmokingGroup_Kai" ...
.. ..- attr(*, "assign")= int [1:6] 0 1 2 3 4 5
..$ qraux: num [1:6] 1.02 1.02 1.02 1.02 1.02 ...
..$ rank : int 6
..$ pivot: int [1:6] 1 2 3 4 5 6
..$ tol : num 1e-07
..- attr(*, "class")= chr "qr"
$ smooth : num [1:3186, 1:2] -0.00618 0.00825 -0.00368 0.0094 0.00506 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : chr [1:3186] "1" "2" "3" "4" ...
.. ..$ : chr [1:2] "s(age3)" "s(bmi3)"
$ nl.df : Named num [1:2] 3 3
..- attr(*, "names")= chr [1:2] "s(age3)" "s(bmi3)"
$ df.residual : num 3174
$ var : num [1:3186, 1:2] 0.000846 0.000722 0.000766 0.000798 0.00082 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : chr [1:3186] "1" "2" "3" "4" ...
.. ..$ : chr [1:2] "s(age3)" "s(bmi3)"
$ additive.predictors: Named num [1:3186] 0.0705 0.0954 0.1249 0.171 0.1079 ...
..- attr(*, "names")= chr [1:3186] "1" "2" "3" "4" ...
$ R : num [1:6, 1:6] -56.4 0 0 0 0 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : chr [1:6] "(Intercept)" "s(age3)" "s(bmi3)" "SmokingGroup_Kai" ...
.. ..$ : chr [1:6] "(Intercept)" "s(age3)" "s(bmi3)" "SmokingGroup_Kai" ...
$ rank : int 6
$ family :List of 11
..$ family : chr "gaussian"
..$ link : chr "identity"
..$ linkfun :function (mu)
..$ linkinv :function (eta)
..$ variance :function (mu)
..$ dev.resids:function (y, mu, wt)
..$ aic :function (y, n, mu, wt, dev)
..$ mu.eta :function (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)
..$ valideta :function (eta)
..- attr(*, "class")= chr "family"
$ deviance : num 258
$ aic : num 1065
$ null.deviance : num 264
$ iter : int 2
$ prior.weights : num [1:3186] 1 1 1 1 1 1 1 1 1 1 ...
$ y : Named num [1:3186] 0 0 0 0 0 0 0 0 0 0 ...
..- attr(*, "names")= chr [1:3186] "1" "2" "3" "4" ...
$ df.null : int 3185
$ nl.chisq : Named num [1:2] 0.399 0.691
..- attr(*, "names")= chr [1:2] "s(age3)" "s(bmi3)"
$ call : language gam(formula = asthma3 ~ s(age3) + s(bmi3) + SmokingGroup_Kai + PA_Score + TargetGroup2012, data = x)
$ formula :Class 'formula' length 3 asthma3 ~ s(age3) + s(bmi3) + SmokingGroup_Kai + PA_Score + TargetGroup2012
.. ..- attr(*, ".Environment")=<environment: R_GlobalEnv>
$ terms :Classes 'terms', 'formula' length 3 asthma3 ~ s(age3) + s(bmi3) + SmokingGroup_Kai + PA_Score + TargetGroup2012
.. ..- attr(*, "variables")= language list(asthma3, s(age3), s(bmi3), SmokingGroup_Kai, PA_Score, TargetGroup2012)
.. ..- attr(*, "factors")= int [1:6, 1:5] 0 1 0 0 0 0 0 0 1 0 ...
.. .. ..- attr(*, "dimnames")=List of 2
.. .. .. ..$ : chr [1:6] "asthma3" "s(age3)" "s(bmi3)" "SmokingGroup_Kai" ...
.. .. .. ..$ : chr [1:5] "s(age3)" "s(bmi3)" "SmokingGroup_Kai" "PA_Score" ...
.. ..- attr(*, "term.labels")= chr [1:5] "s(age3)" "s(bmi3)" "SmokingGroup_Kai" "PA_Score" ...
.. ..- attr(*, "specials")=Dotted pair list of 3
.. .. ..$ s : int [1:2] 2 3
.. .. ..$ lo : NULL
.. .. ..$ random: NULL
.. ..- attr(*, "order")= int [1:5] 1 1 1 1 1
.. ..- attr(*, "intercept")= int 1
.. ..- attr(*, "response")= int 1
.. ..- attr(*, ".Environment")=<environment: R_GlobalEnv>
$ data :'data.frame': 3186 obs. of 6 variables:
..$ asthma3 : num [1:3186] 0 0 0 0 0 0 0 0 0 0 ...
..$ TargetGroup2012 : num [1:3186] 2 2 2 2 3 3 3 2 2 3 ...
..$ age3 : num [1:3186] 55 49 58 51 52 62 62 54 49 54 ...
..$ SmokingGroup_Kai: num [1:3186] 4 4 5 4 3 3 3 2 1 3 ...
..$ bmi3 : num [1:3186] 20.9 20.2 30.5 34.1 23.2 ...
..$ PA_Score : num [1:3186] 2 3 3 2 3 2 2 3 3 2 ...
..- attr(*, "na.action")=Class 'omit' Named int [1:77] 28 74 75 151 179 185 242 267 287 297 ...
.. .. ..- attr(*, "names")= chr [1:77] "28" "74" "75" "151" ...
$ offset : NULL
$ control :List of 5
..$ epsilon : num 1e-07
..$ maxit : num 30
..$ bf.epsilon: num 1e-07
..$ bf.maxit : num 30
..$ trace : logi FALSE
$ method : chr "glm.fit"
$ contrasts : NULL
$ xlevels : Named list()
$ preplot :List of 5
..$ s(age3) :List of 5
.. ..$ x : num [1:3186] 55 49 58 51 52 62 62 54 49 54 ...
.. ..$ y : Named num [1:3186] -0.00801 0.01698 -0.01078 0.01461 0.00851 ...
.. .. ..- attr(*, "names")= chr [1:3186] "1" "2" "3" "4" ...
.. ..$ se.y: Named num [1:3186] 0.00837 0.00929 0.00898 0.00865 0.00843 ...
.. .. ..- attr(*, "names")= chr [1:3186] "1" "2" "3" "4" ...
.. ..$ xlab: chr "age3"
.. ..$ ylab: chr "s(age3)"
.. ..- attr(*, "class")= chr "preplot.gam"
..$ s(bmi3) :List of 5
.. ..$ x : num [1:3186] 20.9 20.2 30.5 34.1 23.2 ...
.. ..$ y : Named num [1:3186] -0.0213 -0.0222 0.0244 0.0566 -0.0152 ...
.. .. ..- attr(*, "names")= chr [1:3186] "1" "2" "3" "4" ...
.. ..$ se.y: Named num [1:3186] 0.00737 0.00917 0.00919 0.01372 0.005 ...
.. .. ..- attr(*, "names")= chr [1:3186] "1" "2" "3" "4" ...
.. ..$ xlab: chr "bmi3"
.. ..$ ylab: chr "s(bmi3)"
.. ..- attr(*, "class")= chr "preplot.gam"
..$ SmokingGroup_Kai:List of 5
.. ..$ x : num [1:3186] 4 4 5 4 3 3 3 2 1 3 ...
.. ..$ y : Named num [1:3186] 0.01692 0.01692 0.02758 0.01692 0.00626 ...
.. .. ..- attr(*, "names")= chr [1:3186] "1" "2" "3" "4" ...
.. ..$ se.y: Named num [1:3186] 0.00516 0.00516 0.00842 0.00516 0.00191 ...
.. .. ..- attr(*, "names")= chr [1:3186] "1" "2" "3" "4" ...
.. ..$ xlab: chr "SmokingGroup_Kai"
.. ..$ ylab: chr "partial for SmokingGroup_Kai"
.. ..- attr(*, "class")= chr "preplot.gam"
..$ PA_Score :List of 5
.. ..$ x : num [1:3186] 2 3 3 2 3 2 2 3 3 2 ...
.. ..$ y : Named num [1:3186] -0.000292 0.000539 0.000539 -0.000292 0.000539 ...
.. .. ..- attr(*, "names")= chr [1:3186] "1" "2" "3" "4" ...
.. ..$ se.y: Named num [1:3186] 0.003 0.00553 0.00553 0.003 0.00553 ...
.. .. ..- attr(*, "names")= chr [1:3186] "1" "2" "3" "4" ...
.. ..$ xlab: chr "PA_Score"
.. ..$ ylab: chr "partial for PA_Score"
.. ..- attr(*, "class")= chr "preplot.gam"
..$ TargetGroup2012 :List of 5
.. ..$ x : num [1:3186] 2 2 2 2 3 3 3 2 2 3 ...
.. ..$ y : Named num [1:3186] -0.00786 -0.00786 -0.00786 -0.00786 0.01679 ...
.. .. ..- attr(*, "names")= chr [1:3186] "1" "2" "3" "4" ...
.. ..$ se.y: Named num [1:3186] 0.00289 0.00289 0.00289 0.00289 0.00617 ...
.. .. ..- attr(*, "names")= chr [1:3186] "1" "2" "3" "4" ...
.. ..$ xlab: chr "TargetGroup2012"
.. ..$ ylab: chr "partial for TargetGroup2012"
.. ..- attr(*, "class")= chr "preplot.gam"
..- attr(*, "class")= chr "preplot.gam"
- attr(*, "class")= chr [1:3] "gam" "glm" "lm"
> p$preplot$s(age3)$y
Error: attempt to apply non-function
> p$preplot$s(age3)$ y
Error: attempt to apply non-function
How can I access to value SE and y for each variables inside plot?
The problem is that the names of the kind s(age3) are not valid R names. So when R sees that it thinks you are trying to evaluate the function p$preplot$s on a variable called age3. You can do this to retreive these values:
p[[c("preplot", "s(age3)", "y")]]
Hope this works.