title: Significance level of kwallis2 test in Stata

text: I have a question concerning the use of the `kwallis2` test in Stata. If looking at the multiple comparison output, I receive a probability and (NS) or (S).   What does the probability there mean and to what does the program compare this value to reach the outcome significant or not? I suppose that the standard 5% significance level is used. But I also see some adjusted p-value for significance... I am a little confused about the output given by this test. Another question: How can the significance level be changed? I hope (and am thankful) for good explanations!


title: Difference between Hpi.kfe, Hpi and Hscv in R's 'ks' package?

text: There are three _bandwidth selector functions_ programmed in the _ks_ package: Hpi, Hscv and Hpi.kfe, what is the difference between them? And what is the difference between a _plug-in_ and a _smoothed cross validation_ bandwidth selector?


title: Probability of outputs of binary classification in MATLAB

text: I have a binary classification problem and using neural network and SVM for it. So I choose a threshold (for instance, 0.5) for output of neural network. If output is greater than 0.5 it belongs to class 1 and if it is smaller than 0.5 it belongs to class 2. After training the network, for out sample data how can I calculate probability of outputs? For out sample data I use same criteria (0.5) and find the class out new these new data? Can we say if output of neural network is 1, the probability of belonging to class 1 (greater than 1 is class 1) is higher than for instance 0.55 ? (I used `tansig` transfer function for output layer of neural network). **Another question:** How can I find probability (or possibility) of belonging to each class in SVM? I want know find how much an out sample belong to a specific class. Can I do that? IS any function in MATLAB for it?


title: Glm gaussian Vs Glm Binomial Vs s log-linked GLM Gaussian

text: I am trying to do a study on death of malaria in certain in order to estimate the best way to predict how dangerous is this disease. I don't have a strong background in statistics, I am auto-learner building my knowledge on the online courses. First, I collected the data this way :               Statistics(gendre, age, ..)   |   Number_Death      As long as I know, GLM Binomial: help predicting if this dangerous or no. (In this case, I labled the predictor to be 0 (for no death), 1 for one or more cases). GLM Poisson : helps predicing the number of events based on the predictor. Now, I am confused : for what we use GLM Gaussian, log-linked GLM Gaussian , GLM Gamma.


title: Converting a set of tweets into vectors for clustering

text: I have a large set of tweets to which i plan to use cosine similarity to cluster the tweets. I found NLTK's GAAC to be good but how do i convert the tweets into vectors? In `nltk.cluster.gaac.demo()` the data was passed as a vector `[array([ 0.70710678, 0.70710678]), array([ 0.4472136 , 0.89442719]), array([ 0.89442719, 0.4472136 ]), array([ 1., 0.]), array([ 0.5547002 , 0.83205029]), array([ 0.9486833 , 0.31622777])]` For nltk gaac it says the input is vectors but i have plain text in the form of tweets. How do i convert my set of tweets into vectors like this.


title: testing for heteroskedasticity with Breusch Pagan in time series

text: I want to use a Breusch-Pagan test with time series data, ￼I have regressed the residual on the independent variables and added a lag for the dependent variable: is this the right way to be going about testing Breusch Pagan in time series data? ![enter image description here](http://i.stack.imgur.com/kFcyz.jpg)


title: Calculating a Constrained Mean

text: Assume that I have ten measurements of a physical quantity that must be non- negative (e.g., mass). Because of measurement error, some of the results are nevertheless negative. Can someone suggest a way to calculate the mean of the 10 measurements so as to incorporate the constraint that the true measurements and therefore the true mean must be non-negative?


title: Univariate priors for the parameters of a Beta distribution

text: I need a rather a prior on the parameters of a Beta distribution (i.e. $\alpha$ and $\beta$). I have an external constraint that requires me to use **univariate** priors, one for $\alpha$ and one $\beta$. Ideally I would like to use two univariate priors that together are as close as possible to something like $p(\alpha,\beta)∝(\alpha+\beta)^{−5/2}$ (for anyone interested in this particular choice, see this thread) What priors can I use for them?


title: Distribution of log-likelihood gradient

text: My question is simple: Is there any results regarding the distribution of log- likelihood function gradient? It may be asymptotic results as well.


title: Likelihood of dynamic conditional correlation (DCC) model goes up as model order goes down

text: The question considers the dynamic conditional correlation (DCC) model by R.Engle (2002) (full reference at the bottom). EDIT: It seems that I misinterpreted `rmgarch` output, i.e. I mixed up the sign of the reported log-likelihood. This is likely the root of the problem. (See **** below for details.) On a given dataset, I estimate two DCC models, namely, DCC(1,1) and DCC(1,2). (The estimation is done using `rmgarch` package in `R`.) The model specifications are identical except for the lag order, (1,1) versus (1,2). Thus the latter specification nests the former one. I then look at the log-likelihoods of the models. Let LLH(i,j) denote the log- likelohood of model DCC(i,j). Here comes the surprise: LLH(1,1) > LLH(1,2) Why be surprised? Because we have two nested models, and the richer model has a lower likelihood among the two. Is this some kind of mistake (software bug?) or is there a logical explanation for that? (This is not a one-off issue. I tried more specifications and another dataset, but still get similar results, i.e. a richer models have lower likelihoods than models nested in them. ****Also, I hope I interpreted `rmgarch` output correctly: the sign of log- likelihood in `rmgarch`package must be reversed. The reported log-likelihood is highly positive; this is unusual (although not impossible), thus I add a minus sign to the log-likelihood reported in `rmgarch`. Also, model estimation procedure did converge, so no problems there.) REFERENCES   Engle, Robert. "Dynamic conditional correlation: A simple class of multivariate generalized autoregressive conditional heteroskedasticity models." Journal of Business & Economic Statistics 20.3 (2002): 339-350.


title: Estimating a vector from a rank-one symmetric matrix plus scaled identity

text: I have a problem regarding estimating a $M\times 1$ vector from a given $M\times M$ symmetric matrix. The known matrix is a scaled identity matrix with a rank-one update. I have some idea how to start, but I don't know if this is the best way to go. My model of the $M\times 1$ received vector is $\boldsymbol{y} = \boldsymbol{h}s + \boldsymbol{w}$, where $\boldsymbol{h}, \boldsymbol{w} \in \mathcal{C}^M$ are independent. $s$ is a complex scalar with zero mean and unit variance, also independent of $\boldsymbol{h}$ and $\boldsymbol{w}$. With this model, the covariance of $y$ is $\boldsymbol{Q}= \boldsymbol{h}\boldsymbol{h}^\text{H} + \sigma^2\boldsymbol{I}$ (assuming variance $\sigma^2$ and independence of the elements in $\boldsymbol{w}$). I now want to estimate $\boldsymbol{h}$ from the ML estimate $\boldsymbol{S}$ of $\boldsymbol{Q}$, how do I do this? Is there a big difference if I make the assumption $\sigma^2 = 1$? My first thought was to do it the following way, using SVD: $\text{svd}(\boldsymbol{Q}) = \boldsymbol{V}(\sigma^2\boldsymbol{I} + \boldsymbol{D})\boldsymbol{V}^\text{H}$, where $\boldsymbol{V}$ is the matrix containing the eigenvectors of $\boldsymbol{h}\boldsymbol{h}^\text{H}$ and $\boldsymbol{D}$ is a diagonal matrix with one non-zero element on its diagonal. My idea was to (from the sampled version of this, i.e. $\text{svd}(\boldsymbol{S})$) take the average of the $M-1$ smallest elements on the diagonal of $\boldsymbol{D}$ to estimate $\sigma^2$ with $\hat{\sigma}^2 = \sum_{i = 2}^M d_i$, where $d_i$ is the $i$:th largest element on the diagonal of $\boldsymbol{D}$. I would then construct the matrix $R = S - \hat{\sigma}^2\boldsymbol{I}$ as my estimate of $\boldsymbol{h}\boldsymbol{h}^\text{H}$. Now, I should be able to find an estimate of $\boldsymbol{h}$, but I don't really see it. I can take the SVD of $\boldsymbol{h}\boldsymbol{h}^\text{H}$ and find the eigenvector, say $\boldsymbol{v}_1$ corresponding to the larges singular value. Now I know that $\boldsymbol{h} = c\boldsymbol{v}_1$ for some complex constant $c$. Is there any way of determining this constant with the information I have? Do you guys have any ideas or references?


title: Apply LASSO Model with nominal target in SAS?

text: I'm building a classification model with a pool of independent variables (hundreds of them). I'm in the step of variable selection/feature selection. Now I'm trying to figure out if there are any linear relationships between my Xs and Y. I want to try LASSO with SAS, but cannot make it happen with proc logistic. Any suggestions?


title: Specifying a Laplace prior using a Gaussian random variable with Gamma variance

text: I need to place a Laplace prior on a random variable, however, I want to use a Gaussian distribution whose variance is Gamma(1,1) distributed, i.e., \begin{align} x &\sim N(\mu,\sigma^2)\\\ \sigma^2&\sim\text{Gamma}(1,1) \end{align} My question is: Given this above setup, what is the corresponding form of my Laplace distribution (location and scale) for this Laplace prior?


title: Time Series Promotional Effectiveness

text: I'm trying to model the impact of two promotional tactics that ran in parallel for an year, I have the sales information at a month level and couple of metrics corresponding to the promotional tactics for the year, what can be the best approach for this situation? Should I be using ARIMA models (with appropriate order after testing for stationarity and seasonality) with the promotional metrics entered as covariates OR a GLM model with lag terms for sales and promotional metrics entered as regressors?


title: How to model when actionable data is generated on a daily basis

text: Trying to build a predictive model for attrition prediction of service desk agents using logistic regression. Data available: Daily performance metrics such as call quality,avg. call duration,attendance,hold duration etc for a duration of 6 months. Issue: all the metrics are measured on a daily basis, so how should i train the model ....should it be using the weekly avg of the scores? if yes how many weeks of data should be used to train the model. if you go through the collected data it doesn't show much variance not sure if the model would be able to predict using a data which has very little variance. Please guide.


title: Excel: Data analysis for university rankings

text: I have some university rankings data where I am asked to analyze whether a higher %of one variable leads to a higher ranking in a particular area. What kind of statistical analysis should I use in excel and how should I report the findings?


title: Using a Dictionary for Text Mining

text: I want to use a dictionary to analyze a text with text mining mechanisms. So I'm looking for a dictionary with scores for each words, so for example the word "cool", +1, gets a positive score, the word "and" a neutral score, 0, and the word "bad" a negative score, -1. Where can I find a dictionary with words and scores? Thank you in advance.


title: How can one empirically demonstrate in R which cross-validation methods the AIC and BIC are equivalent to?

text: In a question elsewhere on this site, several answers mentioned that the AIC is equivalent to leave-one-out (LOO) cross-validation and that the BIC is equivalent to K-fold cross validation. Is there a way to empirically demonstrate this in R such that the techniques involved in LOO and K-fold are made clear and demonstrated to be equivalent to the AIC and BIC values? Well commented code would be helpful in this regard. In addition, in demonstrating the BIC please use the lme4 package. See below for a sample dataset...               library(lme4) #for the BIC function          generate.data <- function(seed)     {         set.seed(seed) #Set a seed so the results are consistent (I hope)         a <- rnorm(60) #predictor         b <- rnorm(60) #predictor         c <- rnorm(60) #predictor         y <- rnorm(60)*3.5+a+b #the outcome is really a function of predictor a and b but not predictor c         data <- data.frame(y,a,b,c)          return(data)         }          data <- generate.data(76)     good.model <- lm(y ~ a+b,data=data)     bad.model <- lm(y ~ a+b+c,data=data)     AIC(good.model)     BIC(logLik(good.model))     AIC(bad.model)     BIC(logLik(bad.model))      Per earlier comments, below I have provided a list of seeds from 1 to 10000 in which AIC and BIC disagree. This was done by a simple search through the available seeds, but if someone could provide a way to generate data which would tend to produce divergent answers from these two information criteria it may be particularly informative.               notable.seeds <- read.csv("http://student.ucr.edu/~rpier001/res.csv")$seed      As an aside, I thought about ordering these seeds by the extent to which the AIC and BIC disagree which I've tried quantifying as the sum of the absolute differences of the AIC and BIC. For example,               AICDiff <- AIC(bad.model) - AIC(good.model)      BICDiff <- BIC(logLik(bad.model)) - BIC(logLik(good.model))     disagreement <- sum(abs(c(AICDiff,BICDiff)))      where my disagreement metric only reasonably applies when the observations are notable. For example,               are.diff <- sum(sign(c(AICDiff,BICDiff)))     notable <- ifelse(are.diff == 0 & AICDiff != 0,TRUE,FALSE)      However in cases where AIC and BIC disagreed, the calculated disagreement value was always the same (and is a function of sample size). Looking back at how AIC and BIC are calculated I can see why this might be the case computationally, but I'm not sure why it would be the case conceptually. If someone could elucidate that issue as well, I'd appreciate it.


title: What is the practical array implementation of a regression function on a multidimensional array and reusing it?

text: I'm going through an algorithm that does regression on some training data. My objective is to practically implement this algorithm. The training dataset $X$ consists of $n$ samples, where n = 10. Each sample $x_i$ from $X$ is an array consisting of, say, 3 features. This means that $X$ is a 2D array with 10 rows and 3 columns. Two of the lines in the algorithm say, in the midst of a loop using iterator _j_ : > Fit the regression function $g_j^h(x)$ by weighted least-squares of the > working response vector $z_{ij}$ to $x_i$ with weights $w_{ij}$ on the > training data. And then: > Set $g_j(x)$ --> $g_j(x)$ + $g_j^h(x)$ These are **Steps (2)(a)(ii) and (iii)** , on this algorithm: ![http://i.stack.imgur.com/d5C6P.png](http://i.stack.imgur.com/d5C6P.png) In this case, the index $i$ represents the sample out of the $n$ samples, and $j$ is a loop iterator. My question is - how can $g_j(x)$ be applied to a 2D array? What is the practical application of this theoretical formula in terms of arrays? With 10 samples x 3 features each, I can apply a regression function to obtain ten different sets of 3 coefficients for each of the ten samples. That is:   * I can fit the 3-feature array of **sample # 1** , in other words $x_1$ to a regression function and obtain the 3 coefficients for the model.   * I can fit the 3-feature array of $x_2$ to a regression function and obtain the 3 coefficients for the model.   * ... and so on, until I've obtained 3 coefficients for each of the ten samples. But now what? Should I even get each sample's regression coefficients? How can I use these sets of coefficients for the next step, **Step (2)(a)(iii)**? Am I correct in understanding that the implementation will go as follows:   * let's say $x_1$'s regression function has coefficients `[a b c]`   * and $x_2$'s function has coefficients `[d e f]` And so - at least for the first two samples, $g(x)$ will end up with the following coefficients: `[a+d b+e c+f]`. And likewise we'll continue adding the functions' coefficients for the next samples' cofficients, until we end up with an array of coefficients that have the following values:               - a+d+ sample 3's first coeff + ... + sample 10's first coeff     - b+e+ sample 3's second coeff + ... + sample 10's second coeff     - c+f+ sample 3's third coeff + ... + sample 10's third coeff      So that we finally have a function where we can pass in a test sample that consists of a 3-feature array...? So that in summary, the final $g(x)$ consists of coefficients that have taken into consideration each feature in each training sample (i.e. each column in each row of the training dataset). I'm not sure if this is the right forum to put this question in - but if anyone has experience in regression and in implementing algorithms, I hope they can point me in the right direction if I'm wrong! Or let me know if I'm right. I'm very confused.


title: VIF in GLM model in R

text: Before running or building a model, ho can we check on the multicollinearity between different covariates in GLM model in R? I know that SAS Proc MIXED procedure gives a column for VIF which is very easy but how to check in R?


title: Can component scores be used for further analyses, e.g. cluster analysis?

text: I have done a principal component analysis using SPSS and now have 3 components. 2 components have 4 items in the subscale, and 1 component has 3 items. Component scores using regression for each factor were added as variables to the data set by SPSS. My question is: can the component scores (the new variables that SPSS added to the data set) be used as the variables for K-means cluster analysis to form clusters? My aim is to find out the exact number of cases in each factor, which I intend to do through cluster analysis. Is this the statistically correct method? Thank you for any input.


title: predict function in glmnet

text: What kind of object the predict function in glmnet returns ? I can't access the values of predict.               x=matrix(rnorm(100*20),100,20)     g4=sample(1:4,100,replace=TRUE)     fit3=glmnet(x,g4,family="multinomial")     res <-predict(fit3,newx=x[1:2,],type="response", s=0.01)          res     , , 1                  1         2         3         4     [1,] 0.3166811 0.1829920 0.3516546 0.1486724     [2,] 0.3413679 0.2532931 0.2466071 0.1587319          res[,1] #try to get the first column      Error in res[, 1] : incorrect number of dimensions      I want to access only the first column. What is the proper method?


title: How to transition from using statistical software to understanding mathematical equations in journal articles?

text: ### Context: I'm a Psychology PhD student. As with many psychology PhD students, I know how to perform various statistical analyses using statistical software, up to techniques such as PCA, classification trees, and cluster analysis. But it's not really satisfying because though I can explain why I did an analysis and what the indicators mean, I can't explain how the technique works. The real problem is that mastering statiscal software is easy, but it is limited. To learn new techniques in articles requires that I understand how to read mathematical equations. At present I couldn't calculate eigenvalues or K-means. Equations are like a foreign language to me. ### Question:   * Is there a comprehensive guide that helps with understanding equations in journal articles? * * * ### EDIT Thanks for your comments. I thought the question would be more self explanatory: above a certain complexity, statistical notation becomes gibberish for me; let's say I would like to code my own functions in R or C++ to understand a technique but there's a barrier. I can't transform an equation into a program. And really: I don't know the situation in US doctoral schools, but in mine (France), the only courses I can follow is about some 16th century litterary movement...


title: How to analyze this design? Multiple batches of Single group pre-test post-test designs

text: first of all - thanks for reading this. I have a want to analyze the effect of an intervention, but due to practical circumstances I cannot randomly assign people to it and the interventions differ from company to company in their particular content (but after all it's leadership training for executives). So I have Single Group pre-test post-test data from multiple companies. What to do now? Is it better to do nest/pair them at first before looking for a general effect for pre vs. post test through companies? I know the single group design has major flaws, but this is just how real world for somebody doing his master's thesis is - I can't get a control group. :( If anybody would know the specific analyses that I would have to do - that would be amazing! Thank you again!


title: When would it be appropriate to report variance instead of standard deviation?

text: I have conducted an analysis in which I have modeled different variance components. When reporting the results in a table, it is much more concise to report standard deviations instead of variances. So, this brings me to the question - is there ever a reason to report variance instead of standard deviation? Is it ever more appropriate to report one over the other?


title: Extracting full set of coefficients for multilevel explanatory factor in clmm

text: I am having trouble extracting coefficients for all levels of a explanatory response from clmm() in the ordinal package. My problem seems basic, and feel I am missing something obvious, forgive me if so. I have a dataset where 120 individuals (users) where each asked to rank 8 musical clips on an ordinal scale of 1-5 (overall_rate). The 8 clips are produced by 8 different computer algorithms which edit music in different ways (edit_method). The overall question is which music editing algorithm is most preferred by listeners. Problem: I would like to examine coefficients and confidence intervals around these coefficients for the effect of each edit_method on overall_rate. However, I am only able to extract coefficients and confidence intervals for 7 levels of edit_method (see below for code). Much appreciate any help! Variables:               overall_rate #Ord.factor w/ 5 levels "0"<"1"<"2"<"3"<..: 1 2 4 2 4 2 3 3 3 5 ...     user #Factor w/ 120 levels "0","1","2","3",..: 1 1 1 1 1 1 1 1 2 2 ...     edit_method #Factor w/ 8 levels "0","1","2","3",..: 1 2 3 4 5 6 7 8 1 2 ...      Model:                   Cumulative Link Mixed Model fitted with the Laplace approximation              formula: overall_rate ~ edit_method + (1 | user)               link  threshold nobs logLik   AIC     niter      max.grad cond.H           logit flexible  960  -1330.72 2685.44 1204(4820) 8.76e-04 1.2e+02              Random effects:          Groups Name        Variance Std.Dev.          user   (Intercept) 1.328    1.152            Number of groups:  user 120               Coefficients:                      Estimate Std. Error z value Pr(>|z|)             edit_method1   1.6803     0.2466   6.814 9.46e-12 ***         edit_method2   1.5125     0.2480   6.098 1.08e-09 ***         edit_method3   0.1075     0.2404   0.447   0.6547             edit_method4   1.6177     0.2473   6.540 6.14e-11 ***         edit_method5   1.3114     0.2433   5.390 7.05e-08 ***         edit_method6   1.6869     0.2468   6.836 8.14e-12 ***         edit_method7   0.5952     0.2518   2.364   0.0181 *           ---         Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1              Threshold coefficients:             Estimate Std. Error z value         0|1  -1.6611     0.2193  -7.574         1|2   0.2035     0.2086   0.975         2|3   1.9617     0.2179   9.002          3|4   4.0311     0.2515  16.027          confint(m1)                            2.5 %     97.5 %     0|1          -2.0909141 -1.2311947     1|2          -0.2053824  0.6123330     2|3           1.5346252  2.3888445     3|4           3.5381467  4.5241025     edit_method1  1.1970306  2.1636202     edit_method2  1.0263374  1.9986432     edit_method3 -0.3636003  0.5786398     edit_method4  1.1329199  2.1024980     edit_method5  0.8345087  1.7882086     edit_method6  1.2032721  2.1705896     edit_method7  0.1016843  1.0886530      This seems strange to me. What happened to the coefficient for edit_method level 0? Thanks in advance.


title: How can I test $H_0:\sigma^2_1=\sigma^2_2$?

text: I have a few hundred estimates of a parameter calculated from two different models and I would like to know if these parameters have different variances. What is a straightforward test for comparing the variances of these parameters? (straightforward meaning, least assumptions).


title: Arima function doesn't consider seasonal components

text: Currently trying to fit several models to some data sets in order to find an accurate enough one, I ran into some difficulties with the `Arima` function of the `forecast` package. Since I am willing to use this function to fit a data set with long and multiple seasonalities, I tried using Fourier terms as regressors in the `Arima` function. This is the command lines I wrote in a script (by trying to adapt the code found here) :               n <- 40000     m <- 144     y <- ts(ledartmp[,2])     tmp <- ts(ledartmp[,3])     fourier <- function(t,terms,period)     {       n <- length(t)       X <- matrix(,nrow=n,ncol=2*terms)       for(i in 1:terms)       {         X[,2*i-1] <- sin(2*pi*i*t/period)         X[,2*i] <- cos(2*pi*i*t/period)       }       colnames(X) <- paste(c("S","C"),rep(1:terms,rep(2,terms)),sep="")       return(X)     }          library(forecast)     fit <- Arima(y[1:n], order=c(2,1,5), seasonal=c(1,2,8), xreg=cbind(fourier(1:n,4,m),fourier(1:n,4,1008), fourier(1:n,4,52560)))     plot(forecast(fit, h=14*m, xreg=cbind(fourier(n+1:(14*m),4,m),fourier(n+1:(14*m),4,1008), fourier(n+1:(14*m),4,52560))))      `m` is the shortest period. I got points every 10 minutes, and 3 types of seasonality : daily (144), weekly (1008) and yearly (52560). My problem is when I want to look at the `SARIMA` coefficients of the model `fit`, I get these results :               > fit     Series: y[1:n]      ARIMA(2,1,5)                              Coefficients: (did not bother copying them for clarity)          sigma^2 estimated as 185248:  log likelihood=-299331.2     AIC=598726.4   AICc=598726.5   BIC=599001.5     Warning message:     In sqrt(diag(x$var.coef)) : NaNs produced      Maybe I misunderstood the point of using the Fourier termes as regressors : are they just nullifying the effect of the seasonality ? Is that why the `Arima` function forget about the seasonal components ? If it is the case, how can I try different values for the `SAR and MAR` coefficients ? **Edit :** Also, is there a limit in the number of elements you can put in the `xreg` argument ? Because I have trouble adding several regressors. And how does `cbind` work ? I read the help in `?cbind` and in my case I figured the result given by :               > length(xreg)     [1] 50400      is correct since it's `2016*25` (3 4-terms Fourier series and an exogenous variable). But sometimes adding another variable to the `cbind` causes an error. **Edit 2 :** In order to try a different approach of the problem, is there any way to transpose the coefficients of the Fourier terms used as regressors (called S1, C1, S2, C...) and get the seasonal components of the `SARIMAX` model from it ? This is the results of the `Arima` function for `fit`:               > fit     Series: y[1:n]      ARIMA(2,1,5)                              Coefficients:              ar1     ar2     ma1      ma2      ma3      ma4      ma5         S1          C1          S2        C2           0.1179  0.4293  0.0309  -0.5623  -0.1999  -0.0637  -0.0516  -930.6403  -1818.6924  -1116.1193  1365.037     s.e.  0.0506  0.0380  0.0508   0.0314   0.0129   0.0078   0.0067    65.9562     56.7716     38.9532    37.514                 S3        C3        S4        C4         S1         C1        S2         C2         S3        C3           378.9935  567.3668  439.7607  666.3292  -748.6429  1897.9871  374.2509  -946.8117  -126.1115  473.9771     s.e.   32.3427   32.3188   29.3524   29.3466   321.1483   321.4964  162.8547   162.6636   110.9011  110.9124                  S4         C4         S1         C1         S2          C2         S3         C3         S4         C4           -342.0248  -336.7882  -7042.582   6872.893  -4794.241   -405.1073  -664.8625  -1983.267  -272.6207  -822.1018     s.e.    85.5908    85.5950  21015.175  23497.058  10966.602  10896.9269  7408.5802   6889.151  4555.2547  5406.5997                -0.1374     s.e.  10.5702          sigma^2 estimated as 690209:  log likelihood=-325636     AIC=651338   AICc=651338.1   BIC=651621.7      Any help would be gladly appreciated.


title: Out-of-sample vs. test set

text: Someone asked me if I did out-of-time testing (which I assume is just out-of- sample testing but with a timeline element). But if I have a test set, is that not essentially the same as out-of-sample testing? **Doing an extra out-of-sample testing would basically equal to getting an additional test set and applying my classifier to that again, wouldn't it?** Please explain if and how my interpretation is wrong. _PS: If relevant to your answer, you can assume that approximately 60% was used as training set, 20% as validation set and the last 20 as test set._


title: How is predictor importance in a cluster analysis (in SPSS) affected by dichotomy and multicollinearity?

text: I want to use a cluster analysis (CA) in SPSS to define different profiles in my dataset. I am using different continuous variables for this, including several neuropsychiatric measures. I am new in using clustering analyses, so the background is not completely familiar, but i have read several chapters and articles about it. I bumped into two problems, hopefully you can help me with these. First, when executing a Two-Step CA with continuous as well as dichotomous variables, my dichotomous variables have a high predictor importance (PI = 1.0, as seen in the model viewer from SPSS), but my continuous variables have very low predictor importance (PI < .5). In contrast, when I eliminate the dichotomous variables (and use the same distance measure), the predictor importance of my continuous variables rises. Where lies this problem? Is it a matter of small variance in the dichotomous vars? Second, a small assumption issue. In Mooi & Starstedt's book (2011) “A Concise Guide to Market Research”, a collinearity of < .90 is allowed for two variables in the same cluster analysis. Other sources (browse “collinearity cluster analysis” and you’ll find several) speak of values between .70 and .80. The point is, two of my variables in my CA (Beck Anxiety Inventory en Beck Depression Inventory) correlate with r = .715. When using both these variables in a cluster analysis, they have the largest part of the predictor importance in my model (1.00 and .84, respectively). The other variables have predictor importances of less than .2. Is this a matter of collinearity? Should I eliminate one of these two vars? Or is there another way of dealing with this in a Two-Step CA in SPSS? Thank you very much for your help.


title: Explicit optimization of RIDGE penalization parameter

text: I work with a very simple setup: I have two possibly very correlated features $x_{1,2}$, with normaliazation such as $\text{var}(x_{1,2}) = 1$, $Ey = Ex_{1,2} = 0$ (for the empirical distribution). I want to avoid multireg effects concerning the size of the $\beta$s, so I want to try the RIDGE regression. I do not have the full regression matrix but only the reduced form $\left(\sum_i x_{i,k} x_{i,l}\right)_{k,l = 1,2}$ and the rhs $\left(\sum_i y_i x_{i,k}\right)_{k = 1,2}$. Is it possible to avoid CV to find the right order of magnitude of the penalization parameter of the RIDGE regression in this simple case, let us say as a function of the correlation $\rho$? $$ \min_{\beta_{1,2}} \sum_i(y_i - \beta_1 x_{i,1} + \beta_2 x_{i,2})^2 \+ \eta( \beta_1 ^2 + \beta_2^2) $$


title: How to assess drug effects using annual death rate data?

text: I am planning a study to collect data of patients dying after brain infection versus number of cases for the next 5 years in 2 cities where they manage the disease with different drugs.   1. Will a **_year wise t-test be better or year wise Relative Risk of death and meta analysis_** be better to compare the results of the study?   2. I am a beginner in statistics. A word of explanation would benefit me. The data would look like               Year     City 1 Deaths    City 1 cases  City 2 Deaths   City 2 cases     2011     237                1124             10            1226     2012     228                1030             26           1181     2013     1500              6061             10           1122     2014     528                2320             32           1173     2015     645                3024             11           1232


title: I find very different results using a k-means or two-step clustering method. How is this?

text: I want to use a cluster analysis (CA) in SPSS to define different profiles in my dataset. I am using different continuous variables for this, including several neuropsychiatric measures. I am new in using clustering analyses, so the background is not completely familiar, but i have read several chapters and articles about it. I bumped into two problems, hopefully you can help me with these. Firstly, I try to find a cluster solution for said variables, so I tried using A) a hierarchical CA (Ward's method, Squared Euclidean distance), and following a k-means CA (k derived from the dendrogram). B) I tried using a Two-Step CA, using both the Euclidean distance, and the Log-likelihood (with both the AI Criterion and the BIC). The problem is, that I find a two cluster solution in both the hierarchical CA and the Two-Step cluster solution, but the composition of the clusters is very different (i.e., the cases aren't distributed the same). How can this be? Lies the difference in the distance measure used? And is there a CA I should prefer, with these variables? Secondly, when I perform the Two-Step CA, I get a output with a Cluster Quality measure. How does SPSS compute this measure? When I use a k-means CA, there is no option (in SPSS) for a clustering criterion like this (or like the CCC, Dunn’s validity, etc.). Is there an option, or a plug-in I can use to calculate such a criterion? Without such a criterion, my choice for k depends only on the dendrogram from the hierarchical CA. I hope you can help me with these problems.


title: What is the name of a distribution within $[a, b]$ but with a sloped straight line density?

text: What is the name of a distribution such that it looks like uniform distribution on $[a, b]$, but the density is not a horizontal line, but a line with a certain slope.


title: GLM with categorical predictor on R

text: I need to do a model with a generalized linear model. My data are these: habitat : 0 or 1, group : 1 or 0 , mortality : yes or no, and the numbers of individuals for each case (habitat, group and mortality are factors). My model is `glm(moratily~indiv+habitat+group,data,family=binomial(logit))` I want to see the influence of groupe and habitat on mortality. and I get this message error : Warning message: glm.fit: fitted probabilities numerically 0 or 1 occurred So what's wrong?


title: Several categories using k-nearest neighbour

text: Is it possible to have a training data set for the k nearest neighbour with several categories likes 700 and about 10 attributes?


title: Statistical error and error propagation

text: I have a quantity defined as: $P_{frac} = \frac{F_{max}-F_{min}}{F_{max}+F_{min}}$ I also have the value for $F_{max}$, $F_{min}$, and their statistical errors. How can I calculate the error for $P_{frac}$?


title: Why use the t-distribution for confidence interval for difference of means for unpaired samples when the population variance is unknown?

text: I am trying to understand how to derive the confidence interval for unpaired samples when the population variance is unknown and we assume that they are equal: $$({\bar X} - {\bar Y}) \pm t_{\alpha,n_x+n_y-2}S_p \sqrt{\frac{1}{n_x}+\frac{1}{n_y}} \tag{*}$$ where $S_p$ is the pooled standard deviation and the pooled variance, $S^2_p $ is: $$S^2_p=\frac{(n_x-1)S^2_x+(n_y-1)S^2_y}{n_x+n_y-2}$$ I managed to figure out the following: $$\frac{(n_x+n_y-2)S^2_p}{\sigma^2}=\frac{(n_x-1)S^2_x}{\sigma^2}+\frac{(n_y-1)S^2_y}{\sigma^2} \sim \chi^2_{n_x+n_y-2} \tag{1}$$ Now, I am supposed to use the following: $$P\left( \bar{X} - t_{\frac{a}{2},n-1} \frac{S}{\sqrt{n}} \le \mu \le \bar{X} + t_{\frac{a}{2},n-1} \frac{S}{\sqrt{n}} \right) = 100(1-\alpha)\% \tag{2}$$ $$ \frac{(\bar{X} - \bar{Y}) - (\mu_x - \mu_y)}{\sqrt{\frac{\sigma^2_x}{m}+\frac{\sigma^2_y}{n}}} \tag{3}$$ To deduce $(*)$ Please show me the steps I need to do to deduce $(*)$ as my textbook omitted it.


title: Normal Approximation of the sum of correlated Bernoulli Random Variables

text: Hi I am looking for a result (if it exists !!!) in the direction of Normal approximation for sum of correlated Bernoulli random variables (edit : with the same parameter $p$) where correlation between any pair of distinct Bernoulli variables entering into the summand is constant and equal to $0<\rho<1$. I have googled the question but with no success, if any has reference or result in this particular case I would be grateful. Edit as suggested by Glen_b here is my results for mean and variance calculations : $$E[\sum_{i=1}^{n}X_i]=n.p$$ $$E[(\sum_{i=1}^{n}X_i)^2]-E[\sum_{i=1}^{n}X_i]^2=n.p.(1-p)+ n.(n-1)\rho.p.(1-p)$$ $$=np.(1-p)(1+\rho(n-1))$$ Best regards


title: Sum of lognormal distributions

text: I am to find the expected value, variance and, preferably, the distribution of the variable $z_n$, where $z_n$ is given by \begin{equation} z_n = \exp\left\\{ \Delta t \sum_{i = 1}^n k_i \exp\left\\{a_i x \right\\} \right\\}, \end{equation} and $x$ is a random variable with Gaussian distribution and is the same for all values of $i$. Both $a_i$ and $k_i$ are known values. Is there any way to approach this problem? Are all the $x$'s independent of each other if they do not vary for different $i$'s?


title: R- Analysis of homogeneity of slopes

text: I´d like to analyse the effect of a treatment (`treatment : Factor w/ 2 levels "ambient","elevated"`) in tree diameter increment. Tree diameter is influenced by tree size. To do so, I performed the following lm:               Call:     lm(formula = BAI2013 ~ diameterJul12 * treatment, data = bandNA)          Residuals:          Min       1Q   Median       3Q      Max      -16.6493  -3.1740  -0.3767   3.3631  22.7267           Coefficients:                                      Estimate Std. Error t value Pr(>|t|)         (Intercept)                     -20.49357    2.12883  -9.627  < 2e-16 ***     diameterJul12                     1.24194    0.08876  13.992  < 2e-16 ***     treatmentelevated                10.72336    3.45783   3.101 0.002295 **      diameterJul12:treatmentelevated  -0.54953    0.14795  -3.714 0.000285 ***     ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1          Residual standard error: 6.461 on 153 degrees of freedom     Multiple R-squared:  0.6035,    Adjusted R-squared:  0.5958      F-statistic: 77.63 on 3 and 153 DF,  p-value: < 2.2e-16      How can I interpret the results? I need to figure out whether the slope of BAI~diameter is steeper for elevated trees than for ambient. Thanks


title: Linear Regression Analysis using SPSS (I think)

text: I am trying to conduct a study into whether or not female board presence affects firm performance. My background is accountancy and business NOT maths or statistics and so I don't really understand SPSS or regression analysis. However, I have been told that this is what i must do. I have obtained data on 68 companies, and am planning to use the % of Women on the Board as an independent variable and Net Income as a dependent variable. I want to compare the effect of % of women on Net Income over 68 companies over 5 years. Does anyone know how to do this in SPSS? And whether this is the right kind of analysis to be conducting based on the research parameters? Any help would be greatly appreciated,


title: Factor analysis

text: I did a principal component factor analysis with varimax rotation in SPSS. 69 variables were reduced to 12 components, but 2 components are blank and without any variable. Namely: all of my variables distributed to 12 component but the table has 14 component & 2 component is blank (without variable). Why? I guess there are 14 eigenvalues.


title: Calculate DIC from WinBUGS output

text: Has anyone ever tried to calculate DIC "manually" from WinBUGS output? I tried for the simplest linear model, and could get exact the same Dbar, but always a bit different Dhat. I think Dhat calculation is just -2(log(likehood(posterior mean of parameters))). Any clue about that? Thanks in advance!


title: Finding gaussian distribution given 15th and 50th percentile.

text: My knowledge of statistics are small so I need a little bit of help. I have some data showing the relationship between age and walking speed. The numbers of walking speed are presented only as the 15th and the 50th percentile. For example: Child(0-12 years): 1.07m/s (15th percentile) 1.33m/s (50th percentile). I need to find how many percent is between 1.0m/s - 1.2m/s and 1.2m/s - 1.4 m/s etc. That is the likelihood of observing a child walking between 1.0 and 1.2 m/s etc. Is that possible? Please help me out.


title: Time series data analysis - determining significance between treatment and control

text: We have some time series data for evaluating the effect of a drug in an animal model, with two dose levels and a control. So far we've looked at ANOVA for each time point, and there's a significant effect for the middle time points. However, the drug dose is not constant, as it needed to be reduced at later time points, which caused a convergence in the treatment and control by the end, so at the final time point the effect is not significant with ANOVA. What other approaches might we take with analyzing this data aside from the basic ANOVA? We aren't interested in prediction, just in evaluating significance at this time.


title: Transformation of multivariate normal sum of chi-squared

text: If $A$ is symmetric and $Y\sim\mathcal N(0,V)$, how can I show that $Y'AY\sim\sum_{i=1}^{t}(c_i * \chi^2(1))$ with 1 degree of freedom), where $c_i$ can be any scalar? I multiplied out the canonical case where length($Y$) is two, and got $$ \begin{aligned} Y'AY &= [Y_1 Y_2] \left[\begin{array}{cc} a_{11} & a_{12} \\\ a_{12} & a_{22} \\\ \end{array}\right] \left[\begin{array}{cc} Y_1 \\\ Y_2 \end{array}\right] \\\ &= [Y_1 Y_2] [a_{11}Y_1 + a_{12}Y_2 ; a_{12}Y_1 + a_{22}Y_2] \\\ &= a_{11} * Y_2^2 + 2a_{12}Y_1Y_2 + a_{22}Y_2^2 \end{aligned} $$ obviously the first and last term are chi-squared of degree 1, but the cross term? As $n$ increases (where the above example is $n = 2$), the cross terms increase as well.


title: Prove/counter example: A minimax decision rule is always Bayes wrt some proper prior

text: Not sure whether the claim is true or false. If claim is true, intuitively, it might have something to do with "least favorable priors", but am not able to figure out the connection. If claim is false, one example is when $X_i|\theta \sim $ Poisson$(\theta)$, then $\bar{X}$ is minimax. But a Gamma$(\alpha, \beta)$ prior fails since, that would indicate $\beta = 0$, which is improper. But, how do we know there is no other prior that gives $\bar{X}$ as a Bayes rule?


title: No-simultaneous-events assumption of the Poisson process

text: I am checking the description of Poisson process from Wikipedia. The Poisson process is based on four assumptions, but I am not clear this: > No counted occurrences are simultaneous. Can you explain this more?


title: Numerical sampling in hierarchical Bayesian models (HBM)

text: I am reading chapter 5 of $\textit{Bayesian Data Analysis}$ by Gelman $\textit{et al.}$. There it explains the few steps of data analysis for hierarchical models and if I quote from the book it will be: "We first perform the following three steps analytically.   1. Write the joint posterior density, $p( \theta, \phi|y)$, in unnormalized form as a product of the hyperprior distribution $p(\phi)$, the population distribution $p(\theta|\phi)$, and the likelihood $p(y|\theta)$.    2. Determine analytically the conditional posterior density of $\theta$ given the hyperparameters $\phi$; for fixed observed $y$, this is a function of $\phi$, $p(\theta|\phi, y)$.   3. Estimate $\phi$ using the Bayesian paradigm; that is, obtain its marginal posterior distribution, $p(\phi|y)$. The first step is immediate, and the second step is easy for conjugate models because, conditional on $t/J$, the population distribution for $\theta$ is just the iid model (5.1), so that the conditional posterior density is a product of conjugate posterior densities for the components $\theta_j$. The third step can be performed by brute force by integrating the joint posterior distribution over $\theta$: $$p(\phi|y) = \int{p(\phi,\theta|y)d\theta} \ \ \ \ \ \ \ \ \ \ \ \text{(5.4)}"$$ Now my question is whether we can (in a numerical approach) sample from $p(\phi,\theta|y)$ and ignore the $\theta$s to generate a sample of $p(\phi|y)$?


title: How to calculate p value from ANOVA function for LMM results?

text: I have used ANOVA function so that I can get the overall p value of significant factors:               > anova(lmer81)     Analysis of Variance Table                         Df Sum Sq Mean Sq F value     sex                  1   0.12   0.118  0.0195     vowel3               2 399.96 199.982 33.0859     Language             2 120.41  60.204  9.9604     sex:vowel3           2  89.73  44.865  7.4227     sex:Language         2 166.93  83.463 13.8084     vowel3:Language      4  48.27  12.067  1.9964     sex:vowel3:Language  4  52.76  13.189  2.1821      However, this does not give me p value. Can I ask how I can get p value from this?


title: Calculate rating/score based on multiple criteria to rank data

text: I extract data related to a movie by sentiment analyzing the reviews. Hence,extracted movie data contains average sentiment values (avg_pos and avg_neg) calculated over multiple reviews, total no of reviews, total positive reviews and negative reviews.(Note: total_reviews may not be equal to sum of positive and negative reviews). Now I would need to find a "score" using all these data. so that I can use the score to rank movies. I found that IMDB uses ,               weighted rating (WR) = (v ÷ (v+m)) × R + (m ÷ (v+m)) × C          where:     R = average for the movie (mean) = (avg_pos)     v = number of votes for the movie = (number of positive responses)     m = minimum votes required to be listed in the Top 250 (5)     C = the mean vote across the whole report ( sum of avg_pos/total number of movies)      But I see an issue, with it. For example, if the value "R" is same for all movies. Then, WR is also same irrespective of the number of reviews.                           V   M   R   C   WR     Christopher 1   5   0.8 0.8 0.8     Printania   38  5   0.8 0.8 0.8     Oops        32  5   0.8 0.8 0.8     bluePlanet  100 5   0.8 0.8 0.8      In-fact we need to have higher WR for "bluePlanet" because it has 100 votes with "0.8"(avg_pos) rating. Is there any algorithm, which I can use to solve this problem? Or Am I doing something wrong here?


title: Can I use the chi-squared test of independence with skewed data?

text: I have two variables, both categorical, one with skewed responses. How do you deal with skewed data in the chi-squared test? Are there any other relevant tests? I want to perform the test in SPSS.


title: Moving from conversion rates to sales volume in A/B tests

text: there is a short introduction to AB Tests in this question or here at 20bits. We are currently testing different versions of landing pages and are using the conversion rate (e.g. 4% vs. 5%) to track which version performs better. This is working fine so far. What I would like to do is start calculating whether a version performs better using the sales volume. So I could say the control version sold $\$$10.000 with 100 visitors, while the new version sold $\$$12.000 with 110 visitors. Is the difference statistically significant? I would appreciate any points in the right direction. Specifically I am having problems understanding how I can calculate 95% and 99% percentile intervals for the data above. I can share my Excel sheet for calculating "basic" A/B tests, if that helps. Thank you in advance for any help!


title: Bias towards natural numbers in the case of least squares

text: Why do we seek to minimize `x^2` instead of minimizing `|x|^1.95` or `|x|^2.05`. Are there reasons why the number should be exactly two or is it simply a convention that has the advantage of simplifying the math?


title: Difference between 'ovtest' and ovtest,rhs'

text: After running an OLS regression, I ran the Ramsey RESET test using the stata command 'ovtest' and 'ovtest,rhs'. Two tests gave the opposite results. With the first one, I could not reject the Null, while the latter rejected the Null. The stata says "rhs specifies that powers of the right-hand-side (explanatory) variables be used in the test rather than powers of the fitted values". Which test should I believe? Thank you very much


title: Standard mean and deviation

text: A company supplies specialized high-tensile pins to a customer.It uses an automatic lathe to produce the pins.Due to factors such as vibration, temperature and wear and tear, the lengths of pins are normally distributed with a mean of 25.30 mm and a standard deviation of 0.45 mm. The customer will only buy those pins with lengths in the interval 25+-0.50¨.   1. What percentage of the pins will be acceptable to the customer?   2. Suppose that the mean (25.30 mm) cannot be adjusted but the standard deviation can be reduced. Calculate the maximum reduction in the standard deviation that would make 85% of the pins acceptable.


title: showing a random variable has an exponential distribution

text: Let $X_{1},..,X_{n}$ be independent, each with a exp($\lambda$) distribution. Let $Z=min(X_{1},..X_n)$. Show that $n\lambda Z$ has an exp$(1)$ distribution. I calculate that $P(Z>z)=e^{-n\lambda z}$. Hence the density function of $n\lambda Z$ is $f(z)=(n\lambda)^2e^{-n\lambda z}$. But an $exp(1)$ distribution would have density $e^{-z}$.


title: What is a good oob score for random forests?

text: I have learning data consisting of ~45k samples, each has 21 features. I am trying to train a random forest classifier on this data, which is labelled to 3 classes (-1, 0 and 1). The classes are more or less equal in their sizes. My random forest classifier model is using `gini` as its split quality criterion, the number of trees is 10, and I have not limited the depth of a tree. Most of the features have shown negligible importance - the mean is about 5%, a third of them is of importance 0, a third of them is of importance above the mean. However, perhaps the most striking fact is the oob (out-of-bag) score: a bit less than 1%. It made me think the model fails, and indeed, testing the model on a new independent set of size ~40k, I got score of 63% (sounds good so far), but a deeper inspection of the confusion matrix have shown me that the model only succeeds for class 0, and fails in about 50% of the cases when it comes to decide between 1 and -1. Python's output attached:               array([[ 7732,   185,  6259],            [  390, 11506,   256],            [ 7442,   161,  6378]])      This is naturally because the 0 class has special properties which makes it much easier to predict. However, is it true that the oob score I've found is already a sign that the model is not good? **What is a good oob score for random forests**? Is there some law-of-thumb which helps determining whether a model is "good", using the oob score alone, or in combination with some other results of the model? * * * **Edit:** after removing bad data (about third of the data), the labels were more or less 2% for 0 and 49% for each of -1/+1. The oob score was 0.011 and the score on the test data was 0.49, with confusion matrix hardly biased towards class 1 (about 3/4 of the predictions).


title: Handling overflow warnings in pymc

text: ## Abstract I am getting numerical overflow warnings in pymc that are stalling the sampler. I'll first specify what the context is then ask more directed questions about the solution. The context portion can probably be skipped; I included it in case some feature of it is possibly relevant to the solution. ## Context I've got a mixed effects bivariate logistic AR(1) model that I am fitting to time series binary data in pymc 2.3. The model specification is as follows, where $\mathbf{l}$ and $\mathbf{t}$ are observed and $\mathcal{N}$ is parameterized by precision instead of variance. The random intercepts are $\mathbf{U}$ (subjects) and $\mathbf{v}$ (items) and the random slopes are $\mathbf{d}$. The fixed effects are encoded in $\mathbf{c}$ (conditions) and $f$ maps a subject index to its corresponding fixed effect. $$\tau_{11}, \tau_{12}, \tau_2, \lambda \sim \mathrm{Exp}(1)$$ $$u_{s1} | \tau_{11} \sim \mathcal{N}(0, \tau_{11})$$ $$u_{s2} | \tau_{12} \sim \mathcal{N}(0, \tau_{12})$$ $$d_s | \lambda \sim \mathrm{Exp}(\lambda)$$ $$v_i | \tau_2 \sim \mathcal{N}(0, \tau_2)$$ $$c_j \sim \mathcal{N}(0, 10^{-6})$$ $$l_{si1} | u_s, v_i \sim \mathrm{Bern}(\mathrm{logit}^{-1}(u_{s1} + v_i + c_{f(s)}))$$ $$l_{sij} | u_s, v_i, d_s \sim \mathrm{Bern}(\mathrm{logit}^{-1}(u_{s1} + v_i + c_{f(s)} + d_sl_{si(j-1)}))$$ $$t_{si1} | u_s, v_i \sim \mathrm{Bern}(\mathrm{logit}^{-1}(u_{s2} + v_i + c_{f(s)}))$$ $$t_{sij} | u_s, v_i, d_s \sim \mathrm{Bern}(\mathrm{logit}^{-1}(u_{s2} + v_i + c_{f(s)} + d_st_{si(j-1)}))$$ ## Problem I can run the pymc 2.3 sampler for this model at about 10 samples a second over approximately 120000 (bivariate binary) data points. The sampler will go for a few hundred thousand samples before hitting a`RuntimeWarning: overflow encountered in double_scalars`. The sampler then stalls indefinitely. I have encountered this specific warning when hand-rolling samplers in numpy/scipy. In general, it happens when the MCMC proposal is particularly bad thus making the likelihood extremely small and the negative log-likelihood extremely large. I assume this is happening here since the traceback is telling me StepMethods.py is throwing this warning. My tack in these cases is to catch the warning and generate a new proposal, since I can be nearly certain that these proposals would never be accepted. I'm surprised that pymc doesn't do something like this. (Or maybe it does?) My question: is there some parameter I should be flipping to force pymc to repropose in these cases or do I need to muck with the proposal generation code? (Or is this even a problem arising from particularly bad proposals?)


title: Power analysis for survival analysis

text: If I hypothesize that a gene signature will identify subjects at a lower risk of recurrence, that is decrease by 0.5 (hazard ratio of 0.5) the event rate in 20% of the population and I intend to use samples from a retrospective cohort study does the sample size need to be adjusted for unequal numbers in the two hypothesised groups? For example using Collett, D: Modelling Survival Data in Medical Research, Second Edition - 2nd Edition 2003. The required total number of events, d, can be found using, \begin{equation} d = \frac{(Z_{\alpha/2} + Z_{\beta/2})^2}{p_1 p_2 (\theta R)^2} \end{equation} where $Z_{\alpha/2}$ and $Z_{\beta/2}$ are the upper $\alpha/2$ and upper $\beta/2$ points, respectively, of the standard normal distribution. For the particular values,   * $p_1 = 0.20$   * $p_2 = 1 - p_1$   * $\theta R = -0.693$   * $\alpha = 0.05$ and so $Z_{0.025}= 1.96$   * $\beta = 0.10$ and so $Z_{0.05} = 1.28$,  and taking $\theta R = \log \psi R = \log 0.50 = -0.693$, the number of events required (rounded up) to have a 90% chance of detecting a hazard ratio of 0.50 to be significant at the two-sided 5% level is then given by \begin{equation} d = \frac{(1.96 + 1.28)^2}{0.20 \times 0.80\times (\log 0.5)^2}= 137 \end{equation}


title: The trinity of tests in maximum likelihood: what to do when faced with contradicting conclusions?

text: The Wald, Likelihood Ratio and Lagrange Multiplier tests in the context of maximum likelihood estimation are asymptotically equivalent. However, for small samples, they tend to diverge quite a bit, and in some cases they result in different conclusions. How can they be ranked according to how likely they are to reject the null? What to do when the tests have conflicting answers? Can you just pick the one which gives the answer you want or is there a "rule" or "guideline" as to how to proceed?


title: Variation of binary responses between repeated measures

text: I am bit confused in terms of analyzing my data. I have done an experiment which involves testing of detection systems with controlled footage.The footage has been grouped based on scene properties (or factors). I have about 7 factors and each of the factors has different levels. For example, Distance factor (see below in the data example) has 3 levels (i.e. medium, close,and far), Contrast factor has 9 levels (i.e. 0.1 to 0.2, 0.2 to 0.3) and so on. The analysis of system performance is done based on events that have been correctly detected (e.g. in the table: Event No. 1, 2, 3 and so on). The obtained results have scores of 1 to the correctly detected events and 0 to the un-detected events. To estimate the consistency of experimental results, each event was repeated 10 times (e.g. event was put through the detection system 10 times). There are some small variations on the results between the repeated times - the event was not always detected or un-detected but for some events was varied (e.g. see data table SysA_Response column, which represents the proportion/average of detection for events 1, 2, 3 and so on). I want to find out the most influential scene features (factor levels) that affect the performance of the detection systems. I guess i need to do factorial analysis for binary data but my problem is that my data are now proportions (because of the repeated testing of the same event). Any help in terms even for a bit of direction will be appreciated! Many thanks, Ana ![Example of data ](http://i.stack.imgur.com/FIgUO.png)


title: Is there an intuitive explanation why multicollinearity is a problem in linear regression?

text: The wiki discusses the problems that arise when multicollinearity is an issue in linear regression. The basic problem is multicollinearity results in unstable parameter estimates which makes it very difficult to assess the effect of independent variables on dependent variables. I understand the technical reasons behind the problems (may not be able to invert $X' X$, ill-conditioned $X' X$ etc) but I am searching for a more intuitive (perhaps geometric?) explanation for this issue. Is there a geometric or perhaps some other form of easily understandable explanation as to why multicollinearity is problematic in the context of linear regression?


title: Recognize outliers in a set as data collected

text: Suppose I am aggregating data at multiple granularities, where each key is associated with the number of counts collected during the time interval. EX: 1 minute { x:2, y:5, z:3, a:312 } This is a small example, the real set will have anywhere from 1000-6000 unique keys associated with the number of hits collected in a time interval. Is there a method that will immediately allow me to recognize that a member(s) of the set is an outlier? If not, what are some efficient algorithms to recognize outliers for data of this sort?


title: Peformance measure in the case of interval censored data

text: I am working with real world data which is interval censored data and also with the case of mixed exact observation with right censored data. I would like to know what is the formula of the Hellinger distance and Integrated square error in the case of interval censored data ? Is there any references in this case?


title: AR(1) correlation matrix

text: Suppose we have the following AR(1) correlation matrix: $$ \begin{bmatrix} 1 & \alpha & \alpha^{2} \\\ \alpha & 1 & \alpha \\\ \alpha^{2} & \alpha & 1 \end{bmatrix}$$ How many parameters need to be estimated? Wouldn't it just be $1$ since if we have an estimator for $\alpha$ then the whole matrix is determined?


title: Significance test for the difference of two correlation coefficients

text: Okay, I have two different data sets, A and B, both of which correspond to a third data set, C. I have run least-squares regression on both A with C and B with C. The correlation coefficient of A with C is 72.9%, compared to 66.5% for B with C. How do I know or how do I test that that 6.5% difference is significant?


title: Interpreting odds ratio of an ordinal regression when independent variables are negative percentages

text: I'm trying to express the results of an ordinal regression with a certain "perspective", and I'm confused. My dependent variable is an ordinal representing the progression in a scale of negative outcome (e.g. 0 = ok, 1 = bad outcome, 2 = very bad outcome). My independent variables are negative percentages that represent the percentual change in a measure between two time points (i.e. +X% if there was an increase, -Y% if the measure decreased). I have the following ORs:               Percentual change in var A: 0.895 (95% CI: 0.801; 0.988)     Percentual change in var B: 0.870 (95% CI: 0.559; 1.337)     Other variable:             1.007 (95% CI: 0.995; 1.019)     Age:                        0.970 (95% CI: 0.895; 1.045)        I am interpreting these ORs as follows: Percentual change in var B, Other variable, and age have their CIs touching 1; this means that I cannot refute the hypothesis that they have no effect on $Bad_Outcome. Percentual change in var A has a significant effect on $Bad_outcome (statistically speaking). The 0.895 OR means that an increase in var a (e.g. a change from, say, -7% to -6%) is associated with an odds 0.895 times the odds of passing from an ok outcome (0) to a bad outcome (1), or from a bad outcome (1) to a very bad outcome (2) (I am not sure if this is correct). Given that my previous interpretation is correct, how can I "translate" this results taking as reference the _decrease_ in the percentual change in var A? In plain English: I'd like to express the results as the odds of having a worse outcome (from ok to bad, and from bad to very bad) when var A _decreases_ (say, from -6% to -7%), but I'm confused by this "reversal of perspective", and I'm not sure how to "convert" the OR. I was thinking of 1/OR, but I'm really uncertain.


title: A question on discriminant analysis- Linear discriminant function

text: ![enter image description here](http://i.stack.imgur.com/Japn9.jpg) ![enter image description here](http://i.stack.imgur.com/8Zpvm.jpg) Above is part of an examination paper. I am not sure how to understand this SAS output. Especially what is there in the last table which looks to me like two discriminant functions. Can someone help me understand this? And how do I use this to classify given some observation say $X_1=94$, $X_2=491$? Thanks


title: If a random variable V is independent of two independent random variables X and Y, how to prove that V is independent of X + Y?

text: This is question 3.8.4 of An Introduction to Mathematical Statistics and Its Applications, 5th Edition, by Larsen and Marx. This is not homework for a class I am taking now, but might someday be for a class I'll take in the future. The goal could be expressed in terms of probabilities, to try and show that for every interval A and B, P(V $\in$ A and X+Y $\in$ B) = P(V $\in$ A) $\centerdot$ P(X+Y \in B). But, I've been trying to answer the question in terms of pdf's -- to show that f$_{V,X+Y}$(v,w) = f$_V$(v) $\centerdot$ f$_{X+Y}$(w), for W = X + Y, without success. The other thing I was considering was the equivalence of f$_{V+(X+Y)}$(w), f$_{(V+X)+Y)}$(w), and f$_{(V+Y)+X}$(w), for W = X + Y + V. Is this the wrong tack? Is there some easy way of proving this, say, using graph theory, with random variables as nodes and dependency relationships as edges? Thanks in advance for any hints / help.


title: Why can certain variables in a multiple regression not be included in logarithmic form?

text: 1. I have a multiple regression equation where log(salary) = b0 + b1(ceotenure). What is the purpose of putting the dependent variable in logarithmic form? How would you interpret the change in y for a 1 unit increase in ceotenure?   2. For this multiple regression equation, what are the unit measurements for RMSE?   3. I have a second multiple regression equation where salary = b0 + ln(sales) + ln(marketvalue). How does the fact that the independent variables are now in natural logarithmic form change the interpretation of a 1 unit increase in x leads to a ____ change in y?   4. I am asked to add profits to this second equation. It asks why this variable cannot be included in logarithmic form? I have no idea.


title: Distribution from $\lim{T\to\infty}\int_{-T}^{T}\phi(t) dt = 2$

text: Let X denote a real-valued random variable with distribution function F and characteristic function $\phi$. Suppose that $\phi$ satisfies the following condition: $\hspace{25mm}$$$\lim_{T\to\infty}\int_{-T}^{T}\phi(t) dt = 2$$ What can be said about the distribution? Attempt: The distribution is symmetric across the x=0 axis(?) The distribution is absolutely continuous.


title: Estimating the parameters of a normal distribution of probabilities

text: I know the question sounds weird or perhaps, silly. And my limited knowledge in statistics might have caused ambiguity in presenting the problem. But still, here's the case: Suppose there's a population of bernoulli experiments and the probability of success of these bernoulli experiments follow a normal distribution Suppose I draw, say, 20 sample randomly from those bernoulli experiments, and I run them one by one. Then I should obtain some success and some failures (which add up to 20),right? Say I do this sampling for several times, my question is, is it possible for me to use the results from these sampling to estimate the mean and variance of the distribution of probabilities of the bernoulli experiments? If yes, how can I do so? And if no, what data should I collect to do so?


title: Why does adding statistically significant interaction reduce true positives?

text: I am trying to increase my model accuracy by taking into account interaction effect of relevant variables. I am choosing variables to interact based more on the common sense than on trying every combination. So far most interaction effects with good p value ( p<.0001) and chi-square have resulted in an increase of true positives. However, the last interaction effect, while having good p values and chi square values, resulted in decreased true positives by 10.   * Should this happen?   * How do I interpret this?   * Shouldnt a variable which is significant always increase my true positives?


title: Inequality for trace of product of matrices given norms of the matrices

text: I have to show that $\sup_{n \in \mathbb{N}}|Tr(A_nB_n)| < \infty$. I know that $\sup_{n \in \mathbb{N}}||A_n||, \sup_{n \in \mathbb{N}}||B_n|| < \infty$, where $A_n,B_n$ are $2 \times 2$ random matrices and $ ||X|| = Tr(XX^T)^{1/2}$. Is there any inequality which I can use?


title: Unordered categorical variables and linear regression

text: I am new to this forum and still quite novice with regards to statistics. In my M.A. thesis I use one binary dependent variable and one unordered categorical dependent variable in separate estimations. I am therefore trying to get an oversight over the debate of whether to use linear or logistic regression when handling discrete dependent variables. The debate of whether linear regression (the Linear Probability Model) is suitable to estimate a binary outcome seems to be still ongoing. Some authors argue that the LPM is inferior to logistic regression, while others argue that it functions just as well. I therefore report results using both these models when handling the binary dependent variable. I am however struggling to find any similar debate with regards to unordered categorical variables with more than two choices. The sources Ive read simply recommend using a multinomial logistical model while referring to the fact that discrete dependent variables violate the normality and homoskedasticity assumptions. This violation is however also done by the LPM model. Is there no viable alternatives using linear regression when estimating unordered categorical variables? Any kind of feedback will be greatly appreciated!


title: Paired permutation test for repeated measures

text: I have a data set consisting of interactions between male-female dyads within a group in two conditions. The values range from -1 to 1 and indicate how responsible the female was for maintaining the relationship (Hinde's index). What I want to see is if the female responsibility changes in the two conditions. My understanding is that the Wilcoxon test is inappropriate because I have both paired and repeated measurements: each individual is part of several dyads. Female1-Male1   F1-M2   F1-M3   F2-M1....etc My other option was to run a permutation test, and I've been searching for one in R (I'm not conversant enough with R to create one myself), but so far the only ones I've found want my data to be integers, which is most definitely not the case. So I suppose my two questions are   1) Am I even on the right track here or should I be looking for a different test?   2) Can anybody suggest a permutation test that can deal with decimal places? Thanks


title: How to perform a regression on 5 years of annual observations in eviews

text: I am new to regression analysis so please excuse my ignorance. I have collected 5 years of annual panel data, and I would like to regress this data on eviews but it is saying there is an insufficient amount of observations. Is there any way around this? I currently have 3 independent variables, 1 dependent and 3 control variables. My suggestion is to just copy the same data down 4 times and regress that but I am thinking that may give me incorrect results. Any help would be appreciated even if just to confirm that it isn't possible.


title: What is the distribution of the ratio of two Poisson random variables?

text: I have a question concerning random variables. Let us assume that we have two random variables $X$ and $Y$. Let's say $X$ is Poisson distributed with parameter $\lambda_1$, and $Y$ is Poisson distributed with parameter $\lambda_2$. When you build the fracture from $X/Y$ and call this a random variable $Z$, how is this distributed and what is the mean? Is it $\lambda_1/\lambda_2$?


title: Sum up uncertainty from confidence intervals in R

text: I'm fitting a linear regression model in R:               head(data)          `x1`         `x2`     `y`           4.9711     32.70   0.2810632           6.756072   31.60   0.3115225           5.895213   56.50   0.4853171           1.951329   29.80   0.3010127            EDIT: fit <- glm(y ~ x1 + x2,family = binomial(link = logit), data = data, control= list(epsilon = 0.0001, maxit = 50, trace = F)))      where `y` is a proportion. Here's my code:               fit.l.hat=predict.glm(fit,se.fit=TRUE,type="response")         ci=c(fit.l.hat$fit - 1.96*fit.l.hat$se.fit, fit.l.hat$fit + 1.96*fit.l.hat%se.fit)      For the first fitted value, `x1 = 4.97` and `x2 = 32.7`. So from my above code, I get this output:               Fitted value: 0.35, Lower 95%: 0.01, Upper 95%: 0.68, Standard Error: 0.17      and so on. So my question is how to actually interpret this for "dummies". I would say "For an input of 4.97 and 32.7, 95 out of 100 samples will conclude a proportion between 0.01 and 0.68." Is that correct? What other information can I get out of this, or what other "dummy" ways are there to present this information? And more importantly, **how do I quantify the error?**


title: Transforming data for chi square -- squaring negative value difference scores?

text: I need to compare proportions of difference scores between 2 (unequal n) samples (males vs. females) on 2 different measures. I want to enter the difference scores into contingency table for chi-square type analysis. I subtract male-female scores for each of the 2 measures               R1C1 = Male"Yes"MeasureA - Female"Yes"MeasureA      R1C2 = Male"Yes"MeasureB - Female"Yes"MeasureB      R2C1 = Male"No"MeasureA - Female"No"MeasureA       R2C2 = Male"No"MeasureB - Female"No"MeasureB      When I use this method for getting difference scores, not surprisingly, some values are negative, which prevents me from doing chi-square. Is there a way to transform the data to do away with the negative values but preserve the proportions? For example, I was wondering if it would be acceptable to just square all difference scores, and then do chi-square? So here's a bit more detail. I am investigating a measure of sexual experiences. The original measure asks respondents to indicate whether or not they've experienced a variety of sexual encounters. The survey has 2 parallel versions -- one for females (asking about sexual victimization) and one for males (asking about perpetration). research has shown that, when given the original measure, females indicate ~2/3s increased rates of victimization, than males reported rate of perpetration. I have created a modified version of the survey (for both male and female versions) and I have hypothesized that this modified version will decrease the discrepancy rate between female/victims and male/perpetrators rates of responding. I have an unequal number of males and females. each participant was given both versions of the survey (original and modified), the original was given first. I have collapsed the response data to be dichotimous -- so either "yes" ([female]i have been raped/[male]i have raped someone) or "no" ([female]i have never been raped/[male]i have never raped anyone). So, what i need is a way to determine if the male-female discrepancy ratio on the original measure is significantly different from the male-female discrepancy ratio of the modifed version. further additional info. I have already run paired sample t-tests and determined that male report rates on the modified versin are significantly higher than on the original -- female report rates are not significantly different across versions. So i know that the discrepancy is reduced (because male reports increased and females did not) but I'm looking for a direct way to compare the difference scores/proportions between measures.


title: How to get confusion matrix with 100% precision or 100% recall in Weka

text: Here is the original confusion matrix:                   | Yes |  No  |                   -----------------        The top right corner:   False Negative      Yes|1000 | 200  |        The top left corner:    True Positive      -----------------        The bottom right corner:True Negative       No|200  | 500  |        The bottom left corner: False Positive      I expect to get these :   1. Basically, I am looking for an option or command in the classifier that can enforce WEKA not misclassifying the 200 TP instances to the FN section. In other words, the number of TP instances will be 1000 + 200 = 1200, and number of FN instances will be 0 or 100% precision.                 |  Yes  |  No  |          -------------------          Yes| 1200  |  0   |          -------------------          No | 200   | 500  |             2. Same idea for the second part, all the FP instances are moved to TN section or 100% recall.                 |  Yes  |  No  |         -------------------         Yes| 1000  | 200  |     -------------------     No |   0   | 500  |      Is there any way to implement this idea in WEKA?


title: 3D axial statistics on a surface

text: My data is a sample of 3D axes collected over a curved surface. I want to know whether these axes have preferential directions relative to the surface: i.e. whether their direction is preferentially either within the plane of the surface, or perpendicular to the plane of the surface (as opposed to a uniform spherical distribution of axes). I can think of 2 methods:   1. Apply crudely the Bingham Test of uniformity (based on the sample scatter matrix). Since my surface, although curved, happens to be rather flat, a preferential direction (or a preferential plane) should result in a large value of the Bingham statistic.   2. Compute the angle between each axis and the local vector normal to the surface.Then compare the distribution of these angles to a uniform distribution. Neither method is fully satisfactory however: (1) because it relies on the surface being rather flat; (2) because its conclusion will be very crude (deviation from a uniform distribution) and will not prove the significance of either a planar bias or a perpendicular bias.


title: Average number of failures before success

text: By the title, you must instantly think that it's a geometric probability distribution (or perhaps binomial, by the story below) that I'm looking for, but my stats is failing me. Allow me to put this as a story: An office worker is given 13 phone numbers for technical support. Each phone number will reach one and only one support agent. Each agent specializes in their own field. The office worker has no idea which agent can answer which question in advance. When they have a query, they start phoning the numbers one by one, in some random, non-repeating order. They stop after they have reached the correct number. It is guaranteed that their query will be answered by only one support agent. What is the average number of calls the worker needs to make? I.e. what is the average number of failures before he has a single success? If I try with the geometric distribution, I have p = 1/13. The number of failures before a success is k = {0,1,2,3..12}, since there is a 1/13 chance that they will get success on the first trial. The mean is defined as (1-p)/p = 12. This translate into "The average number of calls the worker is expected to make is 12". This seems excessively high. Is this correct? Thank you.


title: What is the relationship between a p-value and a confidence interval?

text: Following from this question on the difference between confidence intervals and hypothesis testing, I would love to have a simple example to better understand the the relationship between confidence intervals and p-values. For example, if I have $r=.8768 (n=300)$, then it is significant at $0.05$. The confidence interval is (lower) $0.847729 < (r) 0.8768 < 0.900619$ (higher). This is against 0 (i.e. no relationship).   * **How does it rate with the confidence level testing?**   * _What do I look for?_


title: Determining values of correction factor based on x bins in observed vs. actual data

text: I am trying to automate a problem I usually solve by hand. I have a sensor that collects data from the field. Every 6 months or so, I have to do a calibration on that sensor by collecting `x` data points. Using a calibrated device, I measure the same data measured by the field sensor and I have the following table:                   Field  | Actual      #   spd val val   err     1    89 25  28  -10.71%     2   100 32  39  -17.95%     3    54 64  58   10.34%     4    65 32  34   -5.88%     5    70 12  11    9.09%     6   100 32  37  -13.51%     7   110 48  50   -4.00%     8    78 32  32    0.00%     9    95 18  22  -18.18%     10   92 82  81    1.23%     11   78 32  30    6.67%     12   98 38  40   -5.00%     13  103 42  45   -6.67%     14   99 54  57   -5.26%     15   88 61  57    7.02%     16   79 65  62    4.84%     17   83 58  54    7.41%     18  111 39  44  -11.36%     19  112 33  38  -13.16%     20   68 28  26    7.69%      The following is the `Error versus Speed plot`: ![enter image description here](http://i.stack.imgur.com/UxXsJ.png) **How I currently solve this problem:**   * look at the data visually and determine how the data is "trending".   * I can split the data in x bins (depending on how new the sensor is, newer sensors can divide in to multiple bins)   * For the example above, I notice 3 bins:      * 50 to 70 - around 10% with one outlier     * 70 to 95 - around 3-4% with one outlier     * 95 to 120 - around -10% with one outlier.    * I then create the three bins for the sensor   * based on these bins, I determine a `fudge` or `correction` factor for each bin   * I then retest the sensor for more counts and redraw the graph. and repeat the process   * this process is time consuming especially when I have 10 sensors to look after, and that process can take up to 3 hours to complete **Question:** How can I create a statistical or mathematical solution to aid me in solving my problem? **Edit 1:** Some clarification:   * Speed: speed as the object hits the sensor   * Value: read by the sensor   * Final result is to obtain an average error of less than 6% on all reads   * Sometimes I collect more than 20 reads   * I would like to have uniformly calibrated across all ranges of speed as in the difference between actual vs measured to be relatively consistent so that every measured value could be used without significant error **Edit 2:**   * tools I have access to are Excel, Python, and R (in that order of usability for myself) **Edit 3:**   * Each object that's hitting the sensor is different in shape, size, and speed, so therefore the object hitting at any speed is not related to the size of the object. The sensor as it gets older tends to reads values differently at different speed bins.  **Edit 4:**   * Based on @Glen_b 's answer, I created the following MS Excel file for non binned results: `http://97.107.136.148/results_stats_se.xls` . Disregard the `steer` columns in the excel sheet. as it does not impact this problem. **Edit 5:**   * I collected more data as shown below.   * Sub is a partial measurement of the object, so sub is part of the full_measurement.   * There is no ratio relationship between sub and full other than it's sub Data:               Speed    Non Calibrated sensor        Calibrated sensor         Speed    Sub_1 Measurement_full    Sub_1    Measurement_full     sensor 1       88     5500      51800            5610    49010       95     5400      66000            5630    62700     sensor 2       99     4700      34600           43100    35850      103     5400      38200            5160    35360       84     4900      47700            4550    48330      101     5400      44800            5190    43050      106     5200      28400            5290    24390       96     4900      44900            4730    42900      104     5300      35000            5350    35290      102     5400      47700            5290    46810       99     5500      35500            5270    34410       89     5400      35600            5290    35490       98     5300      39300            5240    37580     sensor 3       93     5000      35600            4860    31050       79     5500      41000            5130    37810       92     4600      30700            4750    29720       99     4900      47100            4820    46230       94     5100      48300            4990    48230       93     5300      51600            5460    54850       99     5100      41700            4910    38960       85     5100      42200            4960    42810       81     5500      44100            5120    37960       98     5600      64300            5550    62510      102     5300      33400            5380    31220     sensor 4      106     4800      35000            4430    30440       99     5900      56300            5590    52740       78     5900      56300            5730    54380       96     6300      49100            5620    47790       98     5500      42100            5340    38340       70     5300      45700            4750    43670       77     5100      37900            5100    34470      106     4300      12100            3930     9750      107     5800      35600            5270    28930       90     3800      32300            3310    29390       75     5800      62200            5490    61470       88     4900      39200            5010    35720     sensor 5      104     4500      44200            5000    45120       95     5460      49400            5440    48520      104     4700      33900            5210    33420       86     5200      36400            5280    35140       88     5000      49500            5760    52250       94     5300      44400            5200    43310       95     5200      54000            5370    54460       99     5100      37000            5410    35600       91     5600      62400            5540    61860       95     5400      38000            5310    32510       94     5100      49300            4990    46210      101     4800      44300            5380    45130      100     6600      23700            7500    24630      105     5100      29200            5400    26740       94     5500      42400            5730    40760       89     5200      49600            5270    49570       99     5300      47100            5700    47970      100     4900      36600            5220    36790       90     7200      25500            7340    23480      100     5200      36500            5410    36210


title: How to transform factor scores of a PCA for a regression, in SPSS?

text: I am really a novice in statistics and really need some help! I find here, or in any books answer to my problem. If I missed it, sorry about that and could you share the link with me, please? Here, is briefly described the background of my study: I have collected a set of 20 independent variables (dummy and numerical), let's call them A, B, C...T, which are different potential reasons of getting sick. I have a population of 60 families, who during a year suffered of sickness 1, sickness 2 and sickness 3. I have three dependent variables (S1, S2, S3) telling the percentage of days that the different members of each family were sick during a year. My 20 independent variables can be explain by less factors (diet, living environment, etc.). So, I have done a Principal Component Analysis (rotated factors, varimax, etc), which gave me 5 factors. Everything works perfectly so far, even if some variables loads in different factors at the same time, their strongest loading (.6 and higher) is clearly identify in one particular variable. So I have: _Factor 1 = B, D, E, F, K, N, O and S Factor 2 = C, G and T Factor 3 = H, J and P Factor 4 = A, I and M Factor 5 = L, Q and R_ B, H, M, O, P, S and T load in more than one factor but not in a much less significant way. The total variances explained is 63%, I have good eigenvalues (respectively: 5.1, 2.1, 2.0, 1.9 and 1.6) From this PCA, I saved the factor scores for regression. (In SPSS, I used during my Facto Analysis: Scores --» save as variables --» Regression) I understand that the regression factor scores in SPSS are standardized, with a mean = 0 and Std Deviation = 1. A score of 0 on a factor therefore means that this variable's ratings of the importance of the relevant attributes is close to the average for my sample. I cannot use them directly in my regression. **Here is where I need help...** I have been told to multiply the factor loading with my original variables, and then sum them up to obtain my new variable to use in the regression. So I did the fallowing thing in spss:               *[(factor scores 1) x B] + [(factor scores 1) x D] + [(factor scores 1) x E] + [(factor   scores 1) x F] + [(factor scores 1) x K] + [(factor scores 1) x N] + [(factor scores 1) x   O] + [(factor scores 1) x S] = N-Var1* (new variable 1)          *[(factor scores 2) x C] + [(factor scores 2) x G] + [(factor scores 2) x T] = N-Var2     [(factor scores 3) x H] + [(factor scores 3) x J] + [(factor scores 3) x P] = N-Var3     [(factor scores 4) x A] + [(factor scores 4) x I] + [(factor scores 4) x M] = N-Var4     [(factor scores 5) x L] + [(factor scores 5) x Q] + [(factor scores 5) x R] = N-Var5*      Then, I used the _N-Var1, N-Var2, N-Var3, N-Var4_ and _N-Var5_ as independent variables to explain the model for the degree of sickness of S1 in a first regression, then do different regression for S2 and S3, still using the same independent variables to explain the model ( _N-Var1_ ,...). I want to measure the worst factor (= strongest influence on the dependent). The one that makes my population more easily sick with S1 (the worst) and encourage the development of S2 and S3 (I am expecting to see a strong influence on S1 from many factors). However, when my PCA is working perfectly (I'm sure about it), the results of my regression does not make sense at all (I'm also sure about it). **Does it mean I did something wrong when I transformed my factor scores? Or just that no model can work to explain S1, S2 and S3?**


title: How to analyse data from the same group where some people have opted out and data is anonymous?

text: I have a peculiar situation, I am conducting research and had a group of students participate in it. I had them scored on (property) $X$ and say $n_1$ people participated. After an intervention, I measured the same property $X$ but this time $n_2$ (where $n_2 < n_1$) people participated. Now the catch is that the experiment is so designed that the forms are filled in anonymously and there is no way to compare the scores of the same individual. So now I am left with the group mean and I want to run some tests on these. I want to know -- which tests are most suitable for this kind of data? Please explain the rationale as well or give hints /links so that I am convinced.


title: Min-max vs. standard deviation when plotting error bars

text: I am running some experiments in which I plot the values of some variable _y_ (averaged over 4 runs) against some independent variable _x_. I also compare the effect of _z_. In other words:   * horizontal axis: _x_   * vertical axis: _y_   * plots: one line for each value of _z_ I would like to know if there are any best practices or rules of thumb to decide whether my error bars should reflect (a) min-max values, or (b) standard deviations. What characteristics of the data set should lead me to prefer one over the other? Some things that I think might be relevant:   * The raw values of _y_ are integers bounded between 0 and 60.   * I can only perform a small number (e.g. 4-5) replications due to the time cost of generating each raw data point.   * There can be some pretty wild fluctuations, e.g. _y_ could be in the high 50's for 3 of the runs, but 20 in the last one.   * The data does not follow any model


title: When testing for a difference between effect sizes, is test power affected by their respective magnitude?

text: When comparing two effect sizes with a statistical test, does it matter for the test's power whether these effect sizes themselves are large or small? Imagine you are planning an experiment in which you compute an effect (difference between two conditions) for two groups (A,B). You want to show that the effect is larger in one group than the other, i.e. your _hypothesis is effect(A) > effect(B)_. Assuming that the difference between the groups is equivalent to Cohen's D = 0.5, _what is the required sample size to reach a power of 0.8_? In hypothetical experiment I the effect sizes in both groups are **1 and 1.5** , the difference between them is **0.5**. In hypothetical expeirment II the effect sizes in both groups are **0.1 and 0.6**. Again, the difference is **0.5**. Although the difference between the two groups is the same in both experiments, I would assume that experiment II requires a larger sample size. After all, the effects are located _on a much smaller scale than in experiment 1_. I have been told an a priori power analysis needs to take this into account (the magnitude of the effects that are being compared), but I do not really know how, and it's hard to find answers to this very specific problem on the web. How can I take this into account in my _a priori power analysis_? I'm using G*Power and unfortunately, I can't simply answer this by simulating data (have to use a validated method).


title: Is this a CUSUM limitation?

text: I'm using CUSUM to check structural breaks. The plot you see below pass the **CUSUM test**. How is that possible?               > cs <- sctest(pair.res ~ 1)     > cs              Recursive CUSUM test          data:  pair.res ~ 1      S = 0.9448, p-value = 0.05118      ![enter image description here](http://i.stack.imgur.com/mmZ34.png)


title: Comparing spatial patterns to on-average-distribution by averaging intervals that overlap in space

text: I am analyzing spatial patterns in vegetation structure across vegetation transects. I broke transects up into 3-m intervals that move across the transect (a moving window) 1-m at a time. So, the intervals in the analysis overlap in space. I think this is somewhat comparable to a time-series analysis as phenomena are related to previous phenomena in the time sequence. To be honest, I know very little about time series analyses. I want to compare patterns in the data to what is on average occurring on the transect. I am considering using a random-walk type approach. **My question is can I just average each interval together (so I would be averaging intervals that overlap in space) to compare to my transect data to detect spatial patterns?** Thanks in advance. To explain what kind of data I have, I measured vegetation structure (i.e. shrub heights, grass height, shrub width, etc.) along line transects. I suspect there are patches of taller shrubs and patches of shorter shrubs. To evaluate if this is the case I wanted to compare measurements along transects (broken into moving window intervals to capture variation) to a scenario where shrub heights are randomly distributed (e.g. not patterned in patches). I was thinking a way to do this would be to compare transect measurements to the average shrub height witnessed across intervals.


title: Nested design with pairing

text: I have been asked to help a colleague with analysing the following experimental design. 2 treatments (fixed), 5 time points (fixed) with 4 animals (random) at each time point. The animals within each time point are not the same. Each animal had two samples taken, one from a treated limb and one from an untreated limb so within each time point the samples from each animal are paired. I would like to analyse this data as a nested ANOVA in R. The formula I have currently is: `value ~ treatment + time/individual` Testing the measured value for the treatment at each time with individual nested within time. This doesn't account for the paired nature of the samples though. Also this is the first example of a nested design I've been askedt to analyse and I wonder if someone could advise me as to whether I'm on the right path. Thank you.


title: problem with scaling of data in PCA

text: I'd like to have your opinion about my case study. I'm analyzing a dataset of 9 experiments and 15 variables with the aim to highlight the variables that can majorly explain the variance between the experiments. This is an example with only 3 rows and 5 variables:                           var1     var2    var3    var4    var5     sample5    0.067    0.005   0.008   0.100   0.005     sample6    0.069    0.001   0.011   0.084   0.005     sample7   -7       -5      -1      34       4      My problem is that in some experiments (like in sample7) the measures related to my variables are measured as delta values (initial condition - final condition). In the other cases the variables are measured considering only the absolute values at my final condition. After PCA the model looks like strongly influenced by this difference (even if my data are centered to 0 and scaled to 1) because in the score plot I see with the first PC mainly the separation between experiments with positive and negative values and the second PC is not able to give to me any further information.   * In your opinion, is there a way to compare experiments measured in this different way?    * Alternatively do you think that the Dual Multiple Factor Analysis available with the package `FactorMineR` included in R could be a better way to analyze these data?


title: When is it acceptable to collapse across groups when performing a factor analysis?

text: Imagine that   * responses were collected on a 20 item scale which was designed to measure 4 factors with 5 items on each scale.    * participants were drawn from two groups (Group 1) and (Group 2) with sample size $n_1 = 150$ and $n_2 = 150$.   * a researcher wanted to assess the factor structure of the scale Common scenarios that I see in my consulting:   * Group 1 are first year psychology students and Group 2 is sampled from the general community   * Group 1 are sampled at one period of time and Group 2 is sampled several years later   * Group 1 is a normal population and Group 2 is a clinical population **Question**   1. Under what circumstances would it be appropriate to collapse across groups?   2. How would these circumstances be assessed? **My initial Thoughts:** My own initial thoughts were as follows:   1. **Theoretical assessment** : assess the degree to which the two groups were sampled or measured in ways that would alter the means, sds, or correlation between the items   2. **Empirical assessment:** Examine differences between means, sds, and intercorrelations on the scales and optionally on other relevant variables (e.g., demographics); perform a two-group confirmatory factor analysis to assess the consistency of the factor structure across groups. Essentially, if the empirical evidence suggests that the groups are similar and the theoretical assessment suggests that they are similar, then it should be reasonable to combine. **Conclusion**   * Does the approach above seem reasonable?   * Do you have alternative strategies?   * Are there any references that provide recommendations or examples regarding best practice in this situation?


title: How to prove that preimage of measurable space is measurable?

text: Given a function mapping between two sample spaces $S_1$ and $S_2$, if $S_2$,$F_2$ is measurable how do I show that preimage of $S_2$,$F_2$ in $S_1$ is measurable set?