Una delle assunzioni di base della modellazione lineare è la varianza costante, o omogenea. Cosa significa esattamente? Simuliamo alcuni dati che soddisfano questa condizione per illustrare il concetto.
Di seguito creiamo un vettore ordinato di numeri che vanno da 1 a 10 chiamato x, e poi creiamo un vettore di numeri chiamato y che è una funzione di x. Quando tracciamo x contro y, otteniamo una linea retta con un’intercetta di 1,2 e una pendenza di 2,1.
x <- seq(1,10, length.out = 100)y <- 1.2 + 2.1 * xplot(x, y)
Ora aggiungiamo un po’ di “rumore” ai nostri dati in modo che y non sia completamente determinato da x. Possiamo farlo pescando a caso valori da una distribuzione normale teorica con media 0 e una certa varianza impostata, e poi aggiungendoli alla formula che genera y. La funzione rnorm in R ci permette di farlo facilmente. Di seguito peschiamo 100 valori casuali da una distribuzione normale con media 0 e deviazione standard 2 e li salviamo come un vettore chiamato noise. (Ricordiamo che la deviazione standard è semplicemente la radice quadrata della varianza.) Poi generiamo y con il rumore aggiunto. La funzione set.seed(222) vi permette di ottenere gli stessi dati “casuali” nel caso vogliate seguirli. Infine ri-tracciamo i dati.
set.seed(222)noise <- rnorm(n = 100, mean = 0, sd = 2)y <- 1.2 + 2.1 * x + noiseplot(x, y)
Ora abbiamo dati che sono una combinazione di una componente lineare e deterministica (\(y = 1,2 + 2,1x\)) e rumore casuale tratto da una distribuzione \(N(0, 2)\). Queste sono le ipotesi di base che facciamo sui nostri dati quando adattiamo un modello lineare tradizionale. Di seguito usiamo la funzione lm per “recuperare” i “veri” valori che abbiamo usato per generare i dati, che sono tre:
- L’intercetta: 1,2
- La pendenza: 2.1
- La deviazione standard: 2
m <- lm(y ~ x)summary(m)## ## Call:## lm(formula = y ~ x)## ## Residuals:## Min 1Q Median 3Q Max ## -5.5831 -1.2165 0.3288 1.3022 4.3714 ## ## Coefficients:## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 1.27426 0.44720 2.849 0.00534 ** ## x 2.09449 0.07338 28.541 < 2e-16 ***## ---## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## ## Residual standard error: 1.926 on 98 degrees of freedom## Multiple R-squared: 0.8926, Adjusted R-squared: 0.8915 ## F-statistic: 814.6 on 1 and 98 DF, p-value: < 2.2e-16Le stime (Intercept) e x (cioè la pendenza) nella sezione Coefficienti, 1,27 e 2,09, sono abbastanza vicine a 1,2 e 2,1, rispettivamente. L’errore standard residuo di 1,926 è anche abbastanza vicino al valore costante di 2. Abbiamo prodotto un modello “buono” perché sapevamo come y è stato generato e abbiamo dato alla funzione lm il modello “corretto” da adattare. Anche se non possiamo fare questo nella vita reale, è un esercizio utile per aiutarci a capire i presupposti della modellazione lineare.
E se la varianza non fosse costante? Cosa succederebbe se moltiplicassimo la deviazione standard di 2 per la radice quadrata di x? Come vediamo nel grafico qui sotto, la dispersione verticale dei punti aumenta all’aumentare di x.
set.seed(222)noise <- rnorm(n = 100, mean = 0, sd = 2 * sqrt(x))y <- 1.2 + 2.1 * x + noiseplot(x, y)
Abbiamo moltiplicato 2 per sqrt(x) perché stiamo specificando la deviazione standard. Se potessimo specificare la varianza, avremmo moltiplicato 4 solo per x.
Se adattiamo lo stesso modello usando lm, otteniamo un errore standard residuo di 4,488.
m2 <- lm(y ~ x)summary(m2)## ## Call:## lm(formula = y ~ x)## ## Residuals:## Min 1Q Median 3Q Max ## -15.0460 -2.4013 0.5638 2.8734 10.2996 ## ## Coefficients:## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 1.313 1.042 1.26 0.211 ## x 2.096 0.171 12.26 <2e-16 ***## ---## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## ## Residual standard error: 4.488 on 98 degrees of freedom## Multiple R-squared: 0.6051, Adjusted R-squared: 0.6011 ## F-statistic: 150.2 on 1 and 98 DF, p-value: < 2.2e-16Sappiamo che non è giusto, perché abbiamo simulato i dati. Non c’era una deviazione standard costante quando abbiamo creato il “rumore”. Ogni valore casuale è stato estratto da una diversa distribuzione normale, ognuna con media 0 ma una deviazione standard che variava in base a x. Questo significa che la nostra assunzione di varianza costante è violata. Come potremmo rilevarlo nella vita reale?
Il modo più comune è tracciare i residui rispetto ai valori montati. Questo è facile da fare in R. Basta chiamare plot sull’oggetto modello. Questo genera quattro diversi grafici per valutare le ipotesi di modellazione tradizionali. Si veda questo post sul blog per maggiori informazioni. I grafici che ci interessano sono il primo e il terzo, che possiamo specificare con l’argomento which.
plot(m2, which = 1)
plot(m2, which = 3)
Nel primo grafico vediamo la variabilità intorno a 0 aumentare man mano che ci spostiamo più a destra con valori montati maggiori. Anche nel terzo grafico vediamo aumentare la variabilità man mano che ci spostiamo verso destra, anche se questa volta i residui sono stati standardizzati e trasformati nella scala della radice quadrata. La traiettoria positiva della linea rossa liscia indica un aumento della varianza.
Ora che abbiamo confermato che la nostra ipotesi di varianza non costante non è valida, cosa possiamo fare? Un approccio è quello di trasformare logicamente i dati. Questo a volte funziona se la vostra variabile di risposta è positiva e altamente asimmetrica. Ma non è questo il caso. y è solo leggermente asimmetrico. (Chiamate hist() su y per verificare.) Inoltre sappiamo che i nostri dati non sono su una scala logaritmica.
Un altro approccio è quello di modellare la varianza non costante. Questo è l’argomento di questo post del blog.
Per fare questo, useremo le funzioni del pacchetto nlme che è incluso nell’installazione base di R. La funzione cavallo di battaglia è gls, che sta per “minimi quadrati generalizzati”. La usiamo proprio come faremmo con la funzione lm, eccetto che usiamo anche l’argomento weights insieme a una manciata di funzioni di varianza. Dimostriamolo con un esempio.
Di seguito adattiamo il modello “corretto” ai nostri dati che presentano una varianza non costante. Carichiamo il pacchetto nlme in modo da poter utilizzare la funzione gls. Specifichiamo la sintassi del modello come prima, y ~ x. Poi usiamo l’argomento weights per specificare la funzione varianza, in questo caso varFixed, anch’essa parte del pacchetto nlme. Questo dice che la nostra funzione di varianza non ha parametri e una sola covariata, x, che è esattamente come abbiamo generato i dati. La sintassi, ~x, è una formula unilaterale che può essere letta come “varianza del modello in funzione di x.”
library(nlme)vm1 <- gls(y ~ x, weights = varFixed(~x))summary(vm1)## Generalized least squares fit by REML## Model: y ~ x ## Data: NULL ## AIC BIC logLik## 576.2928 584.0477 -285.1464## ## Variance function:## Structure: fixed weights## Formula: ~x ## ## Coefficients:## Value Std.Error t-value p-value## (Intercept) 1.369583 0.6936599 1.974431 0.0511## x 2.085425 0.1504863 13.857900 0.0000## ## Correlation: ## (Intr)## x -0.838## ## Standardized residuals:## Min Q1 Med Q3 Max ## -2.8942967 -0.6293867 0.1551594 0.6758773 2.2722755 ## ## Residual standard error: 1.925342 ## Degrees of freedom: 100 total; 98 residualIl risultato produce un errore standard residuo di 1,925 che è molto vicino a 2, il valore “vero” che abbiamo usato per generare i dati. Anche i valori di (intercetta) e pendenza, 1,37 e 2,09, sono molto vicini a 1,2 e 2,1.
Ancora una volta possiamo chiamare plot sull’oggetto modello. In questo caso viene generato solo un grafico: residui standardizzati rispetto ai valori montati:
plot(vm1)
Questo grafico sembra buono. Finché modelliamo la nostra varianza come una funzione di x, il modello adattato non si adatta né si sottomette in alcun modo sistematico (a differenza di quando abbiamo usato lm per adattare il modello e assunto una varianza costante. L’idea è che più alto è il peso di una data osservazione, più piccola è la varianza della distribuzione normale da cui è stata tratta la sua componente di rumore. Nel nostro esempio, all’aumentare di x aumenta la varianza. Quindi pesi più alti dovrebbero essere associati a valori x più piccoli. Questo può essere realizzato prendendo il reciproco di x, cioè \(w_i = 1/x\). Così quando \(x = 2\), \(w = 1/2\). Quando \(x = 10\), \(w = 1/10\). Più grandi x ottengono pesi più piccoli.
Infine, per assicurare che pesi più grandi siano associati a varianze più piccole, dividiamo la varianza costante per il peso. Detto matematicamente,
\
o prendendo la radice quadrata per esprimere la deviazione standard,
\
Così, più grande è il denominatore (cioè, più grande è il peso e quindi, più piccolo è il x), più piccola è la varianza e più precisa è la y osservata.
Incidentalmente possiamo usare anche lm per pesare le osservazioni. Anch’essa ha un argomento weights come la funzione gls. L’unica differenza è che dobbiamo essere più espliciti su come esprimiamo i pesi. In questo caso, dobbiamo specificare il reciproco di x. Notate che il risultato è quasi identico a quello che otteniamo usando gls e varFixed.
m3 <- lm(y ~ x, weights = 1/x)summary(m3)## ## Call:## lm(formula = y ~ x, weights = 1/x)## ## Weighted Residuals:## Min 1Q Median 3Q Max ## -5.5725 -1.2118 0.2987 1.3013 4.3749 ## ## Coefficients:## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 1.3696 0.6937 1.974 0.0511 . ## x 2.0854 0.1505 13.858 <2e-16 ***## ---## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## ## Residual standard error: 1.925 on 98 degrees of freedom## Multiple R-squared: 0.6621, Adjusted R-squared: 0.6587 ## F-statistic: 192 on 1 and 98 DF, p-value: < 2.2e-16La potenza del pacchetto nlme è che permette una varietà di funzioni di varianza. La funzione varFixed che abbiamo appena illustrato è la più semplice e qualcosa che può essere fatto con lm. Le altre funzioni di varianza includono:
varIdentvarPowervarExpvarConstPowervarComb
Esaminiamo ognuna di queste usando dati simulati.
varIdent
La funzione varIdent ci permette di modellare diverse varianze per strato. Per vedere come funziona, simuleremo prima i dati con questa proprietà. Di seguito usiamo set.seed(11) nel caso qualcuno voglia simulare gli stessi dati casuali. Poi impostiamo n uguale a 400, il numero di osservazioni che simuleremo. x viene generato come prima. In questo esempio includiamo un predittore aggiuntivo chiamato g che può essere pensato come il genere. Campioniamo casualmente 400 valori di “m” e “f”. Poi simuliamo un vettore di due deviazioni standard e lo salviamo come msd. Se g == "m", allora la deviazione standard è \(2 volte 2,5). Altrimenti è \(2 \) per g == "f". Usiamo questo vettore nella prossima riga per generare y. Notate dove msd è inserito nella funzione rnorm. Notate anche come generiamo y. Includiamo un’interazione tra x e y. Quando g == "f", l’intercetta e la pendenza sono 1,2 e 2,1. Quando g == "m", l’intercetta e la pendenza sono rispettivamente (1,2 + 2,8) e (2,1 + 2,8). Infine mettiamo i nostri dati simulati in un data frame e tracciamo con ggplot2.
set.seed(11)n <- 400x <- seq(1,10, length.out = n)g <- sample(c("m","f"), size = n, replace = TRUE)msd <- ifelse(g=="m", 2*2.5, 2)y <- 1.2 + 2.1*x - 1.5*(g == "m") + 2.8*x*(g == "m") + rnorm(n, sd = msd)d <- data.frame(y, x, g)library(ggplot2)ggplot(d, aes(x, y, color = g)) + geom_point()
Nota le diverse varianze per ogni livello di g. La varianza per “m” è molto più grande della varianza per “f”. Ha molta più dispersione di “f”. Abbiamo simulato i dati in questo modo. Abbiamo impostato la varianza di “m” a 2,5 volte quella di “f”.
Ora usiamo la funzione gls con varIdent per cercare di recuperare questi valori veri. Usiamo lo stesso modo di prima: definiamo la nostra funzione varianza nell’argomento weights. Di seguito specifichiamo nell’argomento form che la formula per modellare la varianza è condizionata a g. L’espressione ~ 1|g è una formula unilaterale che dice che la varianza differisce tra i livelli di g. Il 1 significa solo che non abbiamo covariate aggiuntive nel nostro modello per la varianza. È possibile includere una formula come ~ x|g, ma sarebbe errato in questo caso poiché non abbiamo usato x nel generare la varianza. Guardando il grafico si nota anche che mentre la variabilità in y è diversa tra i gruppi, la variabilità non aumenta all’aumentare di x.
vm2 <- gls(y ~ x*g, data = d, weights = varIdent(form = ~ 1|g))summary(vm2)## Generalized least squares fit by REML## Model: y ~ x * g ## Data: d ## AIC BIC logLik## 2081.752 2105.64 -1034.876## ## Variance function:## Structure: Different standard deviations per stratum## Formula: ~1 | g ## Parameter estimates:## f m ## 1.00000 2.58127 ## ## Coefficients:## Value Std.Error t-value p-value## (Intercept) 0.9686349 0.3237724 2.99172 0.0029## x 2.1222707 0.0525024 40.42239 0.0000## gm -1.9765090 0.9352152 -2.11343 0.0352## x:gm 2.9957974 0.1553551 19.28355 0.0000## ## Correlation: ## (Intr) x gm ## x -0.901 ## gm -0.346 0.312 ## x:gm 0.304 -0.338 -0.907## ## Standardized residuals:## Min Q1 Med Q3 Max ## -2.74115039 -0.67013954 0.01619031 0.69793776 2.72985748 ## ## Residual standard error: 2.005397 ## Degrees of freedom: 400 total; 396 residualPrima si nota in fondo all’output che l’errore standard residuo stimato è 2.0053974, molto vicino al valore “vero” di 2. Notate anche che nella sezione “Funzione varianza” otteniamo un valore stimato di 2,58127 per il gruppo “m”, che è molto vicino al valore “vero” di 2,5 che abbiamo usato per generare la varianza diversa per il gruppo “m”. In generale, quando si usa la funzione varIdent per stimare diverse varianze tra livelli di strati, uno dei livelli sarà impostato come base, e gli altri saranno stimati come multipli dell’errore standard residuo. In questo caso, poiché “f” viene prima di “m” in ordine alfabetico, “f” è stato impostato su baseline, ovvero 1. L’errore standard residuo stimato per il gruppo “f” è \(2,005397 \volte 1\). L’errore standard residuo stimato per il gruppo “m” è \(2,005397 \times 2,58127\).
È importante notare che queste sono solo stime. Per avere un’idea dell’incertezza di queste stime, possiamo usare la funzione intervals del pacchetto nlme per ottenere intervalli di confidenza al 95%. Per ridurre l’output, salviamo il risultato e visualizziamo gli elementi selezionati dell’oggetto.
int <- intervals(vm2)L’elemento varStruct contiene l’intervallo di confidenza al 95% per il parametro stimato nella funzione varianza. Il parametro in questo caso è il moltiplicatore dell’errore standard residuo per il gruppo maschile.
int$varStruct## lower est. upper## m 2.245615 2.58127 2.967095## attr(,"label")## "Variance function:"L’elemento sigma contiene l’intervallo di confidenza al 95% per l’errore standard residuo.
int$sigma## lower est. upper ## 1.818639 2.005397 2.211335 ## attr(,"label")## "Residual standard error:"Entrambi gli intervalli sono abbastanza piccoli e contengono il “vero” valore che abbiamo usato per generare i dati.
varPower
La funzione varPower ci permette di modellare la varianza come potenza del valore assoluto di una covariata. Ancora una volta, per aiutare a spiegare questo, simuleremo i dati con questa proprietà. Di seguito la cosa principale da notare è l’argomento sd della funzione rnorm. È dove prendiamo la deviazione standard di 2 e poi la moltiplichiamo per il valore assoluto di x elevato alla potenza di 1,5. Questo è simile a come abbiamo simulato i dati per dimostrare la funzione varFixed sopra. In quel caso abbiamo semplicemente assunto che l’esponente fosse 0,5. (Ricordate che prendere la radice quadrata di un numero equivale ad elevarlo alla potenza di 0,5). Qui abbiamo arbitrariamente scelto una potenza di 1,5. Quando usiamo gls con varPower tenteremo di recuperare il “vero” valore di 1,5 oltre al 2.
set.seed(4)n <- 1000x <- seq(1,10,length.out = n)y <- 1.2 + 2.1*x + rnorm(n, sd = 2*abs(x)^1.5)d <- data.frame(y, x)ggplot(d, aes(x, y)) + geom_point()
Vediamo i dati che mostrano la classica forma di varianza che aumenta all’aumentare del predittore. Per modellare correttamente questi dati usando gls, gli forniamo la formula y ~x e usiamo l’argomento weights con varPower. Notate che specifichiamo la formula unilaterale proprio come abbiamo fatto con la funzione varFixed. In questo modello, tuttavia, otterremo una stima della potenza.
vm3 <- gls(y ~ x, data = d, weights = varPower(form = ~x))summary(vm3)## Generalized least squares fit by REML## Model: y ~ x ## Data: d ## AIC BIC logLik## 8840.188 8859.811 -4416.094## ## Variance function:## Structure: Power of variance covariate## Formula: ~x ## Parameter estimates:## power ## 1.520915 ## ## Coefficients:## Value Std.Error t-value p-value## (Intercept) 2.321502 0.4750149 4.887218 0## x 1.582272 0.2241367 7.059404 0## ## Correlation: ## (Intr)## x -0.845## ## Standardized residuals:## Min Q1 Med Q3 Max ## -2.892319514 -0.655237190 0.001653178 0.691241690 3.346273816 ## ## Residual standard error: 1.872944 ## Degrees of freedom: 1000 total; 998 residualint <- intervals(vm3)int$varStruct## lower est. upper## power 1.445064 1.520915 1.596765## attr(,"label")## "Variance function:"int$sigma## lower est. upper ## 1.650987 1.872944 2.124740 ## attr(,"label")## "Residual standard error:"La potenza è stimata essere 1,52 che è molto vicina al valore “vero” di 1,5. L’errore standard residuo stimato di 1,8729439 è anche abbastanza vicino a 2. Entrambi gli intervalli catturano i valori che abbiamo usato per simulare i dati. I coefficienti del modello, d’altra parte, sono stimati un po’ male. Questo non è sorprendente dato quanta variabilità c’è in y, specialmente per \(x > 2\).
varExp
La funzione varExp ci permette di modellare la varianza come una funzione esponenziale di una covariata. Ancora una volta spiegheremo questa funzione di varianza usando dati simulati. L’unico cambiamento è nell’argomento sd della funzione rnorm. Abbiamo un valore fisso di 2 che moltiplichiamo per x che è a sua volta moltiplicato per 0,5 ed esponenziato. Notate come la varianza aumenta rapidamente all’aumentare di x. Questo è dovuto alla crescita esponenziale della varianza.
set.seed(55)n <- 400x <- seq(1, 10, length.out = n)y <- 1.2 + 2.1*x + rnorm(n, sd = 2*exp(0.5*x))d <- data.frame(y, x)ggplot(d, aes(x, y)) + geom_point()
Per lavorare a ritroso e recuperare questi valori usiamo la funzione varExp nell’argomento weights di gls. La formula unilaterale non cambia in questo caso. Dice di modellare la varianza in funzione di x. La funzione varExp dice che x è stato moltiplicato per qualche valore ed esponenziato, quindi gls cercherà di stimare quel valore.
vm4 <- gls(y ~ x, data = d, weights = varExp(form = ~x))summary(vm4)## Generalized least squares fit by REML## Model: y ~ x ## Data: d ## AIC BIC logLik## 3873.153 3889.098 -1932.576## ## Variance function:## Structure: Exponential of variance covariate## Formula: ~x ## Parameter estimates:## expon ## 0.4845623 ## ## Coefficients:## Value Std.Error t-value p-value## (Intercept) 1.198525 1.1172841 1.072713 0.284## x 2.073297 0.4933502 4.202486 0.000## ## Correlation: ## (Intr)## x -0.892## ## Standardized residuals:## Min Q1 Med Q3 Max ## -3.16976627 -0.69525696 0.00614222 0.63718022 2.90423217 ## ## Residual standard error: 2.119192 ## Degrees of freedom: 400 total; 398 residualint <- intervals(vm4)int$varStruct## lower est. upper## expon 0.4562871 0.4845623 0.5128374## attr(,"label")## "Variance function:"int$sigma## lower est. upper ## 1.786737 2.119192 2.513505 ## attr(,"label")## "Residual standard error:"La stima “expon” di 0,4845623 nella sezione “funzione varianza” è molto vicina al nostro valore specificato di 0,5. Allo stesso modo l’errore standard residuo stimato di 2,1191918 è vicino al valore “vero” di 2. Anche le stime dei coefficienti del modello sono vicine ai valori che abbiamo usato per generare i dati, ma si noti l’incertezza nell'(intercetta). Il test d’ipotesi nel sommario non può escludere un’intercetta negativa. Questo di nuovo non è sorprendente dato che y ha così tanta varianza non costante così come il fatto che non abbiamo osservazioni di x uguale a 0. Poiché l’intercetta è il valore che y assume quando x è uguale a 0, la nostra intercetta stimata è un’estrapolazione ad un evento che non abbiamo osservato.
varConstPower
La funzione varConstPower ci permette di modellare la varianza come una costante positiva più una potenza del valore assoluto di una covariata. E’ una parola grossa, ma è fondamentalmente la stessa della funzione varPower tranne che ora vi aggiungiamo una costante positiva. Il seguente codice simula i dati per i quali la funzione varConstPower sarebbe adatta da usare. Notate che è identico al codice che abbiamo usato per simulare i dati per la sezione varPower, tranne che aggiungiamo 0,7 a x nella funzione abs. Perché 0,7? Nessuna ragione speciale. È solo una costante positiva che abbiamo scelto.
set.seed(4)n <- 1000x <- seq(1,10,length.out = n)y <- 1.2 + 2.1*x + rnorm(n, sd = 2*abs(0.7 + x)^1.5)d <- data.frame(y, x)ggplot(d, aes(x, y)) + geom_point()
Il modo corretto di modellare questi dati è usare gls con la funzione varConstPower. La formula unilaterale è la stessa di prima. Il sommario restituisce tre stime per il modello della varianza: la costante, la potenza e l’errore standard residuo. Ricordiamo che i valori “veri” sono 0,7, 1,5 e 2, rispettivamente.
vm5 <- gls(y ~ x, data = d, weights = varConstPower(form = ~x))summary(vm5)## Generalized least squares fit by REML## Model: y ~ x ## Data: d ## AIC BIC logLik## 9319.794 9344.323 -4654.897## ## Variance function:## Structure: Constant plus power of variance covariate## Formula: ~x ## Parameter estimates:## const power ## 0.3554921 1.3593016 ## ## Coefficients:## Value Std.Error t-value p-value## (Intercept) 2.763788 0.7607449 3.633003 3e-04## x 1.456677 0.2915621 4.996112 0e+00## ## Correlation: ## (Intr)## x -0.824## ## Standardized residuals:## Min Q1 Med Q3 Max ## -2.846999539 -0.664551194 0.006709051 0.679895874 3.287366683 ## ## Residual standard error: 2.887761 ## Degrees of freedom: 1000 total; 998 residualint <- intervals(vm5)int$varStruct## lower est. upper## const 0.02658804 0.3554921 4.753063## power 1.14258677 1.3593016 1.576016## attr(,"label")## "Variance function:"int$sigma## lower est. upper ## 1.814757 2.887761 4.595195 ## attr(,"label")## "Residual standard error:"Gli intervalli catturano i valori “veri”, sebbene l’intervallo per la costante sia piuttosto ampio. Se non conoscessimo il vero valore, sembrerebbe che la costante potrebbe plausibilmente essere 2 o anche 4.
varComb
Infine, la funzione varComb ci permette di modellare combinazioni di due o più modelli di varianza moltiplicando insieme le funzioni di varianza corrispondenti. Ovviamente questo può ospitare alcuni modelli di varianza molto complessi. Simuleremo alcuni dati di base che sarebbe opportuno utilizzare con la funzione varComb.
Quello che facciamo di seguito è combinare due diversi processi di varianza:
- Uno che permette alla deviazione standard di differire tra i generi (“f” = \(2\), “m” = \(2 \times 2.5\))
- Uno che permette alla deviazione standard di aumentare all’aumentare di
x, dovexè moltiplicato per 1,5 ed esponenziato
Per aiutare la visualizzazione dei dati abbiamo limitato x ad un intervallo da 1 a 3.
set.seed(77)n <- 400x <- seq(1, 3, length.out = n)g <- sample(c("m","f"), size = n, replace = TRUE)msd <- ifelse(g=="m", 2*2.5, 2) * exp(1.5*x)y <- 1.2 + 2.1*x - 1.5*(g == "m") + 2.8*x*(g == "m") + rnorm(n, sd = msd)d <- data.frame(y, x, g)ggplot(d, aes(x, y, color = g)) + geom_point()
Il grafico mostra una varianza crescente all’aumentare di x, ma anche differenze di varianza tra i generi. La varianza di y per il gruppo “m” è molto più grande della varianza di y nel gruppo “f”, specialmente quando x è maggiore di 1,5.
Per modellare correttamente il processo di generazione dei dati che abbiamo specificato sopra e tentare di recuperare i veri valori, usiamo la funzione varComb come involucro attorno ad altre due funzioni di varianza: varIdent e varExp. Perché queste due? Perché abbiamo diverse varianze tra i generi, e abbiamo una varianza che aumenta esponenzialmente in funzione di x.
vm6 <- gls(y ~ x*g, data = d, weights = varComb(varIdent(form = ~ 1|g), varExp(form = ~x)))summary(vm6)## Generalized least squares fit by REML## Model: y ~ x * g ## Data: d ## AIC BIC logLik## 4431.45 4459.32 -2208.725## ## Combination of variance functions: ## Structure: Different standard deviations per stratum## Formula: ~1 | g ## Parameter estimates:## f m ## 1.000000 2.437046 ## Structure: Exponential of variance covariate## Formula: ~x ## Parameter estimates:## expon ## 1.540608 ## ## Coefficients:## Value Std.Error t-value p-value## (Intercept) -5.996337 6.539552 -0.9169339 0.3597## x 8.829971 4.849272 1.8208860 0.0694## gm -17.572238 16.548540 -1.0618603 0.2889## x:gm 16.932938 12.192793 1.3887661 0.1657## ## Correlation: ## (Intr) x gm ## x -0.972 ## gm -0.395 0.384 ## x:gm 0.387 -0.398 -0.974## ## Standardized residuals:## Min Q1 Med Q3 Max ## -2.74441479 -0.67478610 -0.00262221 0.66079254 3.14975288 ## ## Residual standard error: 1.794185 ## Degrees of freedom: 400 total; 396 residualLa sezione riassuntiva contiene due sezioni per modellare la varianza. La prima stima il moltiplicatore per il gruppo “m” a 2,437, che è molto vicino al vero valore di 2,5. Il parametro esponenziale è stimato essere 1,54, estremamente vicino al vero valore di 1,5 che abbiamo usato quando abbiamo generato i dati. Infine, l’errore standard residuo è stimato a circa 1,79, vicino al vero valore di 2. La chiamata intervals(vm6) mostra intervalli di confidenza molto stretti. I prefissi A e B nell’elemento varStruct sono solo etichette per i due diversi modelli di varianza.
int <- intervals(vm6)int$varStruct## lower est. upper## A.m 2.119874 2.437046 2.801673## B.expon 1.419366 1.540608 1.661850## attr(,"label")## "Variance function:"int$sigma## lower est. upper ## 1.380008 1.794185 2.332669 ## attr(,"label")## "Residual standard error:"Purtroppo a causa della grande variabilità esponenziale, le stime dei coefficienti del modello sono pessime.
Che senso ha?
Perciò perché abbiamo simulato tutti questi dati falsi e poi abbiamo cercato di recuperare i valori “veri”? A cosa serve? Domande giuste. La risposta è che ci aiuta a capire quali ipotesi stiamo facendo quando specifichiamo e adattiamo un modello. Quando specifichiamo e adattiamo il seguente modello…
m <- lm(y ~ x1 + x2)… stiamo dicendo che pensiamo che y sia approssimativamente uguale a una somma ponderata di x1 e x2 (più un’intercetta) con gli errori che sono estratti casualmente da una distribuzione normale con media 0 e una deviazione standard fissa. Allo stesso modo, se specifichiamo e adattiamo il seguente modello…
m <- gls(y ~ x1 + x2, data = d, weights = varPower(form = ~x1))…stiamo dicendo che pensiamo che y sia approssimativamente uguale a una somma ponderata di x1 e x2 (più un’intercetta) con gli errori che sono estrazioni casuali da una distribuzione normale con media di 0 e una deviazione standard che cresce come un multiplo di x1 aumentato di qualche potenza.
Se possiamo simulare dati adatti a questi modelli, allora abbiamo una migliore comprensione delle ipotesi che stiamo facendo quando usiamo questi modelli. Speriamo che ora tu abbia una migliore comprensione di ciò che puoi fare per modellare la varianza usando il pacchetto nlme.
Per domande o chiarimenti su questo articolo, contatta l’UVA Library StatLab: [email protected]
Vedi l’intera collezione di articoli dell’UVA Library StatLab.
Clay Ford
Consulente di ricerca statistica
University of Virginia Library
Il 7 aprile 2020
-
1. Il pacchetto
nlmeè forse meglio conosciuto per la sua funzionelme, che è usata per adattare modelli a effetti misti (cioè, modelli con effetti sia fissi che casuali). Questo post sul blog dimostra le funzioni di varianza usandogls, che non adatta gli effetti casuali. Tuttavia tutto ciò che presentiamo in questo post del blog può essere usato conlme.