title: Different estimated parameters in similar models in R

text: A particular series (std), seems to exhibit a trend-like behavior. According to the ADF test for this series:               Dickey-Fuller = -2.8618, Lag order = 6, p-value = 0.2131      Therefore, I am taking the first difference of std with this code               stddif1<- diff(std)      Here is the tricky part, the acf and pacf suggest that this would be an ARMA process (2,1), with a d=1. But the code shows both different estimates and different AIC values, when (I think) this _shouldn't_ be the case: For std with no difference:               > arima(std, order=c(2,1,1))          Call:     arima(x = std, order = c(2, 1, 1))          Coefficients:          ar1     ar2      ma1       0.5206  0.2697  -0.7638     s.e.  0.1218  0.0552   0.1153          sigma^2 estimated as 0.06355:  log likelihood = -13.3,  aic = 34.6      And, for the differenced std (stddif):               > arima(stddif, order=c(2,0,1))          Call:     arima(x = stddif, order = c(2, 0, 1))          Coefficients:          ar1     ar2      ma1  intercept       0.5188  0.2695  -0.7620    -0.0003     s.e.  0.1223  0.0554   0.1159     0.0157          sigma^2 estimated as 0.06355:  log likelihood = -13.3,  aic = 36.6      The values for the AR1, AR2, MA1 as well as the AIC are different. Why is this? This was all done in R, the relevant package is 'tseries'.


title: Deriving Linear MMSE Estimator

text: I am trying to verify a certain derivation of **linear minimum mean square error estimator** as it appears in [1] (also called _linear a posteriori estimate_ , since the estimate is based on measurement). Linear estimator is assumed $\hat{\mathbf{x}}(\mathbf{z}) = \mathbf{Az} + \mathbf{b}$, where $\mathbf{A} \in \mathbb{R}^{n\times p},\ \mathbf{b} \in \mathbb{R}^{n}$ are constant and $\mathbf{z} \in \mathbb{R}^{p}$ is measurement of $\mathbf{x} \in \mathbb{R}^n$ (both $\mathbf{x}$ and $\mathbf{z}$ are random variables). The goal is to find such $\mathbf{A},\,\mathbf{b}$ that minimize $J = \overline{(\mathbf{x} - \hat{\mathbf{x}})^T(\mathbf{x} - \hat{\mathbf{x}})}$. (Overbar denotes expectation.) The difficulty I am facing is in simplifying this expression to the form presented in the book. The derivation in [1] proceeds as follows: $$ J = \overline{(\mathbf{x} - \hat{\mathbf{x}})^T(\mathbf{x} - \hat{\mathbf{x}})} = \mathrm{Tr}\ \overline{(\mathbf{x} - \hat{\mathbf{x}})^T(\mathbf{x} - \hat{\mathbf{x}})} = \mathrm{Tr}\ \overline{(\mathbf{x} - \mathbf{Az} - \mathbf{b})(\mathbf{x} - \mathbf{Az} - \mathbf{b})^T} = \mathrm{Tr}\ \overline{\left[(\mathbf{x} - \bar{\mathbf{x}}) - (\mathbf{Az} + \mathbf{b} - \bar{\mathbf{x}})\right]\left[(\mathbf{x} - \bar{\mathbf{x}}) - (\mathbf{Az} + \mathbf{b} - \bar{\mathbf{x}})\right]^T} $$ I the last equality the ''smart zero'' $ -\bar{\mathbf{x}} + \bar{\mathbf{x}}$ was added. Here I am still able to follow. Authors then write (and that's where the problem lies) ''After some straighforward work we arrive at'' $$ J = \mathrm{Tr}\ \left[\mathbf{P}_x + \mathbf{A}(\mathbf{P}_z + \overline{\mathbf{zz^T}})\mathbf{A}^T + (\mathbf{b} - \bar{\mathbf{x}})(\mathbf{b} - \bar{\mathbf{x}})^T + 2\mathbf{A}\bar{\mathbf{z}}(\mathbf{b} - \bar{\mathbf{x}})^T - 2\mathbf{AP}_{zx} \right] $$ where $$\mathbf{P}_x = \overline{(\mathbf{x} - \bar{\mathbf{x}})(\mathbf{x} - \bar{\mathbf{x}})^T}, \ \mathbf{P}_{zx} = \overline{(\mathbf{z} - \bar{\mathbf{z}})(\mathbf{x} - \bar{\mathbf{x}})^T} $$ (analogously for $\mathbf{P}_z$ and $\mathbf{P}_{xz} = \mathbf{P}^T_{zx}$) are covariance matrices. * * * **My work so far**   If I multiply out the inner expression under overbar in $$\mathrm{Tr}\ \overline{\left[(\mathbf{x} - \bar{\mathbf{x}}) - (\mathbf{Az} + (\mathbf{b} - \bar{\mathbf{x}}))\right]\left[(\mathbf{x} - \bar{\mathbf{x}}) - (\mathbf{Az} + (\mathbf{b} - \bar{\mathbf{x}}))\right]^T} $$   I get $$ \mathrm{Tr}\ \left[\overline{(\mathbf{x} - \mathbf{\bar{x}})(\mathbf{x} - \mathbf{\bar{x}})^T} - 2\overline{(\mathbf{x} - \mathbf{\bar{x}})(\mathbf{Az} + (\mathbf{b} - \bar{\mathbf{x}}))^T} + \overline{(\mathbf{Az} + (\mathbf{b} - \bar{\mathbf{x}}))(\mathbf{Az} + (\mathbf{b} - \bar{\mathbf{x}}))^T}\right] $$ Working out individual terms:   1. The first term: $\overline{(\mathbf{x} - \bar{\mathbf{x}})(\mathbf{x} - \bar{\mathbf{x}})^T} = \mathbf{P}_x$    2. The second term: (utilizing ''smart zero'' $-\mathbf{A}\bar{\mathbf{z}} + \mathbf{A}\bar{\mathbf{z}}$ and the fact that $\overline{(\mathbf{x} - \bar{\mathbf{x}})\bar{\mathbf{z}}\mathbf{A}^T} = 0, \ \overline{(\mathbf{x} - \bar{\mathbf{x}})(\mathbf{b} - \bar{\mathbf{x}})^T} = 0$)   $-2\overline{(\mathbf{x} - \mathbf{\bar{x}})(\mathbf{Az} + (\mathbf{b} - \bar{\mathbf{x}}))^T} = -2\overline{(\mathbf{x} - \mathbf{\bar{x}})(\mathbf{A}(\mathbf{z}-\bar{\mathbf{z}}) + \mathbf{A}\mathbf{\bar{z}} + (\mathbf{b} - \bar{\mathbf{x}}))^T} = -2\mathbf{AP}_{zx}$   3. The third term: (''smart zero'' $-\mathbf{A}\bar{\mathbf{z}} + \mathbf{A}\bar{\mathbf{z}}$ again! <-- this hit me when writing this post, *sigh*, fuuuu)   $\overline{(\mathbf{Az} + (\mathbf{b} - \bar{\mathbf{x}}))(\mathbf{Az} + (\mathbf{b} - \bar{\mathbf{x}}))^T} = \mathbf{A}\bar{\mathbf{z}}\bar{\mathbf{z}}^T\mathbf{A}^T + \mathbf{A}\mathbf{P}_z\mathbf{A}^T + 2\mathbf{A}\bar{\mathbf{z}}(\mathbf{b}-\bar{\mathbf{x}})^T + (\mathbf{b} - \bar{\mathbf{x}})(\mathbf{b} - \bar{\mathbf{x}})^T$ This way I am able to account for all the terms in the objective function $J$ and, of course, everything is solved. Although, I still have one question:   Authors in [1] have $\mathbf{A}(\mathbf{P}_z + \overline{\mathbf{zz}}^T)\mathbf{A}^T$ whereas what I obtained is $\mathbf{A}(\mathbf{P}_z + \mathbf{\bar{z}}\mathbf{\bar{z}}^T)\mathbf{A}^T$. But I think that $\overline{\mathbf{zz}}^T \neq \mathbf{\bar{z}}\mathbf{\bar{z}}^T$ and therefore, unless authors of [1] used some other devilish approach, there's a typo in the book. What do you think? * * * [1]: F.L. Lewis, L. Xie, D. Popa, Optimal and Robust Estimation: With an Introduction to Stochastic Control Theory, 2008


title: PAC learning theory and lower bound on the amount of input samples

text: I am trying to answer the following question: "How much (binary) data do I need for my learner to have seen every variable of the dataset at least once?" In my set-up I am feeding my algorithm binary vectors (i.e. with all elements equal to either 1 or 0), these vectors have a known 'density' (average amount of ones) which - for the purpose of answering this question - are uniformly constant (ok) or follow a long tailed distribution (better). I have tried looking at it from the perspective of combinatorics but this was harder than expected. I suppose this question must have been asked before, but I have not been able to find any references so far. In "A theory of the Learnable" by Valiant, I read that: > Let L(h,S) be the smallest integer such that in $L(h,S)$ independent > Bernoulli trials each with probability at least $h^{-1}$ of success, the > probability of having fewer than $S$ successes is less than $h^{-1}$. [...] > PROPOSITION: For all integers $S > 1$ and all real $h > 1$. $$L(h,S) \leq 2 > h (S + \ln (h))$$ This can be translated to an upper bound for my question given that each feature is assumed to be drawn from an independent Bernoulli trial, but not a lower bound. Does anyone know of other related work that could point me towards a lower bound?


title: Categorical independent variables and binary dependent variable? GLM?

text: I've done a literature review. The variable `treat` takes only TRUE and FALSE values. It says whether they applied a given treatment in a given article. I want to know if the probability of applying the treatment changes depending on the `journal` and on the `article.type`. By the way, I should say that all journals do not have the same frequency of article type. Here is a data to use as an example:               df = data.frame ( treat=c(T,F,T,T,F,T,T,T,F,T,T,T,T,F,T,T,F,T,F,T,T,T,F,F),     journal=rep(c('a','b','c'),8),     article.type=c('a','b','c','a','b','a','a','b','c','a','a','b','c','a','c',     'a','a','a','b','a','c','a','a','c')     )      **First question:** What statistics test should I use ? I thought of using GLM (generalized linear model) with a binomial error distribution. Does it seem good to you ? **Second question: Already Answered**               summary(glm(treat~journal*article.type, data=df, family='binomial'))          Coefficients: (1 not defined because of singularities)                              Estimate Std. Error z value Pr(>|z|)     (Intercept)             1.957e+01  4.390e+03   0.004    0.996     journalb               -1.997e+01  4.390e+03  -0.005    0.996     journalc                8.429e-09  8.781e+03   0.000    1.000     article.typeb          -3.913e+01  1.162e+04  -0.003    0.997     article.typec           4.134e-08  1.162e+04   0.000    1.000     journalb:article.typeb  3.884e+01  1.162e+04   0.003    0.997     journalc:article.typeb  3.913e+01  1.756e+04   0.002    0.998     journalb:article.typec         NA         NA      NA       NA     journalc:article.typec -1.916e+01  1.388e+04  -0.001    0.999      When running this, I get several p-values per variable. One p-value for "article.typeb" and one for "article.typec" for example. What does it mean ?


title: How to test for heteroskedasticity in Random Effects model?

text: I know how to test for heteroskedasticity using pooled OLS? But how can I do it after I have run a Random Effects model? If I find heteroskedasticity in pooled ols does it mean it is also present in the RE model?


title: Kendall's coefficient of concordance (W) for ratings with a lot of ties

text: I have read that Kendall W should be avoided when it comes to deal with non- rankings especially for rating scales which tend to have a lot of ties. Yet posts here seem to suggest it for ratings. As stated in this post I have a small study of 21 respondents, who rated some items from 0-5, with 0 being unimportant and 5 being very important and I'm looking for measures of agreement for specific respondents. I am not looking for Absolute agreement. Whilst ICC was suggested as a possible solution, there is an issue with the use of the F test in this case, given the small number of respondents. What are your views on Kendall W in this case?


title: Evaluating differences between repeated measurements of one group to another

text: I'm having difficulties understanding what type of statistical testing I need for my questions. My measurements include the following:               Organism 1     Location 1      5 repeated measurements      Organism 1     Location 2      5 repeated measurements      Organism 1     Location 3      5 repeated measurements        .       .       .       .     Organism 10     Location 1      5 repeated measurements      Organism 10     Location 2      5 repeated measurements      Organism 10     Location 3      5 repeated measurements      The questions I would to ask include:   1. What is the expected variation from a single location   2. What is the expected variation from a single organism   3. What is the expected variation between organisms    4. Is the variation within a single organism significantly different from the variation from a single location    5. Is the variation within a single organism significantly different from the variation between organisms  I have been reading up on $p$ and $t$ tests but I cannot get my head around whether these are the tests I'm required to use. I'm also unsure whether I need to compute the mean of means for the first two questions and simply the mean of all 150 measurements for the third question.


title: Reference paper and/or books about spatial data analysis, possibly bayesian

text: It is for my master thesis, I need to go more in depth. Any suggestion about recent works is really welcome.


title: Designing an experiment: Geometric or Bernoulli data

text: I have some process that succeeds or fails with probability $p$. I can do distinct simulations to estimate $p$:   1. Run $N$ simulations of a single process, record $N$ samples of a $\text{Bernoulli}(p)$, do GLM with binomial family to estimate $\log(p/(1-p))$.   2. Run $N$ simulations in which I continue running the process until a success, record $N$ samples of a $\text{Geometric}(p)$, do GLM with gamma family to estimate $p$. What is the theoretical difference? Which way is better? ~~I thought (1) is nicer but then I can't get $p$, only $\text{logit}(p)$, right?~~ (see answer below). UPDATE: I do regression because there is a simulation parameter `t` and I want to determine if `p` depends on `t` and how. Also, the above `N`s are not the same, because each simulation in the geometric design is roughly equivalent to `1/p` simulations in the Bernoulli design. But this is part of the question - is it better to make each trial on its own (Bernoulli) or to make sequences of trials? How does this design affect the power/significance/other statistical measures of both estimation of `p` and estimation of the regression parameters?


title: Confidence interval for median

text: I have to find a 95% C.I. on the median and other percentiles. I don't know how to approach this. I mainly use R as a programming tool.


title: Subset data by month in R

text: I am working with a time series of meteorological data and want to extract just the summer months. The data frame looks like this:               FECHA;H_SOLAR;DIR_M;DIR_S;VEL_M;VEL_S;VEL_X;U;V;TEMP_M;HR;BAT;PRECIP;RAD;UVA;UVB;FOG;GRID;     00/01/01;23:50:00;203.5;6.6;2.0;0.5;-99.9;-99.9;-99.9;6.0;-99.9;9.0;-99.9;-99.9;-99.9;-99.9;-99.9;-99.9     00/01/02;23:50:00;235.5;7.5;1.8;0.5;-99.9;-99.9;-99.9;6.1;-99.9;8.9;-99.9;-99.9;-99.9;-99.9;-99.9;-99.9     00/01/03;23:50:00;217.4;6.1;1.4;0.5;-99.9;-99.9;-99.9;7.0;-99.9;8.9;-99.9;-99.9;-99.9;-99.9;-99.9;-99.9     00/01/04;23:50:00;202.5;8.6;1.8;0.5;-99.9;-99.9;-99.9;6.4;-99.9;8.8;-99.9;-99.9;-99.9;-99.9;-99.9;-99.9     00/01/05;23:50:00;198.5;7.1;1.8;0.5;-99.9;-99.9;-99.9;5.4;-99.9;8.8;-99.9;-99.9;-99.9;-99.9;-99.9;-99.9      I have found some examples of time subsetting in R but only between an starting and end date. What I want is to extract all data from a month for all years to create a new data frame to work with. I can create a zoo time series from the data but how do I subset? zoo aggregate? Thanks in advance for your help.


title: Pre-test, post-test multiple comparisons

text: This feels like it should be a very straightforward stats analysis but I'm struggling to find a solution. I have a data set made up of pre- and post-test continuous data values taken from multiple groups i.e.: Group 1 (pre-test, post-test) Group 2 (pre-test, post-test) Group 3 (pre-test, post-test) etc.. I want to see in which groups there is a significant difference across the pre- and post-test measures. The most simple analysis would be a paired t-test, however I have many groups and am concerned that multiple paired t-tests will leave my results vulnerable to familywise errors. An alternative might be to use a repeated measure GLM. I ran such an analysis on SPSS and sure enough came up with a between(group)*within(test-effect) interaction. I'm now stuck as to how to find in which groups the within-subject comparison is significant. I can only think of splitting the groups up and running multiple comparisons in the estimated marginal means, with a Bonferroni correction. Only the Bonferroni correction does not work with just a single within-subject comparison. Are there any kinds of family-wise error corrections I can do for the follow- up within-subject comparisons when there is only one-level of comparison? Thanks!


title: Cross-validation with neural networks yielding worse results than a standard neural network

text: **Summary:** when using a 10-fold cross-validation procedure where each training set is used to generate N bootstrap samples for processing with NNs. How do I provide my NN with correct sequence and elements of my response variables series to match the length of my input time series so that I follow proper cross validation procedure using bootstrap and generate performance, similar (hopefully better), than using a non CV ANN. The purpose of using a cross-validation bootstrap method is to provide an estimate of bias and variance so I can create confidence intervals for my prediction.I am quite confident that my mistake is theory based. I am not sure if I am handling the response variable correctly when a) partitioning my response variable according to the partitioning of my input series using CV, and b) when generating my input and target series to be input into the NN. The problem continued.... I'm working on a project with time series data to fit a hydrological time series using Neural Networks. I am trying to fit my one explanatory variable (X) to a response variable (T), both series are of size 429x1. I have ran a simple NN using a script generated by the nftool and achieve quite nice results, with low MSE and R values above 0.9 for training, validation, and testing sets. I'm new to using ANNs so the only variation in model architecture I am performing is altering the lags in the input time series as well as the amount of hidden layers, for both I range between 1-9, therefore 81 different model architectures. I chose the best ANN architecture based off a combination of coefficient of determination, mse, mae, and pers index. I now wish to run the best ANN architecture using a 10-fold cross-validation (CV) procedure ensembled with a bootstrap approach. Now I am running into poor performance issues when running my bootstrapped cross-validation sets through the neural network, I am receiving extremely poor results in all training, validation, and testing sets. The worst performing sets are validation and testing, providing R values around -0.1, and the training sets providing R values around 0.2. It is obvious to me I have a problem as a cross-validation bootstrap ensemble NN should almost always perform as good as, or at least, the same as the NN with the same architecture. One common occurrence I have noticed amongst my output responses from the ensemble NN is that there is very small variance in the outputs compared to the inputs and the approximated output bias at each point in the time series hovers around the median value of my input time series. I have also confirmed this by taking a histogram of my CV training sets versus the outputs provided by the NN. I have ran through my code quite a bit to try and figure out what could be responsible for such poor performance and I am unsure if I am feeding my ensemble ANN the correct response variable for each training/test set case. My basic approach is to generate my 10 different input training/test sets using the cvpartition fcn. When using this function I create a simple for loop to place each of the 10 different train/test sets into a cell array. For every unique train/test set of my input series I use their indexes with respect to the original time series and create unique train/test sets of my response time series to mirror the input time series. Is this an incorrect action? Should I always use the same response series when testing each individual train/test set of the input series provided by cvpartition? I ask this because in my 10-fold CV bootstrap procedure I resample each of my unique training sets N times to create my bootstrap samples. I then find the set difference of elements between each unique test set and each of the bootstrap samples created from it, creating N sets of elements to be used for validation in the Neural Networks. So as you can imagine my validation sets with always vary in size, while my training sets will always be of consistent size because they are sampled at the same length as the training sets. I then pair each unique test along with the bootstrapped training sets, and the validation sets created by their set difference and use this as my input series to the NN. For my response variables I create target data for the NN by combining the unique train sets (again not sure if this is proper procedure), with the unique test set, and then indexing my response time series with the set difference between the bootstrap samples and unique train sets of the input series to determine the elements of the response time series to include in the validation set for the target data. I then input both these series into the NN so it can fit my data. I am quite confident that my mistake is theory based. I am not sure if I am handling the response variable correctly when a) partitioning my response variable according to the partitioning of my input series using CV, and b) when generating my input and target series to be input into the NN. Has anyone ran into a similar problem before, or can anyone be so kind as to point me in the correct direction to fixing these poor NN results? I greatly appreciate any help, and if code is needed do not hesitate to ask. * * * @Douglas Zare I do not alter the weights in any direct way in my code, I let the network randomly initialize and update the weights as it operates. The only architecture of the NNs I am altering is the amount of input and hidden nodes. I made a mistake above, I do not have upto 9 hidden layers, I only have one hidden layer, but with variable node sizes from 1-9. If I have not touched the weights in any form, could this still be the reason why my network is training poorly? I sincerely think I am making a procedure mistake when I am trying to create my target vector to feed the NN. My original time series is only 429 values, but with the 10-fold CV bootstrap technique, my input vector to the NN will always be larger due to including the variably sized validation set, as this is calculated as the set difference between a particular bootstrap resample of each 10-fold CV's training sets. So what I am doing to compensate for these additional data points, is indexing the elements in my input validation set, and then using that index in my target series to extract those elements and place in my validation set of my target values. I was always under the impression that you are never to alter your response variable, but I am unsure how I am suppose to create a valid target data set when the amount of data points in the data set exceeds the original time series.


title: Interaction terms and BUGS

text: It occurred to me that how one writes an interaction term in BUGS depends on the types of predictors that are "interacting". Suppose data that are long (rather than wide). Then, a model with an interaction between two factors, A and B, as in two-way ANOVA, may look something like               y[i] ~ dnorm(mu[i], tau)     mu[i] <- A[a_level[i]] + B[b_level[i]] + A.B[a_level[i],b_level[i]]     A.B[,] ~ dnorm(0, ab.tau)      However, a model with an interaction between a factor A and a continuous predictor C may look something like               y[i] ~ dnorm(mu[i], tau)     mu[i] <- A[a_level[i]] + betaC*C[i] + betaX[a_level[i]]*C[i]     betaX ~ dnorm(0, betaX.tau)      So, a few questions:   1. Are these interactions written correctly?   2. What other important interactions are there? Obviously, two continuous predictors is another case.   3. Are there differences for how one would write this for linear regression vs. logistic regression?   4. Are there differences in the priors that should be specified for different types of interactions? This may need to be a community wiki, not sure.


title: ISOMAP for classification

text: I implemented ISOMAP and I plan to use it as a feature extraction technique for a classification task. My problem is that although I can map the training data into a lower dimensional space, how can I use the same mapping for a separate test dataset?


title: How to test for multicollinearity among non-linearly related independent variables

text: I have a health outcome (measured as a rate of cases per 10,000 people in an administrative zone) that I'd like to associate with 15 independent variables (social, economic, and environmental measures of those same administrative zones) through some kind of model (I'm thinking a Poisson GLM or negative binomial if there's overdispersion). From scatterplot investigation I know there is multicollinearity among the variables that I need to investigate and deal with. I am not against removing variables, and within any problematic combinations can justify choosing one variable over the other(s) for reasons of interest or cost/ease of collection. In the past I used the correlation matrix to detect multicollinearity, but I've been reading around on this site and discovered VIF and the condition index/number, which seem to be generally accepted as better options. My question is how any of these measures work with non-linear correlations, which is what I have between nearly all of my variables (determined graphically). What are my options for evaluating multicollinearity besides Spearman correlations? Aside: I realize this is probably a separate question, but in case it's relevant here, I also have a lot of non-monotonic correlations among my variables that I don't know what to do with... Thanks for any help!


title: Two-dimensional distribution and related variance computation

text: I have a simple statistical question (forgive me if I use statistical terminology in a wrong way) Suppose I have a random vector with two components, $(x_1, x_2)$, where $x_1$ can take values from $1$ to $n$ and $x_2$ can take values from $1$ to $m$ and I have a two-dimensional distribution $P$ as a $n \times m$ matrix $(i, j)$-th element of $P$ is the probability to observe $(i, j)$ sequence. Now, I have observations $x^{(k)}_1, x^{(k)}_2$ and can estimate $\widehat{P}$ from the observations (as statistical frequencies)). Now an interesting part. I have a (fixed) vector $v$ and I need to estimate, how well **the matrix-by- vector product** of $\widehat{P}$ by $v$ is estimated: $$ E || (P - \widehat{P})v ||^2 = ?$$


title: How to compute the Standard Normal homogeneity test, Buishand range test, Von Neumann ratio test and the Pettitt test

text: I am not very new to R as I use it for my daily data analysis needs. I have been searching for an R package that can compute Standard Normal Homogeneity Test, Buishand range test, Von Neumann ratio test and the Pettitt test for daily precipitation series. The aim is to capture breaks in years where there is inhomogeneity. Is there a way these 4 methods can be done in R presumably within an R package? XLSTAT does it but it is expensive and time consuming.


title: Multiple covariates for each fixed effect

text: I'm analyzing data from a classical intervention design. Subjects were divided into groups, undertaking different interventions. Each subject was measured using the same tests before and after the intervention. As always, the very reason to have multiple groups is that we want to isolate the effect(s) of interest. For example we have a three-groups study where we want to estimate treatment and placebo effects:   1. Retest group: `measured gain = retest` (no intervention, just pre- and post-test).   2. Placebo group: `measured gain = retest + placebo`   3. Treatment group: `measured gain = retest + placebo + treatment` # Basic model This is easy to analyze in `R` using e.g. lme4::lmer with the mixed effects model `y ~ 1 + retest + placebo + actual + (1|id)` where the three fixed effects are binary (0/1) indicators about which effects this subject is experiencing and the random part captures subject-specific variance. Or as an equation: $y_i =\beta_{pretest} + \beta_{retest}x_{retest_i} + \beta_{placebo}x_{placebo_i} + \beta_{treat}x_{treat_i} + u_i + \epsilon_i$. ... with $i$ indicating a specific subject and $u_i$ being subject-specific offset (random effect). Note that for the pretest, all $x=0$ for the pre-test, which then is modelled as the intercept $y_i = \beta_{pretest} + \epsilon_i$. # Adding covariates Now, these effects might be moderated by covariates like age ($x_{age_i}$), gender ($x_{gender_i}$) etc. I guess it's possible that the effects are moderated in different ways by these covariates. I.e. the placebo effect might increase with age whereas the treatment effect decreases with age. Adding covariates would look like this equation: $y_i =\beta_{pretest} + (\beta_{retest} + \beta_{age_{retest}}x_{age_i} + \beta_{gender_{retest}}x_{gender_i})x_{retest_i} + ... + u_i + \epsilon_i$ So the $x_{effect}$s are on/off switches for everything related to an effect. (thoughts about alternative ways to do it are welcome). # Question: what lmer formula would do this? What would be the `lmer` equivalent? Would `y ~ (retest + placebo + treat)*(1 + age + gender) + (1|id)` do the trick? I fear that the effect-betas would be multiplied with the covariates rather than added to them using this R formula.


title: SVM hyperplane equation

text: I am trying to understand the hyperplane of the SVM algorithm. My problem is that different sources state different equations. Which equation is the right one and if both are right, why?   * Wikipedia: $w\cdot x-b = 0.$   * Intelligent Data Analysis and Probabilistic Inference notes: $ w\cdot x+b = 0.$ Furthermore, I thought the normal form of a plane looks like this: $$ w\cdot x + w\cdot a = 0.$$


title: Derivative of a linear model

text: **Question**   Is there such concept in econometrics/statistics as a derivative of parameter $\hat{b_{p}}$ in a linear model with respect to some observation $X_{ij}$?   By derivative I mean $\frac{\partial \hat{b_{p}}}{\partial X_{ij}}$ - how would change parameter $\hat{b_{p}}$ if we changed $X_{ij}$? **Motivation**   I was thinking about a situation when we have some uncertainty in data (e.g. results from a survey) and we have enough money to obtain precise results only in a single observation, which observation should we choose?   My intuition is saying that we should choose the observation that might change parameters the most, which is equivalent to the highest value of the derivative. If there are any other concepts feel free to write about them.


title: Which test to use for paired, nonparametric, categorical data?

text: We asked a group of subjects to tell us their preferences for a given procedure, of which there were 4 choices. We then provided them with educational material and asked them again of their preference. I would like to compare the two groups and see if there is a statistically significant difference between the two groups. As far as I know, I am making the assumption that the distribution of choices is non parametric since I assume that these subjects had a pence to begin with. Additionally I am polling the same group so the test I would need to use would have to be a paired test. Which test should I use?


title: how to handle occlusion for face recognition?

text: If there some occlusion on face (such as sunglasses, mask, scarf), the recognition rate will decrease steeply? How to handle this case? I have do a survey, the PCA reconstruction method have been mentioned by some literature. I don't know whether or not this method is the state-of-the-art method. Is there other method for address occlusion problem?


title: Does it make sense to partially scale the data matrix X in regression?

text: For some reason my supervisor wants me to centre only the independent variables that are used in interaction terms. I have never heard of such a practice. Does it make sense to partially centre data matrix $X$, that is, centering only some of the columns of data matrix $X$, not all of them? Thank you.


title: Predict binary outcome with R

text: I posted my original question here, but was told that I should have posted at Cross Validated, below is the link to the original question: http://stackoverflow.com/questions/20943584/predict-binary-outcome- with-r?noredirect=1#comment31453858_20943584 thanks!


title: Odd timing behavior with glmnet and n_lambda

text: I'm using the `glmnet` package in R for a relatively large dataset (~50k observations, ~1000 predictors). I'm observing weird timing behavior for different values of the `n_lambda` parameter. Using `cv.glmnet` (5-fold CV done in parallel), smaller values of `n_lambda` correspond to larger execution times, as measured by 'elapsed' in `system.time()`. This is counter-intuitive to me because it seems that fewer values tried for the shrinkage parameter should correspond to shorter execution times. Am I missing something here? Here's the function call:               system.time(cv.fit <- cv.glmnet(x=X, y=Y,                                      type.measure="class",nfolds=5,                                     parallel=TRUE,standardize=FALSE,                                     nlambda=nlambda,family="multinomial"))      Here's the output:               [1] "Fitting the model..."     [1] "Number of lambda values: 10"          user    system   elapsed          68544.863 2.764   31072.312           [1] "Fitting the model..."     [1] "Number of lambda values: 100"      user    system   elapsed       7061.27 4.360    24240.089           [1] "Fitting the model..."     [1] "Number of lambda values: 500"      user    system   elapsed       19164.12 26.773  8712.563


title: Correction for multiple comparisons

text: I have a 2x3x4 repeated measured anova. I have a significant 3-way interaction, and I want to make sure that I am using the correct post hoc comparisons and not violating any key statistics theory. I have run the statistics in SPSS and have adjusted for multiple comparisons using SIDAK but I want to make sure I understand how many adjustments are being made and if the p values are being corrected appropriately. My post hoc tests have analysed the following: AxB at each level of C AxC at each level of B BXC at each level of A I am just trying to determine what correction factor is appropriate. For example b x c at each level of A compares 2 means 12 times. Here am I not adjusting because it is 2 means, or I am dividing the alpha value by 12 as there are 12 comparisons. Similarly for a x c at each level of B I am comparing 3 means (as b has 3 levels) at each combination of a x c (8 combinations). Here would I be dividing the alpha value by 3 or by the combinations x the different means hence 24). Similarly, for axb at each level of c am comparing 4 means (5 comparisons total) at each combination of a x b (6 total).Here would I be dividing the alpha value by 5 or by the combinations x the different means hence 30. Alternatively would i be dividing the alpha value by the sum of all the above combinations? Thanks


title: Combining Exploratory Factor Analysis with Random Forest for Binary Logistic Regression Feature Selection

text: For those of you familiar with Exploratory Factor Analysis (EFA) and Random Forest (RF), I have recently had an idea of combining these two methods to reduce the number of potential predictor variables for use in a parsimonious binary logistic regression model. For the purposes of this post, assume large _n_ (200k or more) and 1000 potential predictor variables. To employ this idea, the first step would be to perform an EFA with all potential predictor variables using `proc varclus`. Additionally, using `randomForest` to rank all potential predictor variables by `IncNodePurity` (Gini Index). After these two methods are independently used, I propose retaining the variable with the largest `IncNodePurity` (from RF) within each factor (from EFA). Does anyone have any thoughts/concerns with this methodology (or lack thereof) for feature selection? I am aware that this "picking and choosing" of methods may be complete garbage, but I had this random thought and wanted to share. Thanks!


title: Hierarchical clustering with NCD as the distance metric

text: I am trying to cluster a bunch of executable files, and I want to use NCD ( normalized compression distance) as the distance metric. Is there any software package which lets me do that? Update: I'm using a python script which produces a distance matrix for example for 4 instances like:               [[1.0004239875837764, 1.0004587405627736, 1.0012468878082528, 1.0004537475969288],      [1.0001477300117407, 1.0003889794346572, 1.0007857384580077],      [1.0005753695194108, 1.0009174811255472],      [1.0008124881847325],      []]      Thanks


title: Performing T-Test with elements with multiple attributes

text: For an archaeology project, I have two assemblages of items, each measured in a variety of ways (mass, area, thickness, etc.). I want to test if the differences between the assemblages are significant. I am pretty sure I have to use the T-Test, but I am not sure how to apply it, taking into account all of the various attributes of each item. Do I perform the T-Test for each attribute?


title: Interpreting results from Sobol sensitivity analysis in R

text: I'm trying to use the sobol2007 model in the R sensitivity package. I'm doing runs on a model with 26 parameters, and using 2 sets of 500 monte-carlo samples to seed the analysis, and nboot=500. This results in 14k runs When I generate the sensitivity results, I'm expecting numbers in 0..1, but I get both * negative numbers * very high numbers (e.g. 10) When I plot the model, I can see that the confidence intervals are gigantic: ![Sobol sensitivity model output](http://i.stack.imgur.com/YfOd5.png) It feels like 14k runs should be reasonable to estimate this number of parameters, but maybe I'm being overly optimistic. Can anyone answer:   * am I totally unrealistic, and/or how many runs should I be using?   * can I read anything from these outputs, or are they essentially so noisy as to be worthless?   * how should I be selecting nboot - there's no guidance in the docs?


title: Confidence intervals for exponential smoothing

text: I'm using exponential smoothing (Brown's method) for forecasting. The forecast can be calculated for one or more steps (time intervals). Is there any way to calculate confidence intervals for such prognosis (ex-ante)?


title: How do I analyze the monthly variance of discrete series?

text: My company gathers all sorts of series for page views, unique visitors and many other user data (collected by our web servers). As I see it, these sets of series are discrete since they are all collected by computer systems and software and there is no 'intermediate' value between the collected numbers. Also because of that I don't have a way to measure its 'estimation error'. These series are compared monthly and my company moves forward according to the percentage increase/decrease in the values. But my managers never take in consideration its standard deviation nor its 'possible error'. The only statistical analysis they do is the percent variation between series. What would be a sounder, proper way to analyses the monthly variation between each series? I aim to derive the 'real' variation in the number of visitors, page-views, etc.


title: Elementary treatment of the Weibull distribution

text: I am looking for a reference to an elementary (or at least fairly simple) treatment of the Weibull-distribution. I have a bright high school student doing a project on wind mills and it turns out that wind speeds follow a Weibull distribution. Ideally I am looking for a derivation of the mean and standard deviations of a Windmill. A reference a proof of the fact, that the distribution in fact is a distribution would also be nice to see.


title: How to test the randomness of a sequence of real m-tuples?

text: 1) Let $\{a_i\}_{i = 1, ...,n}$ be a real sequence generated by a pseudorandom number generator. How can I check the randomness of such real numbers using, for instance, Runs test? 2) Assume three independent sequences of real numbers $\{a_i\}_{i = 1, ...,n}$, $\{b_i\}_{i = 1, ...,n}$ and $\{c_i\}_{i = 1, ...,n}$ are generated by the same PRNG using different seed keys. How can I infer the randomness of a sequence of _m_ -tuples, for instance 3-tuples $\{(a_i,b_i,c_i)\}_{i = 1, ...,n}$? 3) How can I translate the randomness of each coordinate to the randomness of 3-tuples?


title: AB testing with R

text: I'm new with AB testing and i want to do this with R (free and easy to learn software). Can anyone help with sample code in R for AB testing? Thanks Ariel


title: What is the UMVUE for $\sigma^2$ in $\mathcal N(0, \sigma^2 )$?

text: By using the exponential class factorization theorem, I came up with $Y = \sum (x_i)^2$ to be the complete and sufficient statistics for $\sigma^2$ . Using this sufficient statistic as a condition, would I need to get the MLE to derive the MVUE for $\sigma^2$ ? I don't know how to proceed to get the MVUE.


title: Wilcoxon-Mann-Whitney Output

text: I proposed a question early on infection rate data in 6 different hospitals for two years (N=12). After speaking with several people I decided that a Wilcoxon-Mann-Whitney (Mann-Whitney U, Wilcoxon rank-sum) test is the best bet. What I am concerned about is whether I can use rates and what exactly the output means (I've never had to run this test in SAS so I'm not even sure I'm putting it in correctly). Two questions arise: 1.) Will this test show that there is no difference/difference between the 2011 infection rate and the 2012 infection rate? 2.) I was told that the exact Wilcoxon test is the best with small sample sizes so in running this along with the test I get a One-sided p-value of .03, but a two-sided p-value of .06. Do I interpret this the same way as any other test? And is there a CI associated with this test or do you strictly report the p-value and leave it at that? MORE INFO: Rate = number of infections/ total patient days I feel this can also be run as a paired test which would lead me to the difference between the sign and signed rank tests that appear when you run Proc Univariate in SAS.


title: How to calculate expected win percentage where only 2 of 5 values are known?

text: The idea here is that you have a 5 vs. 5 game where each player is using a unique character (henceforth 'hero'), and thousands of matches of this game have been played. The goal of the analysis is to look for the best 2-hero combinations- those that perform better than their expected value. Some heroes work better when combined with other specific heroes, and some get in each others' way and perform worse. Without normalization, the best two- hero combination is unsurprisingly simply the #1 and #2 heroes in terms of solo win percentages. This doesn't yield any meaningful result, so we need to normalize/standardize the data. Here is the data set in question: Spreadsheet. It's 97 choose 2, but doubled due to duplicates. In short, my question's goal is to determine the correct formula to be used for **Column E**. Keep in mind that it's a 5 vs. 5 format, so the other 3 heroes that are teammates of the combination and the enemy 5 heroes are undetermined (but can probably be assumed at a value of 50% win). So for example, what is the expected normalized win percentage for a team which contains heroes which when used alone have 61% and 59% win percentages (along with any other three unknown heroes)? My initial thought was .59+(1-.59)*.61, but that gives an overly high expectation of 84% win, likely because it doesn't account for the other 3 team members bringing the expectation down. I'm unsure of how to proceed- I've tried a half dozen other methodologies but none come very near the actual figures. Thanks.


title: Assumptions of generalized linear models

text: On page 232 of "An R companion to applied regression" Fox and Weisberg note > Only the Gaussian family has constant variance, and in all other GLMs the > conditional variance of y at $\bf{x}$ depends on $\mu(x)$ Earlier, they note that the conditional variance of the Poisson is $\mu$ and that of the binomial is $\frac{\mu(1-\mu)}{N}$. For the Gaussian, this is a familiar and often checked assumption (homoscedasticity). Similarly, I often see the conditional variance of the Poisson discussed as an assumption of Poisson regression, together with remedies for cases when it is violated (e.g. negative binomial, zero inflated, etc). Yet I never see the conditional variance for the binomial discussed as an assumption in logistic regression. A little Googling did not find any mention of it. What am I missing here? EDIT subsequent to @whuber 's comment: As suggested I am looking through Hosmer & Lemeshow. It is interesting and I think it shows why I (and perhaps others) are confused. For example, the word "assumption" is not in the index to the book. In addition, we have this (p. 175) > In logistic regression we have to rely primarily on visual assessment, as > the distribution of the diagnostics under the hypothesis that the model fits > is known only in certain limited settings They show quite a few plots, but concentrate on scatterplots of various residuals vs the estimated probability. These plots (even for a good model, do not have the "blobby" pattern characteristic of similar plots in OLS regression, and so are harder to judge. Further, they show nothing akin to quantile plots. In R, plot.lm offers a nice default set of plots to assess models; I do not know of an equivalent for logistic regression, although it may be in some package. This may be because different plots would be needed for each type of model. SAS does offer some plots in PROC LOGISTIC. This certainly seems to be an area of potential confusion!


title: Omitted dummy variable in panel gls?

text: In Stata, I have panel data for the years 2008-2012, across 7 firms. I assign a value corresponding to `id` for all firms. (`id` 1,3 are specialized firms, which is of interest) I have two dummy variables. One is an indicator for the 2008 global financial crisis and which is 1 for year 2008. The other is to indicate specialized firms and I give 1 only to `id` 1,3 and zero to others. But when I use `xtgls` (panel gls) the omitted dummy variable problem occurs. I really need those two dummies, so how can I solve this? Total number of observations is 8000, each firm has around 1100 observations in each year.


title: Interpreting percentage units regressions

text: I am using a panel of 2249 schools with data from 2002-2008. Some of the schools are single sex schools whilst others are mixed sex. Some background on my regression: Consider the determinants of school performance (as measured by the **percentage** of pupils in the school who get 5 A*-C grades at GCSE: **sch5ac** ). Start by using the following explanatory variables: - **stdschk2aps** : a measure of the academic performance of pupils when they enter the school. The variable has been standardised to have mean zero and standard deviation one.  - **schfsm** : the **percentage** of pupils in the school who are eligible for free school meals; this can be thought of as a proxy for the proportion of pupils from low income families.  - **schsen** : the **percentage** of pupils in the school classified as having special educational needs (i.e. requiring additional support). Schools receive extra funding to help support these pupils.  I am asked to estimate the effect of each of these variables on the dependent variable (sch5ac) for the single sex schools, using the within estimator. sch5ac=beta_0+beta_1*schsen+beta_2*schfsm+beta_3*stdschk2aps I get the following coefficients: beta_1= 0.215 beta_2= -0.357 beta_3= 6.935 Both schsen and schfsm are measured as a percentage, as is the dependent variable sch5ac. So would a 1% increase in schsen lead to a 0.215% increase in sch5ac? Since stdschk2aps is not measured as a %, but the dependent variable sch5ac is, would a unit increase in stdschk2aps increase sch5ac by 6.935 **percentage points?** Any help would be great, Thanks Kai


title: IT-AIC approach for mixed effects model

text: I was reading around the Information Theoretic-AIC approach of model selection where AIC values are used to select the candidate set of models. I am quite clear on this. My question is this: for mixed effects model, how to calculate the AIC values. For e.g. I have several random intercept and fixed slope models like this:               mdl1<-lme(y ~ x1,random=(~1|state),data=data,method="ML")      mdl2<-lme(y ~ x1 + x2,random=(~1|state),data=data,method="ML")      mdl3<-lme(y ~ x2,random=(~1|state),data=data,method="ML")      mdl4<-lme(y ~ x3,random=(~1|state),data=data,method="ML")      mdl5<-lme(y ~ x3 + x1,random=(~1|state),data=data,method="ML")      .     .     .     .     mdlx<-lme(y ~ xi,random=(~1|state),data=data,method="ML")       I can compare all these models to generate AIC value for each model and then select the best models accordingly. However I came across this paper which says that this AIC value in mixed effects model is not correct. It talks about a conditional AIC value for mixed effects model http://biomet.oxfordjournals.org/content/92/2/351.abstract I am not from a strong statistics background and this paper has too much mathematics and equations. I was hoping if someone could explain me how to best understand this and proably how to calculate conditional AIC values in R. Thanks a lot


title: Regression with dependent observations

text: Suppose I have the following dataframe:               id    time    sale_yn number  journey_nr  "nVariables with observation specific variables"     abc   10      0      1        1     abc   11      0      2        1     abc   12      1      3        1     abc   13      0      1        2     def   10      0      1        3     def   15      1      2        3     def   16      0      1        4     def   17      0      2        4     def   18      1      3        4      Every row is an observation from a particular person (distinguished by ID). Every person visits a shop a number of times (see number). After some while it is possible that a customer makes a purchase (in that case sale_yn == 1). The regression model needs to take into account the different journeys a customer has made in order to predict whether a sale happens or not. With logistic regression, each observation should be independent of eachother, but the visits that occurred earlier, also contribute to the purchase. What form of regression can take this into account?


title: Calculate combined odds ratio between two factors

text: I am using synergy factors to test synergy between two variables. The equation is simple enough, where $OR_{12}$ is the combined odds ratio for variables 1 and 2: $ SF = \dfrac{OR_{12}}{(OR_1 \times OR_2)}$ How do I calculate the combined odds ratio for variables 1 and 2? I have all of the data and will calculate odds ratios by $\exp(\beta_j)$. As a side note, the odds ratios will be adjusted for other variables so I can't just use the $4 \times 2$ table as described the linked paper. I'm guessing this is fairly simple but have not found a clear solution. Several questions on the forum address meta-analysis, but I do not want to make any mistakes or poor assumptions. I appreciate your help! Partial LR output:                              Estimate Std. Error z value Pr(>|z|)     var1          -14.09922  769.21550  -0.018  0.98538     var2           -3.38497    2.19931  -1.539  0.12378     var1:var2      -1.61902    0.57816  -2.800  0.00511


title: Probability of collision (two bivariate normal distributions)

text: I am trying to solve this problem on and off for the past couple of months but to no success. This was supposed to be a very small part of my PhD thesis in navigation but I guess I underestimated the problem. It sounded trivial at the beginning, but now I am not so sure. Lets say we have two ships, each with its own nominal position in 2D coordinates (mean). Due to errors in positioning systems we can only be certain that the ships are within 1 mile of the mean with 95% probability (normal distribution). Given these 2 positions and this probability distribution, what is the probability that the ships are within 5 miles from each other? Also, same question if the ship's probable position is an ellipse, not a circle. I asked some people and they told me that there are no analytic solutions. If that is really the case, please explain how to solve it numerically. As you can already tell, I come from engineering background, therefore my math is more than a bit rusty. I apologize in advance if the question is too vague or too trivial for this forum. I will be more than happy to explain in more detail if needed. I found this, but it is only for univariate case, and besides I don't know how to implement it in my case where I need to find the probability that the distance between two ships is less than 5 miles. I imagine this problem as a plane with two hills that intersect and the solution is the volume under the circle with diameter of 5 miles that is located somewhere between the two peaks of hills (means). Am I on the right track? Thanks


title: Time series behavior based on previous data

text: Apologies in advance for incorrect use of terminology - I do not have much of a background in statistics. Suppose I have the following behavior for daily stock returns: When $r_{i+1} > r_i$, then $r_{i+2} < r_{i+1}$ with probability $>50$%. Conversely, when $r_{i+1} < r_i$, then $r_{i+2} > r_{i+1}$ with probability $>50$%. Describe what is happening/define this phenomenon, and explain how one can profit from knowing this information. My knowledge/intuition in dealing with time-series is limited, so any help would be greatly appreciated!


title: What are the (best) methods for multiple comparisons correction with bootstrap for multiple glm models?

text: See the related, but old question: Correcting p values for multiple tests where tests are correlated (genetics). Multiple comparison methods based on bootstrap have the advantage of taking account of dependence structure of pvalues. Regression models are just a little more difficult to bootstrap, though. There have been several methods proposed. I've come across at least three ways of adjusting p-values in such problems. As far as I see it, there is no single "best" solution. The question is: **what methods are available and what are the advantages/disadvantages of each?** I propose the following definitions: $\vec{\theta_0}$ is vector of the hypothesized values of the parameters, which is assumed under the complete null hypothesis. $\hat{\vec{\theta}}$ is an estimator of the parameters, calculated on the original sample. $\hat{\vec{\theta}^*}$ is an estimator of the parameters, calculated on one of the bootstrap samples. $T$, $T^*$ are pivot statistics calculated on original data and bootstrap sample. All pivot statistic are assumed to be Wald, i.e. $T = \frac{\theta}{\operatorname{SE}\theta}$. $m$ is number of terms of interest in all regressions.


title: Total probability problem

text: I currently have log files all off different sizes, and I am trying to calculate the total probability of an event that occurs being an error. For example: > Log file #1: Has a total of 1000 entries and 20 of them entries are errors > so the conditional probability for that file would be 20/1000 = 0.02 > > Log File #2: Has a total of 2500 entries and 35 of them entries are errors > so the conditional probability for that file would be 35/2500 = 0.014 > > Log File #3: Has a total of 5000 entries and 60 of them entries are errors > so the conditional probability for that file would be 60/5000 = 0.012 I was wondering how I would get the total probability of an event that occurs being an error for all the log files combined? Any help would be much appreciated, and if you could use the example I gave to explain to me it would be great. The events in the log file are completely different, the log files are for completely different systems. Would it be `(amount of errors/total number of events) = 115/8500 = 0.0135`?


title: Which ANOVA is appropriate?

text: Bacterial grow rates were measured in 4 flasks, 3 had been inoculated with different antibiotics, 1 flask was kept as a control. Measurements were taken every 20 minutes for 3 hours. The experiment was repeated 6 times. I understand two way ANOVA would be inappropriate because it treats each time point as independent, but I'm not sure which other ANOVA would be suitable


title: How exactly does Wilks' Lambda distribution generalize the Hotelling distribution?

text: According to Wikipedia, Wilks' Lambda distribution generalizes Hotelling's distribution. I am having some problems seeing how this works. I can see how Hotelling's distribution generalizes Student's t-distribution (a RV distributed as Hotelling's law with $p=1$ is just the square of a RV distribution as Student's t), but cannot see how to get to Wilks' Lambda. Is there some setting of the parameters $p,m,n$ such that a RV distributed as Wilks' lambda is some transform of a Hotelling RV?


title: Independence assumptions in the classical regression model and higher moments

text: I refer here to a classical linear regression whose true representation is given by the equation: $y_i=x_i'\beta+u_i$, where as usual $x_i$ is a $K\times 1$ vector of independent explanatory variables, $\beta$ is a $K\times 1$ vector of parameters and $u_i$ the error term of the ith observation and by construction of the true representation $u_i$ is considered to be independent of $x_i$. My first question on that is basic: If we assume independence between $x_i$ and $u_i$ does this imply $E[u_i^d|x_i]=E[u_i^d]$ for $d=1,2,..$ - so to speak, can we generalize the Law of iterated expectations consequently to all moments? Secondly, what is the basis of many textbooks to claim that $x_i$ and $u_i$ are independent in the classical regression model? By independence, I don't mean uncorrelatedness, rather $f(u,x)=f(u)f(x)$. I am aware that the two key standard assumptions on which this small sample properties model is based on is strict exogenity $E[u_i|X]=0$, and normally distributed error terms $u_i|x_i \thicksim N(0, \sigma^2)$. Where does the claim of the independence between $x_i$ and $u_i$ originates from in textbooks? Is it implied by these two statements jointly together or is this a separate assumption? If the latter, why is it not mentioned exclusively in the set of assumptions for classical linear regression models which only assumes strict exogenity and normally distributed errors?


title: How to implement multiple GP submodels in PYMC

text: I'm hoping someone can give me some guidance on implementing Gaussian processes (GP) with PYMC. In particular, I'm not sure how to use multiple GP submodels properly within a single pymc model. More specifically, what I'd like to implement is a classifier for which the input is a time-varying signal (an audio snippet, in this case). For concreteness, the goal is to determine from the input signal whether or not a whale is present (see the Kaggle whale detection competition). So, there are two categories, 'whale' or 'no-whale', which I encode as 1, 0, respectively, and model with a Bernoulli random variable. On 'whale' trials, there is a whale call embedded somewhere in the audio signal. On any given trial we don't know when the call starts or ends, but we assume that the call is fairly stereotyped, and therefore varies little from trial to trial (but on each 'whale' trial we observe a noisy version of it). I'm trying to model the background (no-whale) signal using GP submodels which share a common prior mean function and covariance function, and the 'whale' signal using a single GP submodel with a different covariance function (and possibly different prior mean, but in what I've coded here all submodels are zero mean). On 'whale' trials, the call is embedded into the background signal, starting at 't_start' and lasting for 'call_dur' time units. 't_start' and 'call_dur' are both random variables. My question is whether or not I really need separate submodels for each no- whale signal or whether I can use a single no-whale submodel and tell it to draw a new realization for each trial. (see comments in code, below). The code below runs very very slowly -- it seems to take a long time to instantiate each GP submodel -- too slowly to be useful. Is there any way around this? Thanks for any suggestions! ---- code follows ----               from __future__ import division     import numpy as np     import pymc as pm     import pymc.gp as pmgp          # for concreteness, here are the values     # passed in through the params object     #input_T = 4000     #params.t_start_mean = 0.25 * input_T     #params.t_start_var = 0.1 * input_T     #params.call_dur_mean = 0.4 * input_T     #params.call_dur_var = 0.05 * input_T     #params.n_fmesh = 25     #whale_prior = 0.5          def make_hbayes_model(params, whale_prior, sigs, obs_wtf=None):              # The data in 'sigs' are a set of 1-dimensional signals (audio snippets) passed in as a         # numpy array of size sig_N x sig_T, i.e., sig_N trials, each of duration sig_T         # Trials can be one of two categories -- 'whale'(1), or 'no-whale'(0).         # On 'whale' trials, there is a whale call         # embedded somewhere in the audio snippet, on 'no-whale' trials there is not.         # When training the model, the true categories are passed in in 'obs_wtf',         # i.e., whale true/false, 0 for no-whale, 1 for whale.         sig_N, sig_T = np.shape(sigs)              fmesh_nowhale = np.linspace(0, sig_T, params.n_fmesh) # fmesh for the no-whale GP         # A guess at what fraction of the total snippet the embedded whale call lasts         call_dur_frac = params.call_dur_mean / sig_T          # fmesh for 'whale' GP         fmesh_whale = np.linspace(0, params.call_dur_mean, int(params.n_fmesh*call_dur_frac))              if obs_wtf is None: # when predicting             obs_wtf_flag = False         else: # when training the model             obs_wtf_flag = True              # The random variable for generating the category (whale or no-whale)         # 'whale_prior' is the known value (from the training data) of Pr(whale)         # 'wtf' is, of course, observed when training, but is what we want to predict when, uh,         # predicting         wtf = pm.Bernoulli('wtf', p=whale_prior, value=obs_wtf, observed=obs_wtf_flag)         # we need to know how many whales so we know into how many sound signals we need         # to embed a whale call         N_whales = np.sum(wtf.value)              # On 'whale' trials, we know there is a whale call embedded, but we do not know exactly where         # it starts or exactly how long it lasts. The full snippet is only 2 secs long, and the call is around         # 1 sec long, so we know it should start ('t_start') toward the beginning of the snippet. Note that         # 't_start' and 'call_dur' represent indices, not actual time (they are rounded below).         t_start = pm.TruncatedNormal('t_start', mu=params.t_start_mean, tau=1.0/params.t_start_var,                                      a=1,b=sig_T-2,                             value=np.zeros(N_whales)+params.t_start_mean)         call_dur = pm.TruncatedNormal('call_dur', mu=params.call_dur_mean, tau=1.0/params.call_dur_var,                                       a=1,b=sig_T-t_start,                              value=np.zeros(N_whales)+params.call_dur_mean)                                  # Just hard-code this for now, later make it learnable         sig_obs_var = 1              # Prior parameters of C. Straight out of one of the pymc examples         nu_0 = pm.Uniform('nu_0', 1., 3, value=1.5)         phi_0 = pm.Lognormal('phi_0', mu=.4, tau=1, value=1)         theta_0 = pm.Lognormal('theta_0', mu=.5, tau=1, value=1)                      nu_1 = pm.Uniform('nu_1', 1., 3, value=1.5)         phi_1 = pm.Lognormal('phi_1', mu=.4, tau=1, value=1)         theta_1 = pm.Lognormal('theta_1', mu=.5, tau=1, value=1)              # The covariance function for the no-whale cases         @pm.deterministic         def C_no(eval_fun = pmgp.matern.euclidean, diff_degree=nu_0, amp=phi_0, scale=theta_0):             return pmgp.NearlyFullRankCovariance(eval_fun, diff_degree=diff_degree, amp=amp, scale=scale)             # Mean function -- prior mean for no-whale cases         # mean zero is probably a good-enough guess for now         @pm.deterministic         def M_no(eval_fun = lambda x: 0 * x):             return pmgp.Mean(eval_fun)              # The GP sub-model(s) for the no-whale cases         # I'm not really sure how to implement this.         # I think all no-whale cases should share a common prior mean and covariance function,         # but each individual instance (trial) should have a different underlying signal         # (a draw from the prior mean), around which noise is added to get the observed data.         # Do I need a separate sub-model (as I've coded here) for each no-whale trial to get this         # behavior? I'd rather just tell a single sub-model to draw a new sample from the prior mean         # for each new trial, but I'm not sure how.              # use the container trick to hold all no-whale the sub-models         sm_no = np.empty(sig_N, dtype=object)          for si in range(sig_N):             sm_no[si] = pmgp.GPSubmodel('sm_no_%i' % si, M_no, C_no, fmesh_nowhale)                   # Cov function and prior mean for 'whale' cases         @pm.deterministic         def C_whale(eval_fun = pmgp.matern.euclidean, diff_degree=nu_1, amp=phi_1, scale=theta_1):             return pmgp.NearlyFullRankCovariance(eval_fun, diff_degree=diff_degree, amp=amp, scale=scale)                  @pm.deterministic         def M_whale(eval_fun = lambda x: 0 * x):             return pmgp.Mean(eval_fun)              # GP submodel for 'whale' cases. Here we'll pretend that it is sufficient to have a single signal         # for all whale calls around which we add noise to get the observed data. In other words,         # a single draw from the prior mean is sufficient for all 'whale' trials, in contrast to         # the no-whale trials above. Remember, though, that on 'whale' trials the call is embedded         # in a longer signal, so it is surrounded by a 'no-whale' signal. This is          # instantiated in get_sig_means(), below.         sm_whale = pmgp.GPSubmodel('sm_whale', M_whale, C_whale, fmesh_whale)                  # Given the underlying signals from the 'whale' and 'no-whale' GPs, embed         # a whale call on each 'whale' trial         @pm.deterministic()         def get_sig_means(wtf=wtf, sm_no=sm_no, sm_whale=sm_whale,t_start=t_start,call_dur=call_dur):                     sig_means = np.zeros_like(sigs)             w_ctr = 0             for si in range(sig_N): # loop over trials                 #print si                 if wtf[si] == 1: # this is a 'whale' trial                     this_start = int( t_start[w_ctr] )                      this_dur = int( call_dur[w_ctr] )                      chunk0_t = range(this_start) # indices for the pre-call chunk of signal                     chunk1_t = range(this_start,this_dur) # indices for the call chunk                     chunk2_t = range(this_start+this_dur,sig_T) # indices for post-call chunk                          # build the signal, i.e., embed a call in the background ('no-whale') signal                     # ?? It seems that calling, e.g., sm_whale.f(x) returns                      # values deterministically. That is what we want, given the code here,                     # but it seems inefficient to have all these separate no-whale submodels.                     # It would be better, it seems to me, if we could have a single no-whale                      # submodel from which we draw a new realization for each trial. How?                     sig_means[si,chunk0_t] = sm_no[si].f(chunk0_t)                     sig_means[si,chunk1_t] = sm_whale.f(chunk1_t)                     sig_means[si,chunk2_t] = sm_no[si].f(chunk2_t)                     w_ctr += 1                 else: # this is a no-whale trial                     # Nothing to embed, just use the whole no-whale signal                     sig_means[si,:] = sm_no[si].f(range(sig_T))             return sig_means              # To get the observed signals we just add iid Gaussian noise to each datum, given the underlying         # signal from the GP. Note, we always observe 'sigs'. When training we also observe 'wtf', so         # we know whether or not there is a whale. When predicting we want to predict         # whether or not there is a whale, given the 'sigs' and the model.         signals = pm.Normal('signals',mu=get_sig_means,tau=1.0/sig_obs_var, value=sigs, observed=True)              return locals()


title: About finding testing Errors

text: I have 63 rows and 17 columns in the cocomo81 dataset (see the information here). The first 16 columns are the inputs to the network and the 17th column is the estimated outputs. I took 2/3rd of the data for training and validation (i.e 42 rows), 1/3rd for testing (i.e from 43rd to 63rd rows). I sent 1st row for training with random weights and then 2nd row with updated weights (got when 1st row sent) and so on up to 41 rows (1-1000 iterations) and averaged the error got in each row and took as training Error-1. To get error i used (estimated output-actual output). After this I sent 42nd row for validation and got error in the same way as before and updated weight so got validation Error-1(no iterations). Next I sent rows 43 to 63 for training with the updated weights got in validation and got the error in the same way as before (no iterations here). I used back propagation for updating weights and sigmoid activation function. What I did is the right thing or not? My sir said that the validation errors and the testing errors should be less compared to training errors but I'm getting large values. So please tell me, whats wrong?


title: What is a good and efficient algorithm for a content based recommender?

text: I want to build a content based recommender in a restricted environment regarding cpu power and memory (to be specific: a mobile device, but it is not acceptable to build the recommender on a remote server). Which algorithm or heuristic is feasible in such a situation ?


title: How to interpret parameters in GLM with family=Gamma

text: I have a question regarding parameter interpretation for a GLM with a gamma distributed dependent variable. This is what R returns for my GLM with a log- link:               Call:     glm(formula = income ~ height + age + educat + married + sex +      language + highschool, family = Gamma(link = log), data = fakesoep)          Deviance Residuals:           Min        1Q    Median        3Q       Max         -1.47399  -0.31490  -0.05961   0.18374   1.94176            Coefficients:               Estimate Std. Error t value Pr(>|t|)         (Intercept)  6.2202325  0.2182771  28.497  < 2e-16 ***     height       0.0082530  0.0011930   6.918 5.58e-12 ***     age          0.0001786  0.0009345   0.191    0.848         educat       0.0119425  0.0009816  12.166  < 2e-16 ***     married     -0.0178813  0.0173453  -1.031    0.303         sex         -0.3179608  0.0216168 -14.709  < 2e-16 ***     language     0.0050755  0.0279452   0.182    0.856         highschool   0.3466434  0.0167621  20.680  < 2e-16 ***     ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1          (Dispersion parameter for Gamma family taken to be 0.1747557)          Null deviance: 757.46  on 2999  degrees of freedom     Residual deviance: 502.50  on 2992  degrees of freedom     AIC: 49184      How do I interpret the parameters? If I calculate `exp(coef())` of my model, I get ~ 500 for the intercept. Now I believe that doesn't mean the expected income if all other variables are held constant, does it? Since the average or `mean(age)` lies at ~ 2000. I furthermore have no clue how to interpret the direction and value of the covariates' coefficients.


title: Feature selection for classifier

text: I'm using a supervised machine learning algorithm on some big data. There is much more features than observations. To reduce the number of features, I would like to do some feature selection. However, there is one thing I'm confused about. May I do the feature selection _first_ , using data from _all_ observations, and then do the training and testing of the classifier?


title: How do I fit a set of data to Burr Distribution in R?

text: For example I have a data set as: (2042.044,2218.153,3670.893,5149.684,5533.429,7111.183,11041.569,15459.771,783.477,1701.520,40810.770,67905.857) How do I fit the above data in Burr distribution to compute its parameters in R? fitdist does not provide Burr distribution. Can I explicitly define my probability distribution function for the required computation? Please Help as I am a newbie to R.


title: Goodness of fit in a GLM with scaled deviance

text: On this page, I am interested in the section “goodness of fit”, which is near the bottom of the page and contains the table of deviance functions. The author states that the scaled deviance, i.e. $D^*= \frac{D(y, \mu)}{\phi} $ has a limiting $ \chi^2_{n -p} $ distribution, where $ n $ is the number of observation and $ p $ is the number of predicted parameters. He then goes on to say that in the case $ \phi $ is unknown, it can be predicted as $ \hat{\phi} = \frac{D}{n - p} $. If this is the case, wouldn't the scaled deviance be equal to $ \frac{D(y, \mu)}{\frac{D}{n - p}} = n - p $? I think I am misunderstanding the discrepancy between $ D(y, \mu) $ and $ D $ (without arguments). A similar thing occurs with the scaled Pearson's chi-squared statistic. Could someone elaborate on how to calculate scaled deviance in the case that $ \phi $ is unknown and how to proceed with the g.o.f. test?


title: How to allocate sample size to gain most power on a two stage analysis?

text: In stage 1 we test many SNPs, in stage 2 we try to confirm these SNPs. Usually we start with less sample, but then it's hard to detect a signal. For example, if we can afford 3000 samples, how to allocate these 3000 people in two stages?


title: Testing whether two groups formed from a continuous variable differ on a binary outcome

text: 18 people, who have been divided into 2 groups based on the semiquantitative expression of a protein (low vs high). Group 1 (low protein levels) has 6 people, Group 2 (high protein levels) has 12 people. In Group 1, 5 of the 6 people developed a specific disease, whereas only 2 of 12 in Group 2 developed this disease. I want to determine if there is a statistically significant correlation between low protein levels and the development of this specific disease. If someone could point me towards the tests I should use, it would be great.


title: Non-integer dependent variable in negative binomial models

text: I have non-nested count data that I've interpolated from one area to another based on the proportion of the area that lays in each. This is ZIP codes to counties, so most nest cleanly, with a few around the edges needing to be distributed in this way. I'm aware that it would be better to propagate this uncertainty through the model (for instance Bannerjee Carlin and Gelfand, Chapter 6, Section 2), but that's way, way out of scope for this project for now, and I believe the error will be quite small. This apportionment in proportion to area leaves me with a reasonable distribution of "counts" :                        0          1      (1,5]     (5,10]   (10,250]            0.30       0.35       0.22       0.06       0.06       The NBD models I've run on this have successfully converged, and report parameter estimates in the same basic ballpark as a Poisson regression on the same data. It does return a warning message about non-integer values, though. Is there any bias that results from running a negative binomial regression on count data that is non-integer? Specifically, I'm doing this in R, using `MASS::glm.nb`.


title: Joint PDF change of variables

text: I now understand how to conduct a change of variables for a marginal PDF. Now, given two functions that define parameter's spatially: $C_A(x)$ and $C_B(x)$, is it possible to construct the Joint PDF, $f(C_A,C_B)$? Note: if C_A and C_B are independent, I understand that $f(C_A,C_B)=f(C_A)f(C_B)$.


title: Interpretation lin-log regression where the covariate is log(x1 + 1) transformed

text: I have a lin-log regression model like $$Y = b_0 + b_1 \log(x_1 + 1) + e.$$ The distribution of $x_1$ is very skewed, thus I use the natural logarithm to get a more Gaussian like distribution. Because 3 out of 100 values have zero as entry I add a constant c, in my case plus 1, to avoid -Inf. The resulting estimation of $b_1$ is about -0.14. Without the constant the interpretation is clear: a 1% change in $x$ results in a $0.01\cdot b_1$ change in $y$. I struggle with the constant. How can I account for it in my interpretation? If I change the value of c I get, of course, other estimates. I have chosen + 1 because this results in positive log values (the values of $x_1$ are originally positive too). Or should I add a small value just to the three 0s? Many thanks in advance and, please, a non mathematical answer ;-) Marco


title: Problem using auto.arima() in R

text: I am using `auto.arima()` for forecasting. When I am using any in built data such as "AirPassengers" it is capturing seasonality. But, If I am entering data in any other format (in vector form or from an excel sheet) it is not detecting seasonality. Is there any specific format in which it detects seasonality or I am doing some thing wrong? Does data have to be entered in a specific format?


title: Can the OLS residual variance suggest a polynomial relationship?

text: I am trying to figure out whether from the following graph of the OLS residuals that the linear relationship does not hold, and that probably a cubic relationship would do better? Since both in the beginning as in the middle the variance is larger, and two bends in a regression could lower the variance for these areas. Or is that something that is not determinable from this graph? ![enter image description here](http://i.stack.imgur.com/PLSos.jpg) [edit after comments] The mean structure does suggest a quadratic term, and possibly a cubic term, however this last one seems small. ![enter image description here](http://i.stack.imgur.com/Ndwch.jpg)


title: Comparison of two regression models, relative importance of individual symptoms

text: We have a metric dependent variable (impairment, scale 0-40), and a larger number of symptoms of a mental disorder, each on a scale of 0 (nonexistant) - 3 (very severe). We try to find out whether symptoms are associated with different levels of impairment. We test this in two steps. First, we use a linear regression, predicting impairment by individual symptoms, in which symptoms explain about 45% of the variance of impairment. Second, we use the LMG metric of the R package RELAIMPO (Grömping, 2006; "relative importance") to allocate unique R^2 shares to each symptom. The result is that individual symptoms vary drastically in their explained variance of impairment. The estimate range between 1% and 20% explained. This fits nicely with **our hypothesis** : > Some symptoms are more impairing than others, because they are different > things (and not simply interchangeable indicators of the same latent). > Impairment depends on the _nature_ of symptoms (thus, sum-scores of symptoms > might obfuscate information). However, a reviewer put forward an **alternative hypothesis** : > Impairment depends on the _severity_ of the symptoms, irrespective of the > _nature_ of symptoms. This hypothesis would state that all symptoms that > explain a lot of variance are symptoms that have a high mean, and that > symptoms that explain moderate amounts of variance are average severe > symptoms, and that symptoms that explain only little variance are only > slightly debilitating. I don't quite know how to test this idea, and would very much appreciate input. ( _EDIT_ : A colleague of mine suggested to test a null hypothesis model that uses the score sum only (i.e. forcing all betas to be identical) against the model with individual betas (13 degree of freedom test), but frankly, I don't understand what that would do and how it would help).


title: Panel regression modelling in econometrics

text: I am trying to run a regression to model sovereign default risk of EU countries.   My data includes credit ratings (numerical value assigned 0-20), CDS spreads, Bond Yield Spreads, Bond yields, Inflation and GDP per capita. However, I am a bit confused in terms of setting up the regression. I have to do a panel regression. But I was also thinking of doing individual time series for each country as well. Does this make sense? Or should I just do the panel alone?   Also I know I have to take logs for GDP per capita, and possibly inflation. Is this correct?   Would I have to do that for any other variables? When do you take logs? what is the principle behind this? My lecturer did not explain this to us properly. Also when do we do lags or differentiating?   What is the significance and purpose of doing this and when are we required to do this?   Do I need to do any of these for my variables? I ask this because I keep getting heteroscedasticity when running the individual country time series analysis but when I look at the correlation matrix, it is not showing any high figures or weird-looking ones so say it's a problem.


title: Are we correctly choosing our sample size?

text: We have a set of 60000 files, and a we have a piece of software that attempts to look at each file and determine what date the contents pertain to. We would like to be 99% confident that our software is selecting the correct date 95% +/- 5% of the time. We manually created 658 tests that each take a file and the correct date, and verify if the software selected the correct date. If 95% of our tests our passing have we met our 99% confidence mentioned above. We got the number 658 using this calculator http://www.surveysystem.com/sscalc.htm Are we correctly applying statistics to our problem? Thanks


title: What are "starting values" in glm() function?

text: What are parameters `start`, `etastart`, `mustart` in the glm() function? I have been looking in the docs and the internet but I have found no clear explanation what this means. It resembles bayesian "initial values" for the chains, but I doubt this is related, as glm() function in R is frequentist statistics...


title: What are pros and cons of Atlas.ti 7 versus Nvivo 10?

text: When I talk with colleagues about qualitative data analysis tools, Atlas.ti and Nvivo are the two names often popping out. However, few persons seem to have really tried both. After checking companies websites, one important point for me is that Nvivo will soon offer a native Mac OSX version. My university provides Web based solutions for both tools, but I would prefer to run them directly on my machine for better privacy and responsiveness.


title: Grouping similar subjects based on slopes

text: I am studying the responses of 500 subjects to temperature increase using a linear regression on each degree of temperature (from 10°C to 28°C). Thus, it was possible for me to compute the intercept and slope for each one of the 500 subjects. I am interested in identifying three groups of individuals that are "statistically identical":   * Group 1: those who are robust and respond independently of the temperature effect. These are the subjects which have a slope of approximately zero.   * Group 2: Those having negative / downward slopes    * Group 3: Those having positive / upward slopes. I have computed the slope for each subject and afterwards sorted the data from lowest to highest slope. Can you please help me to use the most relevant statistical parameters and tools that allow me to identify individuals (statistically homogeneous) of each of these three groups? I considered using the t-test, clustering techniques, PCA and discrimination but I need your help for the choice and how to do it using SAS or R. For better illustration, here is part of the data:               Sub. C15    C16     C17     C18     C19     Slope     1    30,55  30,05   29,56   29,07   28,58   -0,49     2    22,22  21,83   21,44   21,05   20,67   -0,39      3    20,16  19,78   19,39   19,01   18,63   -0,38     4    61,07  60,69   60,31   59,93   59,55   -0,38     5    49,29  48,92   48,55   48,18   47,81   -0,37     6    52,54  52,18   51,81   51,44   51,08   -0,37          238  18,19  18,18   18,18   18,18   18,18   -0,0017     239 -10,23 -10,23  -10,23  -10,24  -10,24   -0,0010     240 -14,44 -14,44  -14,44  -14,44  -14,44   -0,0006     241  19,76  19,75   19,75   19,75   19,75   -0,0006     242  13,55  13,55   13,55   13,55   13,55    0,0010     243  19,93  19,93   19,93   19,93   19,94    0,0012     244  55,69  55,69   55,69   55,69   55,69    0,0016          495 -28,70 -28,43  -28,16  -27,90  -27,63    0,27     496 -9,71  -9,40   -9,10   -8,80   -8,49     0,30     497 -12,29 -11,98  -11,67  -11,35  -11,04    0,31     498 -43,85 -43,48  -43,11  -42,74  -42,37    0,37     499 -29,41 -28,97  -28,52  -28,07  -27,62    0,45     500 -8,54  -7,98   -7,43   -6,87   -6,31     0,56      Columns `C15-C19` (for example) are predicted values of the trait of interest at each temperature. These values were the results obtained using a reaction norm model. The slope for each subject is computed here using linear regression and afterwords the data were sorted by slope values. In this example:   1. Subjects (1 to 6) have the lowest negative slopes and remain of interest in our study according to their decrease in response to increasing temperature.   2. Subjects (238-244) are here as example of subjects having a slope of approximately (zero) and are more interesting in our study because they represent the robust ones and where their performance are not affected by increase of the temperature   3. Subjects (495-500) are those having the highest positive slope. We are interested of them because their performance increase even in harsh conditions and high temperature. **Now our question is:** How can you identify more precisely (based on a statistical tool) the subjects that should correspond to each of the 3 groups? In a first step, we think we could use a t-test to identify the subjects of each group. Some people advised us to apply clustering or discriminant approaches. But we were not able to figure out how to do it the right way. We have some background in SAS. We asking how to solve our problem and help us with an example of how to solve the problem in SAS.


title: Concerns about housing data mining

text: Can anyone suggest me good papers related to housing data mining. I have googled for a while and couldn't find any recent paper. I am assuming that people have stopped doing any mining related to housing data. Isn't anyone working on it anymore? I am thinking of applying data mining techniques to predict housing prices, clustering and kind of ranking of houses based upon several factors like spatial location, crime rate information, migration rate etc etc. I am pretty sure it's not a fancy thing. There isn't much to do with it. But I was thinking of having something publishable. Can anyone suggest me good datasets that has spatial information, pricing information and other demographics. I found one at http://archive.ics.uci.edu/ml/datasets/Housing. Are there any other ones?


title: Distribution fiting on time spent of web pages

text: I'm trying to fit distributions of time spent on web pages. My goal is to cluster users based on the time spent distribution. My data set looks like following. (p1,p2 and p3 are web pages and time(in seconds) spent on those pages by each user is listed).               user p1 p2 p3     1    10 20 0      2    0   5 7     3    25  0 2      If the distribution of user '1' is equivalent to user 'n' they are assigned to same cluster. How to achieve this? I'm using python `scipy`. Original data of one user's time spent on 452 web pages.               289 0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   1177    0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   1180    0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   1687    1627    0   0   0   1469    0   0   0   0   0   1470    851 1340    0   1546    0   98  0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   1   775 0   0   815 849 0   1623    0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   1177    0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   1635    1621    0   0   0   0   0   0   0   0   0   601 1274    0   1554    0   87  0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   630 0   0   0   15  0   632 0   0   0


title: NPS confidence intervals with bootstrap

text: I am looking for a method to find confidence intervals for NPS (net promoter score). One method was already suggested here. But I am thinking of using bootstrap. Do you think the following procedure makes sense? Is it much better than what was suggested by whuber?               NPS=function(x){       if(sum(!x%%1==0)>0){stop("Non-integers found in the scores.")}       if(sum(x>10|x<0)>0){stop("Scores not on scale of 0 to 10.")}       sum(ifelse(x<7,-1,ifelse(x>8,1,0)))/length(x)*100     }          NPSconfInt=function(x,confidence=.9,iterations=10000){       quantile(bootstrap(x,iterations,NPS)$thetastar,c((1-confidence)/2, 1-(1-confidence)/2))     }               library(bootstrap)     npsData=c(1,5,6,8,9,7,0,10,7,8,6,5,7,8,2,8,10,9,8,7,10,10)     hist(npsData,breaks=11)          NPS(npsData)     NPSconfInt(npsData,.7)


title: Uniform distribution probability?

text: > Let $X$ and $Y$ be independent random variables, each of which is uniformly > distributed between $0$ and $1$. Find the probability that > $(X−1/2)^2+(Y−1/2)^2≤1/9$. Give at least $8$ correct digits after the > decimal point. I do not know where to start with this question. Can someone please get me on the right track?


title: Standard error of difference of estimates

text: I have two (non-independent) OLS-parameter estimates each with its own standard error. I'm trying to find out what the standard error of the difference of the estimates should be. Can anyone help? Is this too general?


title: What values can the process $X_n$ take when $P(Y = k) = a_k = 0$ for $k \geq 2$

text: > A newspaper uses one ton of newsprint every day. It buys its newsprint from > a local distributor. This ditributor supplies the newsprint in one-ton rolls > at the cheapest price, but unfortunatley its deliveries are erratic. Hence, > the number of rolls of newsprint in the newpaper warehouse varies randomly. > If on a particular day, there are no rolls in the warehouse, the newspaper > buys and uses a rold from an emergency supplier. Let $Y_n$ be the number of > rolls received from the local distributor on day $n$. The supply arrives in > the evening, whereas the one tone is consumed in the morning. Let $X_n$ be > the number of rolls in the warehouse at the beginning of day $n$ before the > one ton demand has to be satisfied and the supply for day $n$ has arrived. > Thus $X_{n=1} = X_n - 1 + Y_n$ if $X_n>0$ and $X_{n+1} = X_n$ is $X_n = 0$. > > Suppose that the random variables $\{Y_n, n \geq 0\}$ are idependent of > $X_{0}$. Besides, suppose that they are naturally independent and > identically distributed and $P(Y_n = k) = a_k, k = 0,1,...$. > >   * Suppose $a_k = 0$ for $k \geq 2$ and $X_0 \in \{0, 1\}$, what values can > the process $X_n$ take in this situation? > From my transition matrix, I said if we are only looking at the values for when $i, j$ are $0$ and $1$, we would get the little $2 \times 2$ matrix of the whole one step transition matrix and so this would just be $$\pmatrix{a_0 & a_1 \\ a_0 & a_1}$$ and so the only values that $X_n$ can take in this situation are $a_0$ and $a_1$. But in the answers it says it should be $0$ and $1$. Why is this?


title: Methods for evaluating predictive models in two-outcome systems

text: I'm quite an amateur data analyst/modeller trying to teach myself some new techniques. I have a system in which, on the occurrence of an event, two outcomes can occur (either 0 or 1), and a number of prediction models designed to give the probability of the outcome 1, based on previous events. Some of these models are obviously better than others, but I'm looking for a way to quantitatively assess the performance of those where it is less obvious. For example, given two models (p1 and p2) and two events (out1, out2):               p1 = 0.70     p2 = 0.60     out1 = 1          p1 = 0.65     p2 = 0.55     out2 = 0      In the first event, I would expect most to consider p1 to be better, as it has predicted a higher probability of the outcome being 1. In the second event, I would expect most to consider p2 to be better as, even though both models have predicted in the wrong direction, p2 predicted a smaller probability of the outcome being wrong. This suggests to me that simply scoring each model based on something like               if (p > 0.5 && out = 1)  score++;     if (p <= 0.5 && out = 0)  score++;      would be insufficient to accurately evaluate these models. My questions are:   * What methods could be used to evaluate these models?   * How do these methods differ in what they tell us about the models? e.g. are some methods skewed towards models which are often correct, but rarely have probabilities outside a range 0.45 < p 0.55, or skewed towards models which are less often correct but regularly have probabilities in the ranges 0 <= p < 0.2 and 0.8 < p <= 1?   * (Slightly less relevant but interesting) How do these methods translate into a system which has more than two outcomes?