title: Problems with calculation of numerical identification w.r.t. ANOVA smooth for large scale matrices

text: Suppose we have two (centered) Spline-matrices $\boldsymbol{B_1}$, $\boldsymbol{B_1}$. Then $\boldsymbol{X_1} = [\boldsymbol{B_1},\boldsymbol{B_2}]$ contrains lower order smooths and $\boldsymbol{X_2}$ the corresponding tensor product smooth, ${\boldsymbol{X_2}}_i = {\boldsymbol{B_1}}_i\otimes{\boldsymbol{B_1}}_i$ is the i-th row of ${\boldsymbol{X_2}}$. To check which columns of $\boldsymbol{X_2}$ have to be cancelled in order to remove dependencies between columns of the joined matrix $\boldsymbol{X} = [\boldsymbol{X_1},\boldsymbol{X_2}]$ we first perform a QR decomposition of $\boldsymbol{X_1}$, i.e., $\boldsymbol{X_1} = \boldsymbol{Q_1}\begin{Bmatrix} \boldsymbol{R_1}\\\ \boldsymbol{0} \end{Bmatrix}$ and then perform an additional QR decomposition of $\boldsymbol{B} = \boldsymbol{Q_1}^T\boldsymbol{X_2}$ with the first $r$ rows removed. According to the zero block of $\boldsymbol{R_2}$ from $\boldsymbol{B} = \boldsymbol{Q_2}\begin{Bmatrix} \boldsymbol{R_2}\\\ \boldsymbol{0} \end{Bmatrix}$ ,taking pivoting into account, we know which columns have to be removed **BUT** in order to do that we need to build $\boldsymbol{Q_1}$ which is of size $n^2$!!! I've tried to build $\boldsymbol{Q_1}$ via the householder reflections but after one iteration, $\boldsymbol{Q_1}$ (with $n=150000$) is already 2GB. How can I get the $\boldsymbol{B}$ matrix in order to check for the zero block in $\boldsymbol{R_2}$ if $n$ is quite large? Here the code of building the matrix via householder reflections               qr.Q2 <- function(a){         if(is(a, "sparseQR")){             Q <- Diagonal(a@V@Dim[1])             for(i in rev(1:a@V@Dim[2])){                 x <- Matrix(a@V[ ,i])                 Q <- Q - (x*a@beta[i])%*%(t(x)%*%Q)             }             return(Q[order(a@p), ][1:a@Dim[1], 1:a@Dim[1]])         }else{                 stop("input has to be sparseQR")         }     }


title: Factor analysis of dyadic data

text: An anonymous reader posted the following question on my blog. **Context:** The reader wanted to run a factor analysis on scales from a questionnaire - but the data was from paired husbands and wives. **Question:**   * Can factor analysis be run on dyadic data? If so, how?   * Would the independence assumption hold for factor analysis?


title: Log graph question

text: So let's say you have a distribution where X is the 16% quantile. Then you take the log of all the values of the distribution. Would log(X) still be the 16% quantile in the log distribution?


title: How does the inclusion of an intercept change the variability of the residual?

text: I want to use the variability of the residual as a measure M and then test whether M is higher or lower after some event. However, I estimate separate regression before and after the event to obtain the residuals and the variability of the residuals. Now I am wondering what the effect of including/excluding an intercept in the first regression to obtain the variability of the residuals is? I'm asking because it seems that this has a strong effect on the pre-post results in my design. In formulas: $Y = b_0 + b_1 X + e$ - estimated separately for pre/post; measure $M = \text{var}(e)$ $M = b_0 + b_1 \text{POST} + e$


title: Unit Root Test and rolling data

text: Please could you kindly advise on the implication of testing for unit root (using the augmented dickey fuller approach) on rolling data. My feeling is that this wouldn't make much sense given the non continuous nature of this data. I will be grateful if someone could provide me with some intuition in this regard.


title: What model for a challenging data set? (hundreds of time series with a lot of nesting)

text: I have quite a complicated data set to analyze, and I cant find a good solution for it. Here is the thing: **1.** the raw data is essentially insect song recordings. Each song is made of several bursts, and each burst made of sub-units. All individuals have been recorded for 5 minutes. The number of bursts and their position in the recording can be very different between individuals, as well as the number of sub-units per burst. **2.** I have the carrier frequency (fundamental frequency) of each sub-unit, and that's what I want to analyze. **My problems:** **1.** The frequencies within a burst are not independent obviously (although it's pretty stable, but the frequency of the sub-unit n-1 will have an influence on the sub-unit n). **2.** The bursts are also not independent, within a recording. **3.** They are even less independent as the frequency drops with time (the individual gets tired of singing so the frequency of the song gets lower and lower). The dropping _seems_ to be linear. **4.** Nesting = I have 3 replicated populations for two locations A and B. So I have A1, A2, A3 & B1, B2, B3. **What I would like to do:** **1.** Characterize the difference in frequency between my two locations (test it statistically) **2.** Characterize the frequency dropping between the two locations (see if it drops faster in one of them) **How to do it:** Well that's why I need help: I don't know. It seems that my case combines problems that are usually not seen together. I've read about mixed models, about GAM, about ARIMA, random and fixed effects, but I cant be really sure of the best way to do it. When I graph it though (frequency ~ sub-unit number _n_ ), the difference is very clear between the two locations. I also have to take other variables into account, like the temperature (makes the frequency higher), etc. I thought about:   * Nesting the individuals within the replicate their are from, and nest the replicate within the location (individual/replicate/location).    * Use a random 'burst' effect, so I take into account the variability within each burst.    * Use a fixed 'burst position in recording' effect, to measure the frequency dropping (hoping it is actually linear).  Would it be correct? Is there a special type of model I could use for this kind of scenario?


title: How does Leave-One-Out-Crossvalidation work?

text: I have some data and I want to build a model (say a linear regression model) out of this data. In a next step, I want to apply Leave-One-Out-Cross- Validation (LOOCV) on the model so see how good it performs. If I understood LOOCV right, I build a new model for each of my samples (the test-set) using every sample except this sample (the training-set). Then I use the model to predict the test-set and calulate the errors $(predicted - actual)$. In a next step I aggregate all the errors generated using a choose function, for example MSE, MAPE. I can use these values to judge on the qualitity (or goodness of fit) of the model. Question: Which model is the model these qualitiy-values apply for, so which model should I choose if I find the metrics generated from LOOCV appropriate for my case? LOOCV looked at $n$ models (where $n$ is the sample count), which one is the model I should choose?   * Is it the model which uses all samples? This model was never calcuated during the LOOCV process!   * Is it the model which has the least error? Have I understood something wrong?


title: How to show linear model corresponds to exponential family

text: I am confronted with the exercise below. I have no given solution, so I hope someone can tell me whether my solution is right or wrong. I want to show the deterministic linear model $\quad \boldsymbol y = \boldsymbol X \boldsymbol \beta + \boldsymbol u,\quad \boldsymbol u \sim N(0,\sigma^2 \boldsymbol I) \quad$ corresponds to an exponential model of the following form: $f(y; \boldsymbol \theta) = C(\boldsymbol \theta)h(y) \mathrm{exp} \left\\{\sum\limits_{j=1}^k \xi_j(\boldsymbol \theta)T_j(y) \right\\} \\\$ with $\boldsymbol \theta = (\beta', \sigma^2)'$. Clearly, $\boldsymbol y \sim N(\boldsymbol X \boldsymbol \beta,\sigma^2 \boldsymbol I) $. So I started as follows: $f(y; \boldsymbol \theta)= \prod\limits_{i=1}^n \frac{1}{{\sigma \sqrt {2\pi } }}\mathrm{exp} \left\\{ -\frac{1}{2\sigma^2} (y_i- \sum\limits_{j=1}^kx_{j}^{(i)} \beta_j)^2 \right\\} $ After some rearrangement, I got: $f(y; \boldsymbol \theta)= \frac{1}{{{(2\pi)^{n/2} } }}\mathrm{exp} \left\\{ -\frac{1}{2\sigma^2} \sum\limits_{i=1}^n (\sum\limits_{j=1}^kx_{j}^{(i)}\beta_j)^2 + \mathrm{ln}(\sigma^{n} )\right\\} \\\ \qquad \quad \quad \mathrm{exp} \left\\{ -\frac{1}{2\sigma^2}\sum\limits_{i=1}^n y_{i}^2 + \frac{1}{\sigma^2}\sum\limits_{i=1}^ny_i\sum\limits_{j=1}^kx_{j}^{(i)}\beta_j \right\\} $ So, I concluded it corresponds to the exponential model with the functions $h(y)=\frac{1}{{(2\pi)^{n/2} }} \\\ C(\boldsymbol \theta)=\mathrm{exp} \left\\{ -\frac{1}{2\sigma^2} \sum\limits_{i=1}^n (\sum\limits_{j=1}^kx_{j}^{(i)}\beta_j)^2 + \mathrm{ln}(\sigma^{n} )\right\\}\\\$ EDIT: I had the wrong solution but found the right one with the help of Dennis. I did not bear in mind that T has to be the sufficient statistic. So, I had $\boldsymbol X$ in $\xi$ but could easily change it with matrix notation: $\xi(\boldsymbol \theta) = (-\frac{1}{{2\sigma^2 } },\quad \frac{1}{{\sigma^2 }} \boldsymbol \beta) \\\ T(y) =( \boldsymbol y'\boldsymbol y,\quad \boldsymbol y' \boldsymbol X ) $.


title: Best method for transforming low discrepancy sequence into normal distribution?

text: I've been using low discrepancy sequences for a while for Uniform Distributions, as I've found their properties useful (mainly in computer graphics for their random appearance and their ability to densely cover [0,1] in an incremental fashion). For example, random values above, Halton sequence values below: ![http://i.stack.imgur.com/gCXrh.png](http://i.stack.imgur.com/T3pEr.png) I was considering using them for some financial analysis planning, but I need different distributions than just uniform. I started off trying to generate a normal distribution from my uniform distributions via the Marsaglia polar algorithm, but the results don't seem as good as with the uniform distribution. Another example, again random above, Halton below: ![http://i.stack.imgur.com/NNcdH.png](http://i.stack.imgur.com/7LTRI.png) My question is: What is the best method for getting a normal distribution with the properties I get from a uniform low discrepancy sequence - coverage, incremental fill-in, non-correlation across multiple dimensions? Am I on the right track, or should I be taking a completely different approach? (Python code for uniform and normal distributions I use above: Gist 2566569)


title: How to score Likert scales and test a hypothesis relating two variables?

text: I am want to assess whether successful implementation of newly developed policies are dependent on involvement of all relevant stakeholders and middle managers. I have a 5-point Likert scale with ten questions (items). I collected data on 85 respondents. I assigned values ('Strongly Agree'=1, 'Agree'=2, 'Undecided'=3, 'Disagree'=4, 'Strongly Disagree'=5). I have calculated percents of valid responses, and mean and standard deviation for each item. I have summed up total agree, total disagree and total undecided (neither agree or disagree). I have put these infos into graphs. **What do I do next?**


title: Power calculation in a Bayesian approach

text: Can somebody give me a understanding of what power means in a Bayesian approach? Or how you would implement it in a Bayesian simulation.


title: How to develop a smart sort ordering of list data using popularity/rating algorithms

text: Here is the problem. We have list of Job orders - each job has certain state associated with it, you can think of state as a situation that job is in, for example,   * a job that is said to be done   * or a job that is awaiting someone's approval   * ... We have figured out there is 36 types state a job could go into, also there can be more than one state a job can have at a given time. What we want to achieve is how smartly we can sort the jobs for the user so he know on what job he immediate attention should be. We are also thinking of coloring job names like if it is red - then it is more "hot" job than others or if it is "green" - the logged in user doesn't need to worry about that job and so on. So, the problem is two -   1. How to go about making a job red or green ? Since we know the state of job it can be used to identify whether user's attention required on that job, it is going to be not hard. I mean, if the job falls into any of that 36 types of states at a given time, that goes red. Am I making sense here ?   2. How to sort the list of jobs to show it in the page? What should be the order? Ideally, more hotness (or green or criticality of the attention) the job is the top it should show up. How to do this?. I am afraid we would end of showing a job that requires immediate attention in low order. We can also track metrics such as a) popularity of the job among other users b) oldness or newness of the job etc., however I am not sure whether this will help on this issue. I tried with several algorithms such as Elo, Bayesian but couldn't arrive to a satisfiable solution.


title: What do "single-step" and "multi-step" mean in post-hoc testing of ANOVAs?

text: From Wikipedia > **Post-hoc testing of ANOVAs** > > Multiple comparison procedures are commonly used in an analysis of variance > after obtaining a significant omnibus test result, like the ANOVA F-test. > The significant ANOVA result suggests rejecting the global null hypothesis > H0 that the means are the same across the groups being compared. Multiple > comparison procedures are then used to determine which means differ. In a > one-way ANOVA involving K group means, there are K(K − 1)/2 pairwise > comparisons. > > A number of methods have been proposed for this problem, some of which are: > > **Single-step procedures** > >   * Tukey–Kramer method (Tukey's HSD) (1951) >   * Scheffe method (1953) > > > **Multi-step procedures based on Studentized range statistic** > >   * Duncan's new multiple range test (1955) >   * The Nemenyi test is similar to Tukey's range test in ANOVA. > >   * The Bonferroni–Dunn test allows comparisons, controlling the familywise > error rate.[vague] > >   * Student Newman-Keuls post-hoc analysis >   * Dunnett's test (1955) for comparison of number of treatments to a single > control group. > I was wondering what "single-step" and "multi-step" mean here? Are they both for pairwise comparisons for every two groups, right? Thanks and regards!


title: Consistent estimator with var=0

text: I'm trying to think of an example of a consistent estimator with var=0 but I'm having a hard time. I pretty sure it exists but I'm unable to find one. Any help\hints would be much appreciated.


title: Finding UMVUE of parameter $\displaystyle a^{-\theta}$, $(0<a<1)$

text: Suppose $X_1,X_2,\ldots,X_n$ is a random sample of $\operatorname{Beta}(\theta,1)$. How can I find UMVUE of parameter $\displaystyle a^{\theta}$, $(0<a<1)$?


title: Calibrating GARCH specifications in R

text: We are looking at a data set and try to use garchSpect o define the Garch parameters and garchSim to simulate the garch model. The original data are as shown in the pic below, have a high sd (around 30000) and a high autocorrelation until lag 5 (not seen) When we now simulate with:               aaa<-garchSim(spec = garchSpec(almost_all_explained), n = 365, n.start = 100, extended = FALSE)     plot(aaa)       The result is a Garch process that has a almost no autocorrelation and the produced values lie in the range of +/- 0.0001. How can we adjust / calibrate the deviation and the AR? All other parameters seem to be fine.Thanks! ![enter image description here](http://i.stack.imgur.com/8KMjK.jpg)


title: Calculate the probability that the difference between two populations is >= d

text: I'm a programmer, not a mathematician. I have tried to use real math terms and notation here, but I may be misusing them. I am doing split-testing, and I want a test to compare different click-through rates between two groups, A and B. Existing tutorials and samples suggest a z-test or a t-test to reject the hypothesis that there was no change between A and B, but I have been trying to work through how to compare effect sizes rather than just "are the groups different?" Let's say my two groups are sampled as ($a_k$ successes out of $a_n$ trials) and ($b_k$ successes out of $b_n$ trials). I am trying to determine differences between $A_p$ and $B_p$, the actual click-through probabilities of the two populations. I have learned recently how to use the Beta distribution to model $A_p$ by setting $\alpha = a_k+1$ and $\beta = a_n-a_k+1$. So, it seems to me that I can answer this question: "What is the probability that $B_p >= A_p + P_{delta}$?" by summing the probability of each candidate $A_p$ times the probability that $B_p$ is greater than $A_p+P_{delta}$: $$\int_{p=0}^1 BetaPDF(p,a_k+1,a_n-a_k+1) (1 - BetaCDF(p+P_{delta}, b_k+1,b_n- b_k+1)) $$ Is this line of reasoning correct? It seems like a much more useful question to answer than what the normal t-test examples give you. If this works, it seems like there must be a standard function/distribution/test for doing this? I have written code to do this using a for loop with 1000 steps, but this is very slow when I want to graph over $P_{delta}$ (each datapoint in the graph means iterating the 1000 step for loop). I would appreciate having a better method.


title: Calculate R-squared with JAGS and R

text: I have the following model that I am running in JAGS from R:                model {             for( i in 1:nData){               y[i] ~ dnorm(mu[i], tau)               mu[i] <- b0 + inprod(b[],x[i,])                             }                  tau ~ dgamma(.01,.01)             b0 ~ dnorm(0,.0001)             for (j in 1:nPredictors){               b[j] ~ dgamma(2,2)               }             }      It is a simple regression model with a gamma prior on the beta coefficients. I can run the model and get some reasonable results, but my boss wants to know whether the Bayesian flavor is better than a normal OLS regression, and to this end he wants me to calculate an R^2 statistic for both models. I am fairly new to Bayesian statistics, so if anyone could share some R code to illustrate how it is done, I would be very grateful!


title: What would be a good paper to cite to explain simulation methods in proc mixed in SAS?

text: In order to analyze data obtained from a repeated measure study, I used the following SAS code:               proc mixed;       class id trial_type trial_seq condition;       model decision_quality = condition|trial_seq;       repeated / subject = user_id type=cs;       lsmeans condition / pdiff CL ADJUST=SIMULATE(SEED=18713 NSAMP=200000);     run;      As you can see, I used ADJUST=SIMULATE to find out pair-wise comparison results. What would be a proper citation to explain this method in an academic paper? I think that our reviewers know about a more common correction methods (e.g., Bonferroni and Tukey), but I am not sure if they are aware of the simulation approach.


title: Using ARMA when data is missing

text: I am using ARMA over a dataset with missing samples. How do I treat them? Would you suggest to make linear/nonlinear interpolation or just keep them out and consider two samples with missing data in between as consecutive samples?


title: efficiency - bias trade-off

text: Under which conditions would a researcher choose optimally when there is a trade-off between the variance and bias of an estimator? I hope this question is not too broad... Any help would be appreciated. As an example, what criteria would I use, if we're referring to Linear Models?


title: How to interpret main effects when the interaction effect is not significant?

text: I ran a Generalized Linear Mixed Model in R and included an interaction effect between two predictors. The interaction was not significant, but the main effects (the two predictors) both were. Now many textbook examples tell me that if there is a significant effect of the interaction, the main effects cannot be interpreted. But what if your interaction is not significant? Can I conclude that the two predictors have an effect on the response? Or is it better to run a new model where I leave out the interaction? I prefer not to do so, because I would then have to control for multiple testing.


title: Adequacy of forecasted volatility vs realized volatility. GARCH model

text: Below I have performed forecasting of volatility, in my simple case of 6 days:                        Date      WIG            r    Forecast     61 2014-03-03 51984.88 -0.003157046 0.255080969     62 2014-03-04 54823.49 -0.007566441 0.011048830     63 2014-03-05 50947.86 -0.003180499 0.012665863     64 2014-03-06 49995.55  0.001131508 0.002204160     65 2014-03-07 50988.59  0.005068452 0.047949976     66 2014-03-08 53665.63  0.001263751 0.001458638                r = based on logarithmic return.      My question is how do I check if the forecasted volatility values are adequate and my forecast is "successful". I guess my main concern is how do I calculate realized volatility(actual) so I could compare it to the forecasted volatility and conclude my findings. I have done the research and there is not one concrete formula for realized volatility. Furthemore, I found stats. tests such as MAPE, MSE, MAE but all of them require actual values (realized volatility) in the formula. If anything is not understandable or must be corrected please let me know and I will gladly do it.


title: Estimating event duration by sampling with a Poisson process

text: Take an event to be a set $E \subseteq {\mathbb R}^{\ge 0}$. We want to estimate the duration of the event (the measure of $E$) by sampling using a Poisson process $N$. That is, for some experiment duration $D$, at each time $t\leq D$ such that $\forall t'>t\colon N(t')>N(t)$, we check whether $t\in E$. See TagTime: Stochastic Time Tracking for Space Cadets for more context. Given $n$, the number of sampled times $t\in E$, and $D$, how can we find a $c$-confidence interval for the duration of the event? In the comment thread in the link above, the following formula is proposed, where $g=1/\lambda$ is the mean interval for the Poisson process, and $Q^{-1}$ is the inverse of the regularized incomplete gamma function (`InverseGammaRegularized` in Mathematica): $$\left\\{g Q^{-1}\left(n,\frac{c+1}{2}\right),g Q^{-1}\left(n,\frac{1}{2}-\frac{c}{2}\right)\right\\}$$ I suspect that the $n$ in the upper bound should be $n+1$, because without this, we get no upper bound if $n=0$. Not understanding the math well enough to determine which formula is correct a priori, I wonder how we might go about testing a formula like this by simulation. My guess is that we should run many experiments, where in each experiment, we choose a random event $E$ and use a new Poisson process. If the formula is correct, the fraction of runs in which the event duration is within the bounds should be close to $c$. But how should we choose $E$?


title: Cross-estimating the independent variables to exclude outliers

text: Purpose: pragmatic data mining and prediction, NOT for publication or science Data: observations from nature, so a high degree of stability is expected in the relationships N: approx. 15 k I am working on a dataset where some of the independent observations have measurement issues (i.e. mixed feet and metric). I'm considering cross- estimating all the variables to exclude outliers before building the actual model. In other words, a method for defining and excluding outliers and issues in the data using ALL variables to estimate the independent variables. From the estimate I can then exclude all observations outside a 99 % confidence interval in the independent variables. How statistically viable would this be? Also, would the standard error of the estimate and R square still be meaningful? I have not seen or heard of this before, and I am a bit skeptical. Still, given that this observations in nature, and not social science the confidence interval of the estimate for the independent should be valid. What can make it fall apart is that the very outliers I am trying to get rid of, is included in the method. However, hoping this will be deflated to marginal influence in the method.


title: Comparing test-retest reliabilities

text: Imagine I have a test for a disease that returns either a positive or a negative result. I administer the test to one group twice in a single session; in another group, I administer it twice with two weeks separating the tests. The hypothesis is that the tests taken at the same time will be more consistent with each other than the tests taken two weeks apart. How can I test this? Bonus question: What if the two groups consisted of the same people?


title: Creating a composite variable from 2 log transformed variables?

text: I have a data set of 1,100 observations, I have two continuious variables (counts of different things) that were positively skewed, so I first log transformed each of them. Now I need to make a composite variable, per literature, that measures a construct. My question, per observation/row, should I now add each of the log transformed numbers together OR create a Z score from each log transformed variable and THEN add the z scores together? Finally, because I am testing a moderating hypothesis using the composite variable noted above, with a seperate categorical variable, I know the composite variable has to be centered before a multiplying with the categorical variable. Does this need for a centered variable change the manner I should add them together? Obviously two Z scores are already mean centered vs. adding two log transforms and subtracting the mean. Which way is best?


title: Testing significance of a random effect glmmADMB model

text:
Below is the output from a model of novel object test scores fit with the `nbinom1` (quasi-Poisson) option in `glmmADMB`. I used this package/method because:   1. the Poisson mean is < 5, so according to Bolker et al. 2009 I should not use `glmmPQL`   2. there was overdipersion with the Poisson model run with `glmer`. The NB1 model fit with `glmmADMB` was the best of all fits I tried. There are 3 fixed factors as below, as well as a `Sex`–`Age` interaction term, and individual (`ID`) is the random factor.               Call:     glmmadmb(formula = Score ~ Sex * Age + Object + (1 | ID), data = PGS,          family = "nbinom1")          AIC: 394.5          Coefficients:                   Estimate Std. Error z value Pr(>|z|)     (Intercept)      0.480      0.195    2.47   0.0137     SexM             0.797      0.262    3.05   0.0023     AgeS             1.224      0.254    4.82  1.4e-06     ObjectPLA       -0.434      0.194   -2.24   0.0248     ObjectSCO       -0.822      0.203   -4.05  5.1e-05      SexM:AgeS       -0.825      0.353   -2.33   0.0196       * * *               Number of observations: total=105, ID=40           Random effect variance(s):     Group=ID
                 Variance StdDev
     (Intercept)  0.03738 0.1933          Negative binomial dispersion parameter: 1.5958 (std. err.: 0.34873)          Log-likelihood: -189.228      **Question:** how can I test if the random effect is significant? I have looked around online, and where there are multiple random effects it seems straightforward, but I know I can't compare the `GLMM` to a `GLM` without my single random effect according to info on http://glmm.wikidot.com/faq. I have tried Ben Bolker's `varprof` function in the reef-fish example:               varprof(PGS.nbinom1)     Error in varprof(PGS.nbinom1) :
       trying to get slot "deviance" from an object (class "glmmadmb") that is not an S4 object      As well as his `simulate.glm(PGS.nbinom1)`:               simulate.glm(PGS.nbinom1)     Error in simulate.glm(M1) : family nbinom not implemented      And also Jeff Evans’ way, posted on archiveorange.com, where I ran the model with an uninformative dummy variable (a column with the word "dummy" in every row) and compared that with my model with ANOVA, but also got a warning message here:               anova(M1,M2)     Analysis of Deviance Table          Model 1: Score ~ Sex * Age + Object
     Model 2: Score ~ Sex * Age + Object
       NoPar  LogLik Df Deviance Pr(>Chi)
     1     8 -189.23
                         2     8 -189.32  0    -0.19        1
          Warning message:     In anova.glmmadmb(M1, M2) :
       something's wrong: models should be nested, increasing complexity should imply increasing log-likelihood      What else can I try? Any help would be much appreciated!


title: Calculate principal components

text: Given the following data:               X          X =                   332         428         354        1437         526         247         427              293         559         388        1527         567         239         258              372         767         562        1948         927         235         433              406         563         341        1507         544         324         407              386         608         396        1501         558         319         363              438         843         689        2345        1148         243         341              534         660         367        1620         638         414         407              460         699         484        1856         762         400         416              385         789         621        2366        1149         304         282              655         776         423        1848         759         495         486              584         995         548        2056         893         518         319              515        1097         887        2630        1167         561         284          [m,n]=size(X)          m =              12          n =               7      I have calculated centered data and covariances from this matrix:               mean_matrix=X-repmat(mean(X),size(X,1),1);          covariance=(mean_matrix'*mean_matrix)/(12-1);      The dimensions of this covariance matrix are:               [m1,n1]=size(covariance)          m1 =               7          n1 =               7       After that I have calculated the eigenvalue decomposition:               >> [V,D]=eig(covariance);     >> [e,i]=sort(diag(D),'descend');     >> sorted=V(:,i);     >> sorted          sorted =             -0.0728    0.5758    0.4040   -0.1140   -0.1687   -0.6737   -0.0678        -0.3281    0.4093   -0.2917   -0.6077    0.4265    0.1828    0.2348        -0.3026   -0.1001   -0.3402    0.3965    0.5682   -0.4320   -0.3406        -0.7532   -0.1082    0.0681    0.2942   -0.2848    0.0011    0.4987        -0.4653   -0.2439    0.3809   -0.3299   -0.0645    0.2076   -0.6503        -0.0911    0.6316   -0.2254    0.4135   -0.2366    0.4390   -0.3498         0.0588    0.1444    0.6599    0.3068    0.5705    0.3005    0.1741          >> e          e =             1.0e+05 *              2.7483         0.2642         0.0625         0.0230         0.0209         0.0034         0.0007      Now I want to know about projection: is it $V'*X$ or $X*V$?


title: how to draw dependent binomial sample

text: Hello I would like to draw samples from two binomial distributions, say B1 ~ binomial(100,p1) and B2 ~ binomial(100,p2) and I want them to be dependent. Any suggestions on how to do it? I am using R.


title: Treating numerals/cardinals in Bag of Words (BOW) model

text: I wish to do topic modeling on text corpus some of which are about company earnings which has lots of numbers in it. It has no sentence structure. I think tagging numbers using nltk.pos_tagging can help me find out if the number is CD (numeral/cardinal). Using numeral feature as one of the many features in the BOW can help me identify the topic related to Statistics or Maths. Intuitively if I drop this information about numbers appearing in the text, I will not be able to figure out the topic of unseen texts. I will really appreciate your views regarding this. Also, if anyone knows about a BOW library in python that can do so internally that will be helpful. I came across scikit-learn library's count vectorizer in which I think if I change token_pattern, it might help me in getting numbers as features but that will not help me in collapsing all numeric features as one feature. I hope I am clear. Really appreciate your time and help. Thanks.


title: How do I calculate the p-value in a gene enrichment analysis, using parametric and non-parametric methods?

text: I have a list of 20 QTL from an experiment involving a certain phenotype related to bone development. Of all the candidate genes within the 20 QTL intervals, 5 are associated with a certain disease in humans. I tested the set of 20 QTL to see if they are enriched for the disease-associating genes. Each of the 20 QTL has a certain length in base pairs. I took that set of 20 lengths, and randomly placed them in the mouse genome. I forced them not to overlap and to land no closer than 500,000 bases on either side. Then I took note of the genes inside each interval and counted how many matched the names contained in the disease-associating list of genes. I repeated the process 10,000 times. This provides the expectations given randomness... how many times can I expect to see disease-associating genes inside a set of 20 randomly placed QTL? I compare that to my actual observed number. Then I calculate the probability. To do that, I simply counted how many times 1 disease-associating gene appeared within each set of 20 randomly placed QTL and divide by 10,000. That gives me the frequency. I did that for 2 disease-associating genes, and 3, and so on up to 5, which is my observed value. I find that the frequency of 5's in the random data is less than 0.05. Is this appropriate? Also, I see online some might suggest a parametric method, like Fishers' exact test, but that is supposed to be used for very small data sets, am I correct? Am I missing something there? How can I implement a parametric calculation for p?


title: Nested, Stratified Ratios

text: I need some advice. I have a stratified sampling design by season and area. However, within any stratum the statistic I'm interested in $\hat{R} = \sum y/\sum x$ is also affected by other factors. To overcome this, I basically take a weighted average across both factors to obtain a single estimate $\hat{R}$ for each stratum. The only problem is that the two factors are nested within one another, such that I need to take a weighted average over two levels. This looks something like this: $ \hat{R} = \frac{1}{W} \sum_{a}\frac{w_{a}}{N_{a}} \sum_{b} \frac{y n_{a,b}}{x} $ where $n$ and $w$ are the weights as fractions. I just need to know how to calculate the variance (if possible), without having to use a bootstrap.


title: Multinomial logistic analysis

text: I am estimating a logistic model for a binary $Y$ (0=controls, 1=cases), and a set of covariates ($x_1, x_2, x_3, \ldots$) including sex (male, female). I know from the data that risk of $y$ is quite different for males and females with certain covariates. What I was trying to do, with no success so far, was to fit one single model for entire data to get estimates for $x_1, x_2, x_3$ for example, and gender specific OR for let's say $x_4, x_5$. I am totally confused with terminology: I have tried nested analysis (which does not apply here because I have one observation per subject), multilevel modeling, and multinomial (0=controls, 1=male cases, 2=female cases). Multinomial seems like what I need, but in Stata that I am using, `mlogit` gives me coefficients for y=1 and y=2 for all covariates, and I can not figure out how to get aggregate coefficients for all male and female cases for certain set of covariates (which are not supposed to be gender-specific). There must be some straightforward way to get OR for both from the same model, without manually calculating the coefficients and CIs. Any suggestions would help.


title: Partitioning the appropriate error term in a repeated measures R

text: Basically I have a repeated measures experiment with two reward levels and two speed levels. Every trial participant can be randomly presented with highReward*slowSpeed or lowReward*fastSpeed etc., resulting in a 2 x 2 design. I am not completely sure when it is appropriate to use either of the models below:               model <- lme(value ~ reward*speed, random = (~reward+speed|participant),data)      adding the participant:reward:speed error term to the model               model <- lme(value ~ reward*speed, random = (~reward*speed|participant),data)      When it is appropriate to add the participant:reward:speed error term to the model?


title: MAD in relation to 95% confidence

text: MAD (Median Absolute Deviation) is: $\text{MAD} = M_i(|x_i-M_j(x_j)|)$ where $M()$ is the median operator ($M_i(x_i) = \text{median}(x_1,...,x_n)$). I'd like to scale the MAD in such a way as to include (say) 95% of a distribution around the median, the way that that 95% of a normal distribution is within $1.96\sigma$ of the mean. That is, if $m = M_i(x_i)$ and $d = \text{MAD}_i(x_i)$, make an interval like $m \pm b\cdot d$ (where $b$ depends on the distribution you are dealing with) that includes 95% of the distribution. Can this be done?


title: Outlier detection in out-sample data for the purpose of classification

text: This is the most widely used method for outlier detection in econometrics and statistical problems. `X` is our data that we're searching for outliers in it (in MATLAB) :                abs(X-mean(X)) >= n*std(X)      So if this inequality was true, that sample is an outlier; otherwise we will keep the sample. I'm using `neural network` and `SVM` for my classification problem. Before normalization of data, I find outlier for each feature in input database. After creating my `MLP` or `SVM` model now I want use it for out-sample data ( new data from **this year** \- model trained by data of **last year** ).   * Should I use outlier detection in this new database?   * If answer to that question is yes, Should I use `mean` and `std` of previous year that was used for outlier detection of main data, or use new `mean` and `std` of out-sample data?


title: Which property of count data make mean-variance dependency?

text: I have read about the fact that, there is dependency of variance on mean of count data.In most of cases they do variance stabilization transfomration as preprocessing step of data modeling. I wonder, why variance depends on mean in count based data ? in other word, which properties of count data make this happen ?


title: questions regarding the order statistics of a normal distribution

text: Assume $x_1, x_2, ..., x_n$ are realizations from a normal distribution with mean $0$ and standard deviation $1$ and assume they are ordered from small to large. Let $y_i=x_i-x_{i-1}$ be the difference of the two adjacent numbers (so $y_i$ is always larger than $0$). My questions are: 1) What is the distribution of $y_i$? and 2) How can I test if $y_i$ is statistically different from $0$?


title: Standard error of mean for a distribution with two dependent variables

text: I want to calculate the mean and standard error of the mean for the amount spent per person visiting a store. Most people don't buy anything, so the distribution looks like this ${\rm P}({\rm purchased},{\rm spend}) = {\rm P}({\rm purchased}) {\space} {\rm P}({\rm spend} | {\rm purchased})$ where purchased and spend are two dependent variables -- dependent in that if purchased = False, then spend is necessarily 0. One way to do this would be by taking the mean and standard error from the collection of amount spent by each person:   $ T = [\$13, \$6, \$0, \$0, \$0, \$0, \$0, \$0, \$0, \$5, \$0, \$0, \$0, ...]$ $ E(T) = \frac{1}{N} \sum_{i=1}^{N} T_{i}$ $ {\rm standard \space error} = \frac{\sqrt{\frac{1}{N-1} \sum_{i=1}^{N} (T_{i} - E(T))^2}}{\sqrt{N}} $ However, I'm worried about missing the inherent binomial distribution, given that some people bought something and others didn't. Is there a more correct way to calculate the standard error of the mean?


title: How to describe a distribution of statistics in a meta-analysis?

text: I'm conducting a meta-analysis using water-quality data collected from multiple lakes in multiple regions. Due to the differences in the sampling programs, each waterbody has been sampled an unequal number of times. I would like to describe the distribution of concentrations at the regional level (rather than the waterbody level). The initial options I came up with are: 1) Create an index for each waterbody (i.e., the mean) and determine the distribution of those means (e.g., describe the grand mean and 95th percentile of the means) 2) Calculate a weighted mean for the region (using the number of samples as the weights). For option 1, I'm finding it difficult to determine what the 95th percentile of a distribution of means really describes? This would also apply to other estimates of dispersal. For option 2, am I limited to calculating a weighted mean only? Can you calculate a weighted percentile? Thoughts on other options and graphical techniques to display these data? Thanks!


title: Accuracy of random prediction with non equal distribution

text: Assume that I want to predict the value of a variable that has three different states: a, b, and c. The chance that these variables have the 3 states is not equally distributed. Out of 10 trials, the distribution is like this: a: 2 b: 3 c: 5 If I have a random predictor that randomly predicts whether the variable is a, b, or c, the chance are intuitively not 1/3. But what are the chances and how can I calculate them? In addition, I have a follow up question. Assume these new distribution: a: 1 b: 0 c: 9 Here one can see, that if my classifier just always predicts c, he is in 90% of the times correct. Is there any way to combine these two concepts? I want to calculate the accuracy of a classifier that predicts randomly, but only if it's not better to just predict always the same.


title: Does dimension reduction always lose some information?

text: Like the title says, does dimension reduction always lose some information? Consider for example PCA. If the data I have is very sparse, I would assume a "better encoding" could be found (is this somehow related to the rank of the data?), and nothing would be lost.


title: Testing for a root larger than 1 in AR(p)

text: My question is the following: How to test for a root larger than 1 in AR(p) process from its observations. Thanks in advance.


title: Bayesian variable selection

text: Chapter 13 of Kevin Murphy's book Machine Learning: A Probabilistic Perspective discusses Sparse Linear Models. After a short introduction on the benefits of sparse models, he introduces the following problem: ![enter image description here](http://i.stack.imgur.com/AciWF.png) How does he derive equation 13.1 above? i.e. why does it take that form, and what is $f$ supposed to represent here?


title: learn R (Statistics)

text: I'm considering using R for my master thesis in order to do some analyses with the mvpart package. The problem is that I only have limited time to learn the program and to do the analyses. I already read through some information concerning R and I'm starting to get the feeling that it is impossible to get the know the program with the limited time I have. I also have very limited background concerning statistics and statistical programs. I only did some limited work on SPSS and SAS. I have no background in programming. Is anyone here familiar with R? What are your thoughts?


title: Why is Mantel's test preferred over Moran's I?

text: **Mantel's test is widely used in biological studies** to examine the correlation between the spatial distribution of animals (position in space) with, for example, their genetic relatedness, rate of aggression or some other attribute. Plenty of good journals are using it ( _PNAS, Animal Behaviour, Molecular Ecology..._ ). I fabricated some patterns which may occur in nature, but Mantel's test **seems to be quite useless to detect them.** On the other hand, Moran's I had better results _(see p-values under each plot)_. > **Why don't scientists use Moran's I instead? Is there some hidden reason I > do not see? And if there is some reason, how can I know (how the hypotheses > must be constructed differently) to appropriately use Mantel's or Moran's I > test? A real-life example will be helpful.** _**Imagine this situation:_** There is an orchard (17 x 17 trees) with a crow is sitting on each tree. Levels of "noise" for each crow are available and you are want to know if the spatial distribution of crows is determined by noise they make. _**There are (at least) 5 possibilities:_**   1. **"Birds of a feather flock together."** The more similar crows are, the smaller the geographical distance between them **(single cluster)**.   2. **"Birds of a feather flock together."** Again, the more similar crows are, the smaller the geographical distance between them, **(multiple clusters)** but one cluster of noisy crows has no knowledge about the existence of second cluster (otherwise they would fuse into one big cluster).   3. **"Monotonic trend."**   4. **"Opposites attract."** Similar crows cannot stand each other.   5. **"Random pattern."** The level of noise has no significant effect on spatial distribution. For each case, I created a plot of points and used the Mantel test to compute a correlation (it is no surprise that its results are non-significant, I would never try to find linear association among such patterns of points). ![enter image description here](http://i.stack.imgur.com/TWQqa.png) * * * **Example data:** _(compressed as possible)_               r.gen   <- seq(-100,100,5)     r.val   <- sample(r.gen, 289, replace=TRUE)     z10     <- rep(0, times=10)     z11     <- rep(0, times=11)     r5      <- c(5,15,25,15,5)     r71     <- c(5,20,40,50,40,20,5)     r72     <- c(15,40,60,75,60,40,15)     r73     <- c(25,50,75,100,75,50,25)     rbPal   <- colorRampPalette(c("blue","red"))     my.data <- data.frame(x = rep(1:17, times=17),y = rep(1:17, each=17),                  c1=c(rep(0,times=155),r5,z11,r71,z10,r72,z10,r73,z10,r72,z10,r71,                  z11,r5,rep(0, times=27)),c2 = c(rep(0,times=19),r5,z11,r71,z10,r72,                  z10,r73,z10,r72,z10,r71,z11,r5,rep(0, times=29),r5,z11,r71,z10,r72,                  z10,r73,z10,r72,z10,r71,z11,r5,rep(0, times=27)),c3 = c(seq(20,100,5),                  seq(15,95,5),seq(10,90,5),seq(5,85,5),seq(0,80,5),seq(-5,75,5),                  seq(-10,70,5),seq(-15,65,5),seq(-20,60,5),seq(-25,55,5),seq(-30,50,5),                  seq(-35,45,5),seq(-40,40,5),seq(-45,35,5),seq(-50,30,5),seq(-55,25,5),                  seq(-60,20,5)),c4 = rep(c(0,100), length=289),c5 = sample(r.gen, 289,                   replace=TRUE))          # adding colors     my.data$Col1 <- rbPal(10)[as.numeric(cut(my.data$c1,breaks = 10))]     my.data$Col2 <- rbPal(10)[as.numeric(cut(my.data$c2,breaks = 10))]     my.data$Col3 <- rbPal(10)[as.numeric(cut(my.data$c3,breaks = 10))]     my.data$Col4 <- rbPal(10)[as.numeric(cut(my.data$c4,breaks = 10))]     my.data$Col5 <- rbPal(10)[as.numeric(cut(my.data$c5,breaks = 10))]      Creating matrix of geographical distances (for Moran's I is inversed):               point.dists           <- dist(cbind(my.data$x, my.data$y))     point.dists.inv       <- 1/point.dists     point.dists.inv       <- as.matrix(point.dists.inv)     diag(point.dists.inv) <- 0      Plot creation:               X11(width=12, height=6)     par(mfrow=c(2,5))     par(mar=c(1,1,1,1))          library(ape)     for (i in 3:7) {       my.res <- mantel.test(as.matrix(dist(my.data[ ,i])), as.matrix(point.dists))       plot(my.data$x,my.data$y,pch=20,col=my.data[ ,c(i+5)], cex=2.5, xlab="",             ylab="", xaxt="n", yaxt="n", ylim=c(-4.5,17))       text(4.5, -2.25, paste("Mantel's test", "\n z.stat =", round(my.res$z.stat,             2), "\n p.value =", round(my.res$p, 3)))            my.res <- Moran.I(my.data[ ,i], point.dists.inv)       text(12.5, -2.25, paste("Moran's I", "\n observed =", round(my.res$observed,             3), "\n expected =",round(my.res$expected,3), "\n std.dev =",             round(my.res$sd,3), "\n p.value =", round(my.res$p.value, 3)))     }          par(mar=c(5,4,4,2)+0.1)          for (i in 3:7) {       plot(dist(my.data[ ,i]), point.dists,pch = 20, xlab="geographical distance",             ylab="behavioural distance")     }      P.S. in the examples on UCLA's statistics help website, both tests are used on the exact same data and the exact same hypothesis, which is **not very helpful** (cf., Mantel test, Moran's I). **Response to I.M.** You have write: > ...it [Mantel]tests whether quiet crows are located near other quiet crows, > while noisy crows have noisy neighbors. I think that **such hypothesis could NOT be tested by Mantel test**. On both plots the hypothesis valid. But if you suppose that one cluster of not noisy crows may not have knowledge about the existence of second cluster of not noisy crows - Mantels test is again useless. Such separation should be very probable in nature (mainly when you are doing data collection on larger scale). ![enter image description here](http://i.stack.imgur.com/w4gav.png)


title: A question on continuous random variable

text: Let say, I have 2 continuous random variables X1 & X2. Both have same location parameters. Other parameters may be same or may not. Now say, the q1-th quantile of X1 is less than the q1-th quantile of x2. But the q2-th quantile of x1 is more than the q2th quantile of x2. My question is, is that possible? Is there any example of x1 & x2 which have that property? I will be really grateful if someone can give me some pointer. * * * Edit: At this point, I realize the question I asked for was not correctly specified. I'm particularly interested in the case where the two quantiles being considered are on the same side of the location parameter.


What is the difference between sample variance and sampling variance? They seem same. Aren't they?


title: Which of the parameters in LibSVM is the slack variable?

text: I am a bit confused about the namings in the SVM. I am using this library LibSVM. There are so many parameters that can be set. Does anyone know which of these is the slack variable? thx


title: Measuring inter-subject agreement in user ratings

text: In my system I have users rating objects (e.g. films) from **1 to 5 stars**. The mean rating clearly tells something about the overall opinion on an object, but I'd like to use a more precise measure of variability. For example, if the stars are set randomly, if half of the users give 1 and the other half is 5, the average is still 3, so it doesn't tell me the agreement among users. I'd like a measure that tells me the **inter-subject agreement** , ranging from 1=all users agree, 0=random ratings. What would be a good and standard measure for this aspect of a set of ratings?


title: Network graph with avatars as nodes

text: I've got a problem. I need to create network graph where avatars are nodes (different sizes). Similar to this one: ![enter image description here](http://i.stack.imgur.com/d6o6k.jpg) Furthermore on this image every avatar has the same size, what I need to do is different size of avatar depending on a chosen statistic (e.g. in-degree). I cannot find any tutorial in R, which explains how to do that.


title: Cohen's d for 2x2 interaction

text: I have a significant 2x2 interaction in a between-subjects design. I would normally report omega-squared as my effect-size measure, but I've been asked to provide Cohen's d instead? Can Cohen's d be computed for an interaction? If so, how?


title: What is the most accurate way of determining an object's color?

text: I have written a computer program that can detect coins in a static image (.jpeg, .png, etc.) using some standard techniques for computer vision (Gaussian Blur, thresholding, Hough-Transform etc.). Using the ratios of the coins picked up from a given image, I can establish with good certainty which coins are which. However, I wish to add to my confidence levels and also determine if a coin that I deduce to be of type-A (from radius ratios) is also of the correct colo[u]r. The problem is that for British coins et al. (copper, silver, gold), the respective colors (esp. of copper to gold) are very similar. I have a routine that extracts the mean color of a given coin in terms of the RedGreenBlue (RGB) 'color-space' and routines to convert this color into HueSaturationBrightness (HSB or HSV) 'color-space'. RGB is not very nice to work with in attempting to differentiating between the three coin colors (see attached [basic] image for an example). I have the following ranges and typical values for the colours of the different coin types: _Note: the typical value here is one selected using a 'pixel-wise' mean of a real image._               **Copper RGB/HSB:** typicalRGB = (153, 117, 89)/(26, 0.42, 0.60).          **Silver RGB/HSB:** typicalRGB = (174, 176, 180)/(220, 0.03, 0.71).          **Gold RGB/HSB:** typicalRGB = (220, 205, 160)/(45, 0.27, 0.86)       I first tried to use the 'Euclidian distance' between a given mean coin color (using RGB) and the typical values for each coin type given above treating the RGB values as a vector; for copper we would have: $$D_{copper} = \sqrt((R_{type} - R_{copper})^{2} + (G_{type} - G_{copper})^{2} + (B_{type} - B_{copper})^{2})$$ where the smallest value of the difference ($D$) would tell us which type the given coin is most likely to be. This method has shown itself to be very inaccurate. I have also tried just comparing the hue of the coins with the typical values of the types provided above. Although theoretically this provides a much better 'color-space' to deal with varying brightness and saturation levels of the images, it too was not accurate enough. **Question:** What is the best method to determine a coins type based on color (from a static image)? Thanks very much for your time. ![Typical Coin Colors](http://i.stack.imgur.com/GpZHm.jpg) **Edit 1** Note: I have tried all of the ideas discussed below and have achieved next to nothing. Variance in lighting conditions (even within the same image) make this problem very tough and should be taken into consideration. **Edit 2 (Summery of Outcome)** Thank you for your answers. Further research of my own (including your answers and comments) has highlighted just how tough this problem is to deal with in the generic case of arbitrary lighting, arbitrary camera (mobile device), fluctuation in coin colour (even for same species/type) etc. I first looked at skin colour recognition (a very active field of research) as a starting point and there are still numerous problems even with the recognition of skin colour for Caucasians alone (see this paper for a review of the current techniques), and the fact that this problem contains three distinct colour objects all of which can have continuous and varying chromacities make this topic of computer vision a very hard one to classify and deal with accordingly (in fact you could do a good Ph.D. on it!). I looked into the Gamut Constraint Method from the very helpful post by D.W. below. This was at first sight very promising as a pre-processing step to transform the image and the separate coin objects to colours that are independent of lighting conditions. However, even this technique does not work perfectly (and involves a library of images/histograms for mappings – which I don’t want to get into) and neither does the much more complex Neural Network Architecture methodologies. In fact this paper states in the abstract that:               "current machine colour constancy algorithms are not good enough for colour-based       object recognition.".      That is not to say that there aren’t much more up-to-date papers on this subject out there, but I can't find them and it does not seem to be a very active research area at this time. The answer by AVB was also helpful and I have looked into L*A*B* briefly.               "The nonlinear relations for L*, a*, and b* are intended to mimic the nonlinear     response of the eye. Furthermore, uniform changes of components in the L*a*b* colour     space aim to correspond to uniform changes in perceived colour, so the relative      perceptual differences between any two colours in L*a*b* can be approximated by      treating each colour as a point in a three dimensional space."      From what I have read, the transformation to this colour space for my device dependent images will be tricky - but I will look into this in detail (with a view to some sort of implementation) when I have a bit more time. I am not holding my breath for a concrete solution to this problem and after the attempt with L*A*B* I shall be neglecting coin colour and looking to sure- up my current geometric detection algorithms (accurate Elliptic Hough Transform etc.). Thanks you all. And as a end note to this question, here is the same image with a new geometric detection algorithm, which has no colour recognition: ![enter image description here](http://i.stack.imgur.com/G4Pnp.jpg)


title: Measure accuracy of Holt-Winters model

text: I'm really confused about measuring the accuracy of Holt-Winters fitted models applying different transformations. How do i compare the accuracy between models when i apply no transformation to the data, the BoxCox power transformation and, for instance, log transformation? Say i have a time series x.               fit1 <- HoltWinters(x)     f1 <- forecast.HoltWinters(fit1, h=12)     accuracy(f1)          fit2 <- HoltWinters(BoxCox(x))     f2 <- forecast.HoltWinters(fit2, h=12, lambda=BoxCox.lambda(x))     accuracy(f2)          fit3 <- HoltWinters(log(x))     f3 <- forecast.HoltWinters(fit3, h=12)     accuracy(f3)      The first 2, i think are comparable because i passed the parameter lambda to the forecast.HoltWinters() to the second one. But what about the other one? Do i apply the inverse if the natural logarithm like this               exp(accuracy(f3))      Any help is appreciated. Thanks guys.


title: Multinomial logistic regression assumptions

text: What are the proper assumptions of Multinomial Logistic Regression? And what are the best tests to satisfy these assumptions using SPSS 18?


title: Analyzing change in categorical data with age

text: I have data in an excel spreadsheet of around 500 participants in a taste experiment. The data includes the age of the participant, in addition to one of five locations on their tongue in which they tasted five separate flavours, e.g. (Sweet, Middle), (Salt, Back) - The image below shows the first 20 entries to demonstrate the formatting. ![http://i.stack.imgur.com/FsqYW.png](http://i.stack.imgur.com/FsqYW.png) I'd preferably like to be able to use R to see if taste location is dependent upon age, finding the point at which it becomes dependent (if this is indeed the case). Any general methods of doing this would be great.


title: Binary classification and ROC curve area less than .5

text: Both my independent and dependent variables are binary. My result for classification table is 72% for predicted, and my ROC curve area is 0.389. Since <0.5 for ROC area is the worst for accuracy model, what should I do?


title: Prediction interval for a fitted log-normal distribution

text: What I am trying to do is to fit a log-normal distribution to a data-set, and then determine confidence and prediction intervals for the fitted distribution - not just for the mean and sd estimates. My final goal is to be able to say that if we repeated a set of measurements, then 95 % of the values would fall below some specific value. I can obviously do that from the fitted distribution itself (or rather the cumulative version), but I'm thinking that will under-estimate the probability because I need to include the variability in the estimates and/or fit itself first? If I use the following code:               require(MASS)     set.seed(123)     x<-rlnorm(100)     fit<-fitdistr(x,"lognormal")      then R will calculate a log-normal distribution fitted to my data. The fitdistr function will return the estimated mean and sd (along with standard errors for these estimates).                  meanlog       sdlog         0.09040591   0.90824033      (0.09082403) (0.06422229)      I understand that these will then allow me to plot the fitted distribution (and histogram) using ggplot2 with the following code:               meanlog<-fit$estimate[[1]]         sdlog<-fit$estimate[[2]]     binwidth<-abs(max(x)-min(x))/20          qplot(x,geom="blank")+geom_histogram(binwidth=binwidth,aes(y= ..density..))+stat_function(fun=dlnorm,arg=list(meanlog=meanlog,sdlog=sdlog),colour="red")      However, what I really want to do is to plot the confidence interval and/or prediction interval of this fitted distribution. Something similar to how ggplot2 does with stat_smooth, something like:               x<-seq(1,100)     y<-x+rnorm(x,sd=10)     qplot(x,y,geom="point")+stat_smooth(method='lm',se=T)      I can use confint(fit) to extract the confidence intervals for the estimated mean and sd, but I think I misunderstand the maths because I can't for the life of me work out how to use those in order to be able to calculate the confidence interval for the actual distribution. So neither can I work out the prediction interval. I've tried writing my own function for the log-normal distribution to input various combinations of the confidence intervals from confint manually - but that doesn't work. Obviously a confidence interval of the estimate does not directly give you the confidence interval of the line. And, therefore, neither the prediction interval. I would really appreciate anyone who can walk me through this please!


title: Low explained variance in Random Forest (R randomForest)

text: I am using randomForest in R for regression, I have many categorical predictors (all of them have the same 3 categories (0,1,2)) and I want to see which of them can predict the response (continuous). I am trying this with many different response variables (one at the time) and all the models have a very low explained variance (basically 0, almost always negative). I checked chi-square between pairs of variables and removed the ones that could be associated (p-value < 0.05), but the result is the same. My questions are: 1 - Is this possible? Am I doing something very wrong without noticing? If no: 2 - In random forest, do I have to throw everything away or can I still use the variable importance for classifying the predictors? (I don't think so, but since I couldn't find anything about this, I still hope I can get something out of it - BTW why does the plot predicted vs observed look good??). If no: 3 - Any suggestion? Also for alternative methods? In the example below I don't divide the data in training and test for simplicity, but I did it in my code - same problem. Also, my original data set is much bigger (>500 observations and almost 100 predictors)               ## predictors     > pred        X1 X2 X3 X4 X5 X6 X7 X8 X9 X10 X11 X12 X13 X14 X15 X16 X17 X18 X19 X20     1   0  0  0  0  0  0  0  0  0   0   0   0   0   0   0   0   0   0   0   0     2   0  1  2  2  2  0  1  2  0   0   1   0   0   1   1   2   2   1   1   2     3   0  1  0  2  2  1  1  2  1   1   2   1   0   0   1   2   2   2   0   0     4   0  0  1  1  1  1  1  1  1   1   1   1   0   0   2   0   2   2   0   1     5   0  1  1  2  2  0  1  2  2   1   2   0   0   0   1   1   0   2   0   1     6   1  1  0  2  2  1  1  1  2   1   0   1   0   1   1   2   2   2   1   2     7   0  1  1  1  1  1  2  1  2   1   2   1   0   1   1   2   1   2   1   1     8   0  1  2  1  0  1  0  2  1   1   1   2   0   0   1   2   1   2   1   2     ........          ## response     > resp          [1]  19.416  46.058  39.496  79.752 301.012 746.377 277.721  13.922  15.598  82.195  86.263     [12]  82.522  30.829 101.369  31.496  39.366 133.510          ## find optimal value of mtry for randomForest     > bestmtry <- tuneRF(pred, resp, ntreeTry=100,     +                    stepFactor=1.5,improve=0.01, trace=F, plot=F, dobest=FALSE)          ## extract optimal value of mtry for randomForest     > ind <- as.numeric(names(which.min(bestmtry[,2])))          ## Random Forest     > RF <-randomForest(pred, , y = resp, mtry=ind, ntree=500,     + keep.forest=TRUE, importance=TRUE)          > RF          Call:          randomForest(x = pred, y = resp, ntree = 500, mtry = ind, importance = TRUE,          keep.forest = TRUE)                     Type of random forest: regression                          Number of trees: 500     No. of variables tried at each split: 9                    Mean of squared residuals: 32713.86                         % Var explained: -6.5          ## Low explained variance (pseudo - r sqaured)          > RF.pr = predict(RF,pred)          ## the plot isn't that bad though... it is if I use the test data set though     > plot(RF.pr, resp)     > abline(c(0,1),col=2)      ![enter image description here](http://i.stack.imgur.com/78z9Z.png)               > varImpPlot(RF)      ![enter image description here](http://i.stack.imgur.com/1MHBJ.png) I have been stuck with this for a while now... any help is extremely appreciated


title: Verifying independence for two Random Variables

text: Could you please give me some hints for the exercise below? Suppose we toss a coin once and let $p$ be the probability of heads. Let $X$ denote the number of heads and let $Y$ denote the number of tails. I have to first show that   * $X$ and $Y$ are dependent    * and afterwards if we let $N\sim $Poisson $(\lambda)$ and toss the coin $N$ times, that $X$ and $Y$ are independent. * * * The first part follows from the fact that two events are disjoint if we toss the coin just once. For the second part, I can see that the vector $(X,Y)$ follows the multinomial distribution with paramemeters $N$, $p$, and $(1-p)$. Therefore: $$P\left( X=x,Y=y,N=n \right)=\binom{n}{x\ \ y } p^x (1-p)^y \times \frac { e^{-\lambda}\lambda^n}{n!}$$ for $x+y=n$ and $x,y \geq 0$ I understand that the marginal distributions are binomial but I do not immediately see how I could proceed. I would appreciate some advice here. Thank you.


title: convert zip codes to latitude and longitude

text: I have a data set that includes a list of zip codes. Is there a way in R to convert these zip codes to latitude and longitude coordinates? Thanks!


title: How to generate survival data with time dependent covariates using R

text: I want to generate survival time from a Cox proportional hazards model that contains time dependent covariate. The model is $h(t|X_i) =h_0(t) \exp(\gamma X_i + \alpha m_{i}(t))$ where $X_i$ is generated from Binomial(1,0.5) and $m_{i}(t)=\beta_0 + \beta_1 X_{i} + \beta_2 X_{i} t$. The true parameter values are used as $\gamma = 1.5, \beta_0 = 0, \beta_1 = -1, \beta_2 = -1.5, h_0(t) = 1$ For time-independent covariate (i.e. $h(t|X_i) =h_0(t) \exp(\gamma X_i) $i generated as follows               #For time independent case     # h_0(t) = 1     gamma <- -1     u <- runif(n=100,min=0,max=1)     Xi <- rbinom(n=100,size=1,prob=0.5)     T <- -log(u)/exp(gamma*Xi)      Can anyone please help me to generate survival data with time-varying covariate.


title: Item Analysis for an R newbie

text: I am trying to assess a 20-item multliple choice test. I want to perform an item analysis such as can be found in this example. So for each question I want the P-value and the correlation with the total, and the distribution of the options selected. I don't know anything about the various statistical software packages out there, but I'd like to use R as I'm comfortable with programming and R is open source. The pseudo-workflow I envision is:   1. prepare data in excel and export to CSV    2. load data in R   3. load a package that does what I need   4. execute that package's commands    5. export and report. I am confident with 1 and 2 but having trouble with 3, probably because I don't have the statistical vocabulary to compare the packages I browsed on CRAN. `ltm` looks like it could be the right package but I can't tell. Whatever package is used, what would the commands be? Side question: in the linked example, what do you suppose MC and MI stand for?


title: When are confidence intervals useful?

text: If I understand correctly a confidence interval of a parameter is an interval constructed by a _method_ which yields intervals containing the true value for a specified proportion of samples. So the 'confidence' is about the method rather than the interval I compute from a particular sample. As a user of statistics I have always felt cheated by this since the space of all samples is hypothetical. All I have is one sample and I want to know what that sample tells me about a parameter. Is this judgement wrong? Are there ways of looking at confidence intervals, at least in some circumstances, which would be meaningful to users of statistics? [This question arises from second thoughts after dissing confidence intervals in a math.se answer http://math.stackexchange.com/questions/7564/calculating- a-sample-size-based-on-a-confidence-level/7572#7572 ]


title: Using Singular Value Decomposition to Compute Variance Covariance Matrix from linear regression model

text: I have a design matrix of p regressors, n observations, and I am trying to compute the sample variance-covariance matrix of the parameters. I am trying to directly calculate it using svd. I am using R, when I take svd of the design matrix, I get three components: a matrix $U$ which is $n \times p$, a matrix $D$ which is $1\times 3$ (presumably eigenvalues), and a matrix $V$ which is $3\times 3$. I diagonalized $D$, making it a $3\times 3$ matrix with 0's in the off- diagonals. Supposedly, the formula for covariance is: $V D^2 V'$, however, the matrix does not match, nor is it even close to R's built in function, `vcov`. Does anyone have any advice/references? I admit that I am a bit unskilled in this area.


title: Simulation study: how to choose the number of iterations?

text: I would like to generate data with "Model 1" and fit them with "Model 2". The underlying idea is to investigate robustness properties of "Model 2". I am particularly interested in the coverage rate of the 95% confidence interval (based on the normal approximation).   * How do I set the number of iteration runs?    * Is it true that larger than necessary replications may result in spurious biases? If so, how is that?


title: Difference between ridge regression implementation in R and SAS

text: I have been reading the description of ridge regression in _Applied Linear Statistical Models_ , 5th Ed chapter 11. The ridge regression is done on body fat data available here. The textbook matches the output in SAS, where the back transformed coefficients are given in the fitted model as:   $$ Y=-7.3978+0.5553X_1+0.3681X_2-0.1917X_3 $$ This is shown from SAS as:               proc reg data = ch7tab1a outest = temp outstb noprint;       model y = x1-x3 / ridge = 0.02;     run;     quit;     proc print data = temp;       where _ridge_ = 0.02 and y = -1;       var y intercept x1 x2 x3;     run;     Obs     Y    Intercept       X1         X2         X3           2     -1     -7.40343    0.55535    0.36814    -0.19163      3     -1      0.00000    0.54633    0.37740    -0.13687      But R gives very different coefficients:               data <- read.table("http://www.cst.cmich.edu/users/lee1c/spss/V16_materials/DataSets_v16/BodyFat-TxtFormat.txt",                         sep=" ", header=FALSE)     data <- data[,c(1,3,5,7)]     colnames(data)<-c("x1","x2","x3","y")     ridge<-lm.ridge(y ~ ., data, lambda=0.02)        ridge$coef     coef(ridge)          >   ridge$coef            x1        x2        x3      10.126984 -4.682273 -3.527010      >   coef(ridge)                        x1         x2         x3      42.2181995  2.0683914 -0.9177207 -0.9921824      >       Can anyone help me understand why?


title: Function to find the quantile in a vector corresponding to constant $x$

text: Suppose I have constant `x=0.1` in the language `R` and I have a vector `vec = rnorm(200)`. Is there a pre-packaged function to find the quantile of `vec` that corresponds the closest to `x`? A solution is as follows:               x = 0.1     vec = rnorm(100)     percentiles = quantile(vec,seq(0,1,by=0.01))     which(abs(x-percentiles)==min(abs(x-percentiles)))      #returns closest match      ... but I would like to know if there's a pre-packaged function.


title: Properties of Battacharyya distance vs Kullback-Leibler divergence

text: What properties do these measures have and how can I determine which one is better for a given purpose? What are extreme cases where they differ a lot?


title: What is currently the more accepted Agreement Between Raters Statistics or most referenced paper?

text: I found a question that is related here, but it doesn't really goes on what I want to know. I found a couple of papers using Kappa Statistic from 2006, and 2010, but afterwards I found other authors disagreeing with it like here and here. I also noticed other tests being listed and suggested like the ICC on the question pointed on cross validated I pointed, and on the wiki. To narrow down the context: I have a set n of participants which consistently were asked to say 'yes/no/maybe' on a set of questions. I've seen authors also changing the yes/maybe to 'evidence' and converting it to yes/no by just summing the quantities of yes and maybe. My questions thus are:   1. What is the most recent/accepted or referenced paper in respect to this discussion?   2. Given my dataset constraints, which test would be more appropriated? I was considering Fleiss but when I saw it was a generalization of the pi-statistic and that the Kappa came out of criticism to pi i got confused. If you could suggest an R function that could do it I would also appreciate. Thank you.


title: How to ´check if two (ergodic) distributions statistically differ (in R)

text: I need to check if two ergodic distributions statistically differ. I know about the Kolmogorov - Smirnov test. However, being the ergodics computed as in Feyrer (2003), _Convergence by parts_ , I do not have two vectors of observations on which performing the test. Is there any way to do this ? Moreover, is there any available function for R ?


title: Reversing Chebyshev inequality argument

text: One way one could state Chebyshev's inequality is > The probability that a realization deviates from the mean more than $k$ > standard deviations is at most $\frac{1}{k^2}$. My question is: Can one rigorously reverse this logic and make a statement about the probability that the actual mean is close to the observation. One immediate technical problem is that one needs to define a probability space on possible probability distributions/means. I'm asking because I think this type of argument (although slightly more convoluted) underlies Vapnik-Chervonenkis theory. In their textbooks this issue is not discussed at all. They prove a large deviation principle and then simply invert all their inequalities. How does this work? Does it?


title: Predicting future survival probability with time dependent covariates

text: I have developed an attrition model for clients (proc phreg) with time dependent covariates in the observation window. How can I predict the future survival probabilities in this case?


title: Are standardized betas in multiple linear regression partial correlations?

text: Since standardized betas are correlation coefficients in bivariate regression, is it the case that standardized betas in multiple regression are partial correlations?


title: Is a chisquare on a (nearly) complete population data necessary?

text: I have data from a survey that was (attempted to to be) administered to all children of particular grades in a certain state. I am getting it after a cleaning step by the survey designers that removed obvious invalid answers (from obnoxious teens). Questions A and B have binary answers, and I'm interested in reporting the percentages of children in the 2x2 categories. Of the approximately 100,000 observations I have, there are 3,500 that have missing data for either A or B, and are not included in the table. There's decent reason to believe that for these 2 questions, non-responses won't be particularly biased one way or the other. What is the proper way to test/summarize any differences between the categories? Is a chi-square test meaningful here? If the non-responses are unbiased, do I just have a really large random sample? Or can I assume that I am just reporting the actual proportions and no statistical testing is needed?


title: What do you call the left hand side variables in a regression/classification?

text: If I have a regression or classification problem               y ~ x1 + x2 + ... + xn      I might call the `xi` the _regressors_ or the _predictors_. But what do you call `y`? The terms _regressand_ and _predictand_ are clunky, don't seem particularly standard, and have the disadvantage that they're very similar to the words for the right-hand-side variables (I'm going to be using these terms in type declarations and documentation for a program, so minimizing the potential for confusion is a priority). Also, using _regressor_ and _regressand_ is inappropriate when you have a classification problem, rather than a regression. Is there a commonly used term?


title: Is bootstrapping a valid method to assess the uncertainty of the median estimate?

text: Bootstrapping works well to access the uncertainty in the mean estimate, however I remember reading somewhere the bootstrap does not do a good job in assessing the uncertainty in quantile estimates (particularly the median). I don't remember where I read this, and I couldn't find much with a quick Google search. Thoughts on this and any references would be greatly appreciated.


title: How to tell if residuals are autocorrelated from a graphic

text: When you do an OLS regression and plot the resulting residuals, how can you tell if the residuals are autocorrelated? I know there are tests for this (Durbin, Breusch-Godfrey) but I was wondering if you can just look at a plot to gauge if autocorrelation could be a problem (because for heteroskedasticity it is fairly easy to do so).


title: Confusion regarding least squares method

text: I am having some confusion regarding least squares method. Actually, least squares method is for minimizing the square of the $L_2$ norm of $Ax-b$ as given in this video lecture. However I am confused if we are given a set of points like $(x_1,y_1),(x_2,y_2)...(x_4,y_4)$. Then our model is something like $$y=ax+b$$ where we have to find the parameters and b. So we need to minimize the square of the $L_2$ norm of $$(y-ax-b)$$ isn't it?


title: Why is my shared variance negative?

text: I have two questions regarding standard multiple regression:   1. Why is my shared variance a negative number?   2. Should I only include the positive semipartial correlations when calculating uniquely explained variance?  I am trying to calculate the amount of shared variance explained in a regression model with four predictor variables, and this number is coming out negative (-.465). According to Tabachnick & Fidell, (2001), uniquely explained variance is computed by adding up the squared semipartial correlations. Shared variance is computed by subtracting the uniquely explained variance from the R square. My R square value = .325 (R = .570, Adj R square = .295, fwiw) F(4, 91) = 10.941, p = .0005. The semipartial correlation values are (significant predictors indicated by*, from the ‘Part’ column in SPSS output): .172* -.174* .465* .164 I calculated that the predictors collectively uniquely explained 80% of the variance (0.801, which is the sum of POSITIVE semipartial correlation coefficients). When calculating the shared variance, the figure comes out at -.48 (computed by subtracting the uniquely explained variance from the R square value; .32 - .80 = -.48). I’m not sure whether it is possible to have a negative value for shared variance, where have I gone wrong (if I have)? Any advice would be greatly received. Many thanks.


title: Help with the Information Extraction plugin in RapidMiner

text: Im trying to apply text mining on a text document with the help of Information Extraction plugin. I dont know which operators to use and how. Please help how do I proceed further? Any good video tutorial available for this plugin? I couldn't find any. ![enter image description here](http://i.stack.imgur.com/YFrhW.png)


title: Predictive posterior in distance dependent Chinese Restaraunt Process

text: I'm reading the paper "Distance Dependent Chinese Restaurant Processes" by Blei and Frazier. On page 8 predictive posterior for the process is described for a new data point $x_{new}$: $$ p(x_{new} | \mathbf{x}) = \sum_{c_{new}} p(c_{new}) \sum_{\mathbf{c}} p(x_{new} | c_{new}, \mathbf{c}, \mathbf{x}) p(\mathbf{c} | \mathbf{x}) $$ Now consider that we have observed data $\mathbf{x}$ and two new points $x_1$ and $x_2$. Using chain rule I can write predictive posterior for these two points as (with parameters omitted): $$ p(x_1, x_2 | \mathbf{x}) = p(x_2 | x_1, \mathbf{x}) p(x_1 | \mathbf{x}) $$ $p(x_1 | \mathbf{x})$ could be computed according to expression (8), but when I try to derive $p(x_2 | x_1, \mathbf{x})$ in a similar way updated posteior $p(c_1, \mathbf{c} | x_1, \mathbf{x})$ appears. It is a bit conter-intuitive because updated posterior could change considerably and thus has to be re- computed. I wonder if anyone could explain me is there a mistake in my reasoning or posterior of already observed data must be re-computed after each new observed data point.


title: Non-stationary series keep close to each other but correlation between growth rates is ~0 - how is this possible?

text: I have 2 (monthly) time-series that look like this: ![cross plot of 2 series](http://i.stack.imgur.com/oewFd.png) Economical intuition suggests that they are positively related and I can see this on the plot but if I compute correlation between their log-returns $\ln x_t/x_{t-1}$ and $\ln y_t/y_{t-1}$ this correlation is -0.04 this is basically zero and not statistically significant for my data size (~60 points). How can it be? One may say that series are cointegrated $y_t = a x_t +\varepsilon_t$, but then returns should follow $\Delta y_t = a \Delta x_t +\Delta \varepsilon_t$ and correlation between returns would also be significant. So if I see zero correlation between returns, there is no cointegration between levels as well - right? So does this zero correlation means that there is no relation between series? If yes - why do they follow each other so closely.. If no - how to quantify this relation if correlation between diff'ed series is ~0 and cointegration tests for original series are inconclusive. EDITS: \-- added cointegration -> correlation link to address @AlecosPapadopoulos question.


title: Are MRFs with edges to all observed data possible?

text: I have been discussing the following issue with a colleague of mine and I can't seem to wrap my head around it. I have a computer vision background, so I'm mostly familiar with 2D MRFs/CRFs for image restoration and segmentation. What is clear to me is that MRFs typically have a simple 2D grid structure for p(x) and each underlying variable x_i has an edge to the corresponding pixel y_i in the graph for the posterior p(y|x). In CRFs you have an edge to each pixel y in the image. I can fully understand that CRFs are trained in a discriminative way and thus they are tractable and that it wouldn't be possible to really learn a generative model based on such a CRF graph. However in all the literature it seems that by definition the MRF model cannot have edges from a x_i to all pixels y, as it then would no longer be a MRF. For example the following states: > If the conditional independence assumption is not used, the posterior will > usually not be a MRF making the inference difficult. Now that it would not be tractable to actually learn such a model I can understand, but why would it no longer be an MRF and also what would it be then? Because I don't see why it would by definition become a CRF, as we do not NEED to condition it on the data. I really hope someone can tell me what I am missing here, or where my mistake is!


title: Resources for learning how to implement ensemble methods

text: I understand theoretically (sort of) how they would work, but am not sure how to go about actually making use an ensemble method (such as voting, weighted mixtures, etc.).   * What are good resources for implementing ensemble methods?   * Are there any particular resources regarding implementation in Python? EDIT: To clear up some based on the discussion on the comments, I'm not looking for ensemble algorithms such as randomForest, etc. Instead, I'm wondering how can you combine different classifications from different algorithms. For example, say someone uses logistic regression, SVM, and some other methods to predict the class of a certain observation. What is the best way to go about capturing the best estimate of the class based upon these predictions?


title: Estimator bias without a closed form?

text: Given a regression loss function $l(Z,\beta)=||Y-Z\beta||_2 + \lambda \beta^TD\beta + r(X,Z)$ where $X$ is the predictor matrix, I would like to estimate a $Z$ that minimizes the above loss in a regression setting. $D$ is a fixed matrix similar to the setting in ridge regression, but $r(.)$ is non- convex. How can I estimate the bias of the minimizer $\hat{\beta}$, especially because, I only have iterative updates that converge to a minimum, and do not have a closed form for both $\hat{\beta}, \hat{Z}$? On the whole my iterative updates have this structure: $Z_n=f(\beta_{n-1})$ $\beta_n=g(Z_{n-1},\beta_{n-1})$ How do I compute $\mathbb{E}(\beta_n)$ as $n \to \infty$? This confuses me! But I can vouch for $||\beta_n-\beta_{n-1}||$ converging to $0$ for some, $n \in \mathbb{Z}^+$. Apart from the question of computing the bias or approximating an estimate for the bias, what are the theoretical/empirical strategies I can use to measure the quality of my estimator?


title: How can you convert Bayesian Information Criterion parameters to a probabilistic interpretation?

text: I'm working with a general bayesian information criteria to determine if this data exhibits a bimodal, trimodal, quadmodal etc. My existing BIC exhibits clear trimodality, but I'd like a hypothesis test of sorts. How can I determine the probability that this model is true?


title: Test a significant difference between two slope values

text: The data I have are a regression slope value of y~time, a standard error, an n value and a p value, for a particular species in two different areas. I want to check whether the the regression slope for one area is significantly different from the regression slope for the other area - is this possible with such data? Does anyone have any suggestions how I could go about this? I can't get access to the raw data unfortunately... Sorry that this is such a simple question!


title: Copy variable value labels

text: I have two datasets, or .dta files. One has certain variables where each value has a corresponding label. The other datafile has no value labels, but the values of the variables are similar to certain variables in the first datafile. What I want is to copy those value labels from the first file to the second file.


title: How to do cross-validation with a Cox proportional hazards model?

text: Suppose I have constructed a prediction model for the occurrence of a particular disease in one dataset (the model building dataset) and now want to check how well the model works in a new dataset (the validation dataset). For a model built with logistic regression, I would calculate the predicted probability for each person in the validation dataset based on the model coefficients obtained from the model building dataset and then, after dichotomizing those probabilities at some cutoff value, I can construct a 2x2 table that allows me to calculate the true positive rate (sensitivity) and the true negative rate (specificity). Moreover, I can construct the entire ROC curve by varying the cutoff and then obtain the AUC for the ROC graph. Now suppose that I actually have survival data. So, I used a Cox proportional hazards model in the model building dataset and now want to check how well the model works in the validation dataset. Since the baseline risk is not a parametric function in Cox models, I do not see how I can get the predicted survival probability for each person in the validation dataset based on the model coefficients obtained in the model building dataset. So, how can I go about checking how well the model works in the validation dataset? Are there established methods for doing this? And if yes, are they implemented in any software? Thanks in advance for any suggestions!


title: Degrees of Freedom Calculation GLMER R

text: I am trying to find out the degrees of freedom in a glmer model in R. I know that their calculation can be tricky, but I can't find any package or clear explanation on how to obtain them. My model is: glmer(Response~A+B+C+(1|D)+(1|E), family="poisson") Where A is a factorial vairable with three levels and B, C, D, E are continuous variables. Thank you very much!


title: Dependent thinning Poisson process

text: If $N_1$ and $N_2$ are independent Poisson processes then the superposition is a Poisson process. Is it possible to construct two dependent Poisson processes such that the superposition is a Poisson process again? How can I do that? Can I use some copula directly? An article that suggests this is possible is Pfeifer and Neslehova (2004) - Modelling and generating dependent risk processes for IRM and DPA. How to do it exactly however is unclear to me. A related question to this topic was asked in: Mixing and dividing point processes


title: Zero-truncated Sichel distribution in R

text: I'm trying to fit the GAMLSS library's Sichel distribution to some zero- truncated data, but the only way to get the function to work is to include the zero-class anyway but give it a frequency of 0, which doesn't take into account the zero-truncated nature of my data. Can anyone suggest a way to properly "redistribute" the zero-class's probability to the remaining probabilities (or some other, better, course of action using Sichel)? If you run the following example, you'll see that `sum(pdf2)` equals 1, but that the zero class that has a probability in my case of 0 is still allocated around 27% of the cum probability:               Counts = data.frame(n = c(0,1,2,3,4,5,6,7,8,9,10),                         freq = c(0,182479,76986,44859,24315,49,100,490,106,0,2))          gamlss(n~1,family=SICHEL, control=gamlss.control(n.cyc=50),data=Counts )          pdf2 = dSICHEL(x=with(Counts, n), mu = 1.610, sigma = 98.43, nu = 3.315)          print( with(Counts, cbind(n, freq, fitted=pdf2*sum(freq))), dig=9)          sum(pdf2)


title: Testing the normality assumption for repeated measures anova? (in R)

text: So assuming that there is a point in testing the normality assumption for anova (see 1 and 2) How can it be tested in R? I would expect to do something like:               ## From Venables and Ripley (2002) p.165.     utils::data(npk, package="MASS")     npk.aovE <- aov(yield ~  N*P*K + Error(block), npk)     residuals(npk.aovE)     qqnorm(residuals(npk.aov))      Which doesn't work, since "residuals" don't have a method (nor predict, for that matter) for the case of repeated measures anova. So what should be done in this case? Can the residuals simply be extracted from the same fit model without the Error term? I am not familiar enough with the literature to know if this is valid or not, thanks in advance for any suggestion.


title: Year fixed effects

text: I am having some confusion with year fixed effects. If I would like to estimate the effect of Policy A on Outcome Y, will the year fixed effects capture the effect of Policy B (also occurring in this time period and affecting Y, but not specifically controlled for in the equation) on Y?


title: Representative sample or not?

text: Please help me out to determine of those are representative samples, not, or either: First example: A conservative talk radio host conducts a poll among his listeners. He asks his listeners to phone in whether they support the President’s policies or not. Will this poll be a representative sample of the nation’s population? Second example: A liberal commentator on a television show asks his listeners to phone in whether or not they support the President’s policies. Ninety percent answer yes. Is this likely to be a representative sample of the U.S. population?


title: CV on training set with feature selection

text: I've got a problem with CV on feature selection. I've used a method, but I don't know it's correct...   1. I split my data into 70% training set and 30% test set    2. I work now with my training set. I do on my training set a 10-fold CV    3. On each fold, I use the training part to search for the minimal subset of features that maximizes the accuracy of the test part (in that fold)   4. Now I've got 10 subsets of features and I rank all the features from this subsets and check which minimal subset maximizes the accuracy of the 30% test set I really don't know if this is correct, can someone help me?


title: Negative binomial regression

text: In a specification where negative binomial regressions are estimated, can the coefficients be interpreted directly as in the OLS case?


title: Bootstrap and var of median

text: I've 2 different questions: **A.** Sample with n variables, while $n-1$ variables are in the range (0-1), and one of the variable is very high: $10^{1000}$. If I'll use bootstrap to estimate the mean, what bias will I get? **B.** Please help me to find the variance of the median (let $n=3$). The distribution is unknown. Please explain the both ways:   1. finding the distribution and then the variance   2. using bootstrap Thanks a lot!


title: A motion model to track a moving object using the condensation algorithm

text: I have implemented the condensation algorithm in order to track a moving object in video sequences, however the predictive step does not work properly, so the samples moves excessively compared to the actual movement of the object. Maybe I need to estimate more accurately the velocity of the object: how to estimate the two components of velocity?


title: How to code values for males when they're only meaningful for females?

text: I'm building a predictive model for a medical condition that happens in both men and women. Physicians have reported that in some cases for women, their menstrual history seems to be a risk factor, though epidemiologists haven't been able to verify this claim. My data includes information about female subjects' menstrual history. My goal is to build a logistic regression model that determines the probability of the existence of this condition. Here's my problem: I'd like to build a model that includes gender and includes menstrual history questions for female subjects, but how do I assign values for male subjects for things like "number of months since last pregnancy"? Is there a general method to deal with this kind of problem?


title: categorizing monthly gross income

text: what is the best way to categorize or even bin monthly gross income? As of now, the variable for income is numeric and continuous. I have seen many examples for annual income. Should I categorize variables into quartiles? The purpose would be to fit the transformed variable into a logistic regression. As of now, I do not feel comfortable with interpreting the continuous variable, as it seems strange to say for every dollar increase? Thanks!


title: Scatter plots of p-value distribution of gene expression data in R

text: How do I distinctly represent two data sets on the same scatter plot , to compare the p-value distributions of the genes using R packages? The gene expression values and p-values (obtained by ANOVA) of the two datasets are:               ### data A              probes     AG1      AG2      LR1       LR2           901  4.991648  6.289011 5.556497 5.698004           902  5.465736  6.036498 7.126263 7.004082           903  8.559740  9.833907 8.567032 8.433505           904  7.068298  6.156696 6.436085 6.793267           905  5.388158  4.911885 5.642148 5.914592          ### data A (p-value)         901           902          903          904          905      8.705207e-01 2.084426e-05 2.313093e-03 2.027818e-01 1.273762e-01                ### data B                  probes  AG1    AG2      LR1     LR2              901 4.562965 4.917767 4.636715 5.223545              902 4.662215 4.845352 4.893514 5.139439              903 5.206248 4.396288 4.274906 5.458326              904 6.027329 6.169694 6.350698 6.317571              905 4.981953 4.980007 4.948762 5.118272          ### data B(p-value)             901           902          903          904          905      3.225798e-01 2.205045e-01 3.957011e-01 7.262735e-01 9.253321e-01       I have around 15000 genes to visualize and would like to produce a scatter plot for comparing the p-value distribution in both datasets.


title: Does fastICA require PCA to run at first?

text: I reviewed an application based paper saying that applying PCA before applying ICA (using fastICA package). My question is that does ICA (fastICA) requires PCA to run at first? This paper mentioned that "it is also argued that pre-applying PCA enhances ICA performance by (1) discarding small trailing eigenvalues before whitening and (2) reducing computational complexity by minimizing pair-wise dependencies. PCA decorrelates the input data; the remaining higher-order dependencies are separated by ICA." Also other papers applying PCA before ICA, e.g., Paper 2. Are there any other pros and cons to run PCA before ICA? Please provide theory with references. Thanks a lot. A.


title: Data classification methods beside LDA

text: I tried both Principal Component Analysis (PCA) and Linear Discriminant Analysis (LDA) in order to classify me data. While LDA is much better than PCA, the results are not good enough. Are there any other alternatives that I can use? Thanks